text:
# TODO: add all exports
module PlasmaTools

include("constants.jl")
include("species.jl")
include("fusion.jl")
include("grids.jl")
include("collisions.jl")
include("gs2ls.jl")
#include("gs2reader.jl")

using .constants
using .species
using .fusion
using .grids
using .collisions
using .gs2ls
#using .gs2plot

export el,me,mp,kB,ep0
export SpeciesData
export Grid,createFVgrid,createAdvecDiffOperator
export slowingDownTime, nu_s, nu_par, collisionOperator, lnLambda, nuhat, gfunc
export createInputFile, GS2species, GS2resolution, GS2params, GS2geometry

end

max_stars_repo_name: gjwilkie/PlasmaTools

max_issues_repo_name: gjwilkie/PlasmaTools

max_forks_repo_name: gjwilkie/PlasmaTools


text:
const ComponentsByType = Dict{DataType, Dict{String, <:InfrastructureSystemsType}}

struct Components
    data::ComponentsByType
    time_series_storage::TimeSeriesStorage
    validation_descriptors::Vector
end

function Components(
    time_series_storage::TimeSeriesStorage,
    validation_descriptors = nothing,
)
    if isnothing(validation_descriptors)
        validation_descriptors = Vector()
    end

    return Components(ComponentsByType(), time_series_storage, validation_descriptors)
end

"""
Add a component.

Throws ArgumentError if the component's name is already stored for its concrete type.

Throws InvalidRange if any of the component's field values are outside of defined valid
range.
"""
function add_component!(
    components::Components,
    component::T;
    skip_validation = false,
) where {T <: InfrastructureSystemsType}
    component_name = get_name(component)
    if !isconcretetype(T)
        throw(ArgumentError("add_component! only accepts concrete types"))
    end

    if !haskey(components.data, T)
        components.data[T] = Dict{String, T}()
    elseif haskey(components.data[T], component_name)
        throw(ArgumentError("$(component_name) is already stored for type $T"))
    end

    if !isempty(components.validation_descriptors) && !skip_validation
        if !validate_fields(components, component)
            throw(InvalidRange("Invalid value"))
        end
    end

    if !skip_validation && !validate_struct(component)
        throw(InvalidValue("Invalid value for $(component)"))
    end

    # TODO: this check doesn't work during deserialization.
    #if has_forecasts(component)
    #    throw(ArgumentError("cannot add a component with forecasts: $component"))
    #end

    set_time_series_storage!(component, components.time_series_storage)
    components.data[T][component_name] = component
    return
end

"""
Removes all components from the system.
"""
function clear_components!(components::Components)
    for type_ in collect(keys(components.data))
        remove_components!(type_, components)
    end
end

"""
Remove all components of type T.

Throws ArgumentError if the type is not stored.
"""
function remove_components!(
    ::Type{T},
    components::Components,
) where {T <: InfrastructureSystemsType}
    if !haskey(components.data, T)
        throw(ArgumentError("component $T is not stored"))
    end

    components_ = pop!(components.data, T)
    for component in values(components_)
        prepare_for_removal!(component)
    end

    @debug "Removed all components of type" T
    return values(components_)
end

"""
Remove a component by its value.

Throws ArgumentError if the component is not stored.
"""
function remove_component!(
    components::Components,
    component::T,
) where {T <: InfrastructureSystemsType}
    return _remove_component!(T, components, get_name(component))
end

"""
Remove a component by its name.

Throws ArgumentError if the component is not stored.
"""
function remove_component!(
    ::Type{T},
    components::Components,
    name::AbstractString,
) where {T <: InfrastructureSystemsType}
    return _remove_component!(T, components, name)
end

function _remove_component!(
    ::Type{T},
    components::Components,
    name::AbstractString,
) where {T <: InfrastructureSystemsType}
    if !haskey(components.data, T)
        throw(ArgumentError("component $T is not stored"))
    end

    if !haskey(components.data[T], name)
        throw(ArgumentError("component $T name=$name is not stored"))
    end

    component = pop!(components.data[T], name)
    prepare_for_removal!(component)
    @debug "Removed component" T name
    return component
end

"""
Get the component of type T with name. Returns nothing if no component matches. If T is an abstract
type then the names of components across all subtypes of T must be unique.

See [`get_components_by_name`](@ref) for abstract types with non-unique names across subtypes.

Throws ArgumentError if T is not a concrete type and there is more than one component with
    requested name
"""
function get_component(
    ::Type{T},
    components::Components,
    name::AbstractString,
)::Union{T, Nothing} where {T <: InfrastructureSystemsType}
    if !isconcretetype(T)
        components = get_components_by_name(T, components, name)
        if length(components) > 1
            throw(ArgumentError("More than one abstract component of type $T with name $name in the system. Operation can't continue"))
        end
        return isempty(components) ? nothing : first(components)
    end

    if !haskey(components.data, T)
        @debug "components of type $T are not stored"
        return nothing
    end

    return get(components.data[T], name, nothing)
end

"""
Get the components of abstract type T with name. Note that
InfrastructureSystems enforces unique names on each concrete type but not
across concrete types.

See [`get_component`](@ref) if the concrete type is known.

Throws ArgumentError if T is not an abstract type.
"""
function get_components_by_name(
    ::Type{T},
    components::Components,
    name::AbstractString,
)::Vector{T} where {T <: InfrastructureSystemsType}
    if isconcretetype(T)
        throw(ArgumentError("get_components_by_name does not support concrete types: $T"))
    end

    components_ = Vector{T}()
    for key in keys(components.data)
        if key <: T
            component = get_component(key, components, name)
            if !isnothing(component)
                push!(components_, component)
            end
        end
    end

    return components_
end

"""
Returns an iterator of components. T can be concrete or abstract.
Call collect on the result if an array is desired.

# Arguments
- `T`: component type
- `components::Components`: Components of the sytem
- `filter_func::Union{Nothing, Function} = nothing`: Optional function that accepts a component
   of type T and returns a Bool. Apply this function to each component and only return components
   where the result is true.

See also: [`iterate_components`](@ref)
"""
function get_components(
    ::Type{T},
    components::Components,
    filter_func::Union{Nothing, Function} = nothing,
)::FlattenIteratorWrapper{T} where {T <: InfrastructureSystemsType}
    if isconcretetype(T)
        components_ = get(components.data, T, nothing)
        if !isnothing(filter_func) && !isnothing(components_)
            _filter_func = x -> filter_func(x.second)
            components_ = values(filter(_filter_func, components_))
        end
        if isnothing(components_)
            iter = FlattenIteratorWrapper(T, Vector{Base.ValueIterator}([]))
        else
            iter =
                FlattenIteratorWrapper(T, Vector{Base.ValueIterator}([values(components_)]))
        end
    else
        types = [x for x in keys(components.data) if x <: T]
        if isnothing(filter_func)
            components_ = [values(components.data[x]) for x in types]
        else
            _filter_func = x -> filter_func(x.second)
            components_ = [values(filter(_filter_func, components.data[x])) for x in types]
        end
        iter = FlattenIteratorWrapper(T, components_)
    end

    @assert eltype(iter) == T
    return iter
end

"""
Iterates over all components.

# Examples
```julia
for component in iterate_components(obj)
    @show component
end
```

See also: [`get_components`](@ref)
"""
function iterate_components(components::Components)
    Channel() do channel
        for comp_dict in values(components.data)
            for component in values(comp_dict)
                put!(channel, component)
            end
        end
    end
end

function iterate_components_with_forecasts(components::Components)
    Channel() do channel
        for comp_dict in values(components.data)
            for component in values(comp_dict)
                if has_forecasts(component)
                    put!(channel, component)
                end
            end
        end
    end
end

function JSON2.write(io::IO, components::Components)
    return JSON2.write(io, encode_for_json(components))
end

function JSON2.write(components::Components)
    return JSON2.write(encode_for_json(components))
end

function encode_for_json(components::Components)
    # Convert each name-to-value component dictionary to arrays.
    new_components = Dict{String, Vector{<:InfrastructureSystemsType}}()
    for (data_type, component_dict) in components.data
        new_components[strip_module_name(data_type)] = [x for x in values(component_dict)]
    end

    return new_components
end

"""
Return an iterable of component types deserialized from JSON.
"""
function get_component_types_raw(::Type{Components}, raw::NamedTuple)
    return propertynames(raw)
end

"""
Return an iterable of components as NamedTuples deserialized from JSON.
"""
function get_components_raw(
    ::Type{Components},
    ::Type{T},
    raw::NamedTuple,
) where {T <: InfrastructureSystemsType}
    return getproperty(raw, Symbol(strip_module_name(string(T))))
end

function get_num_components(components::Components)
    count = 0
    for components in values(components.data)
        count += length(components)
    end
    return count
end

function clear_forecasts!(components::Components)
    for component in iterate_components_with_forecasts(components)
        clear_forecasts!(component)
    end
end

function get_forecast_initial_times(components::Components)::Vector{Dates.DateTime}
    initial_times = Set{Dates.DateTime}()
    for component in iterate_components_with_forecasts(components)
        get_forecast_initial_times!(initial_times, component)
    end

    return sort!(Vector{Dates.DateTime}(collect(initial_times)))
end

function get_forecasts_initial_time(components::Components)
    initial_times = get_forecast_initial_times(components)
    if isempty(initial_times)
        throw(ArgumentError("no forecasts are stored"))
    end

    return initial_times[1]
end

function get_forecasts_last_initial_time(components::Components)
    initial_times = get_forecast_initial_times(components)
    if isempty(initial_times)
        throw(ArgumentError("no forecasts are stored"))
    end

    return initial_times[end]
end

"""
Throws DataFormatError if forecasts have inconsistent parameters.
"""
function check_forecast_consistency(components::Components)
    if !validate_forecast_consistency(components)
        throw(DataFormatError("forecasts have inconsistent parameters"))
    end
end

function validate_forecast_consistency(components::Components)
    # All component initial times must be identical.
    # We verify resolution and horizon at forecast addition.
    initial_times = nothing
    for component in iterate_components_with_forecasts(components)
        if !validate_forecast_consistency(component)
            return false
        end
        component_initial_times = Set{Dates.DateTime}()
        get_forecast_initial_times!(component_initial_times, component)
        if isnothing(initial_times)
            initial_times = component_initial_times
        elseif initial_times != component_initial_times
            @error "initial times don't match" initial_times, component_initial_times
            return false
        end
    end

    return true
end

ext: jl


text:
module CBModelTools

# include -
include("Include.jl")

# export methods -
export load_cobra_model_file
export load_vff_model_file
export write_cobra_model_file
export generate_vff_from_cobra_dictionary
export generate_cobra_dictionary_from_vff
export extract_files_from_cobra_model

# export is_generated_cobra_mat_equal

# export export_reaction_tag_to_ec_mapping_file
# export export_reaction_tag_to_gene_mapping_file
# export export_gene_order_file
# export export_reaction_mapping_file

end # module

max_stars_repo_name: varnerlab/CBModelTools

max_issues_repo_name: varnerlab/CBModelTools

max_forks_repo_name: varnerlab/CBModelTools


text:
### A Pluto.jl notebook ###
# v0.12.17

using Markdown
using InteractiveUtils

# ╔═╡ 695567b6-4307-11eb-3a59-7507eeb658d7
using Pkg;Pkg.activate()

# ╔═╡ 799dec10-4307-11eb-0fd0-6d19c3ac98e3
using Camoforge

# ╔═╡ bc46e572-430b-11eb-2bd5-3fa2afb8d877
Camoforge.get_image("/home/brad/Desktop/plains.jpg")

# ╔═╡ ed9a62ce-4307-11eb-27a0-c34b698ac0c0
Camoforge.most_common_colors("/home/brad/Desktop/plains.jpg", 10)

# ╔═╡ d329324a-430b-11eb-21b2-891efa3f11b6
Camoforge.get_image("/home/brad/Desktop/test.jpg")

# ╔═╡ b712ba46-430a-11eb-0fdc-a30ac392a610
Camoforge.most_common_colors("/home/brad/Desktop/test.jpg", 10)

# ╔═╡ Cell order:
# ╠═695567b6-4307-11eb-3a59-7507eeb658d7
# ╠═799dec10-4307-11eb-0fd0-6d19c3ac98e3
# ╠═bc46e572-430b-11eb-2bd5-3fa2afb8d877
# ╠═ed9a62ce-4307-11eb-27a0-c34b698ac0c0
# ╠═d329324a-430b-11eb-21b2-891efa3f11b6
# ╠═b712ba46-430a-11eb-0fdc-a30ac392a610


text:
using ArgParse
using ProgressMeter

include("../utils/readfiles.jl")

function parse_commandline()
    s = ArgParseSettings(description="read through d670 kernels and write to a text file with each row as: lon, lat, dep, kernel_val")
    @add_arg_table! s begin
        "--nproc"
        help = "the number of slices in the mesh"
        arg_type = Int64
        required = true
        "--base_dir"
        help = "the DATABASE directory"
        arg_type = String
        required = true
        "--output_path"
        help = "the output text file path"
        arg_type = String
        required = true
    end
    return parse_args(s)
end

function xyz2lld(x::Float64, y::Float64, z::Float64)
    R = 6371.0
    r = sqrt(x^2 + y^2 + z^2)
    dep = R - R * r
    rh = sqrt(x^2 + y^2)
    θ = atand(rh, z)
    lat = 90 - θ
    ϕ = atand(y, x)
    lon = ϕ
    return [lon,lat,dep]
end

function write2file(output_path, kernel_d670_all, position_d670_all, NSPEC2D_670::Int64, nproc::Int64)
    open(output_path, "w") do io
        for iproc in 1:nproc
            for ispec in 1:NSPEC2D_670
                for iglly in 1:NGLLY
                    for igllx in 1:NGLLX
                        kernel_val = kernel_d670_all[igllx,iglly,ispec,iproc]
                        lon, lat, dep = position_d670_all[:,igllx,iglly,ispec,iproc]
                        write(io, "$lon $lat $dep $kernel_val \n")
                    end
                end
            end
        end
end
end

function main()
    # parse command line
    command_args = parse_commandline()
    nproc = command_args["nproc"]
    base_dir = command_args["base_dir"]
    output_path = command_args["output_path"]

    # read d670 kernel and boundary mesh
    # for iproc, read the mesh file to locate the location
    boudary_disc_data = sem_boundary_disc_read(base_dir, 0)
    kernel_d670_all = zeros(Float64, NGLLX, NGLLY, boudary_disc_data.NSPEC2D_670, nproc)
    position_d670_all = zeros(Float64, 3, NGLLX, NGLLY, boudary_disc_data.NSPEC2D_670, nproc)
    @showprogress for iproc in 1:nproc
        kernel_d670_all[:,:,:,iproc] = sem_d670_read(base_dir, iproc - 1, boudary_disc_data.NSPEC2D_670)
        mesh_data_iproc = sem_mesh_data()
        mesh_data_iproc = sem_mesh_read(base_dir, iproc - 1)
        for (index_ispec, ispec) in enumerate(boudary_disc_data.ibelm_670_bot)
            for igllx in 1:NGLLX
                for iglly in 1:NGLLY
                    iglob = mesh_data_iproc.ibool[igllx,iglly,5,ispec]
                    x, y, z = mesh_data_iproc.xyz_glob[:,iglob]
                    position_d670_all[:,igllx,iglly,index_ispec,iproc] = xyz2lld(x, y, z)
                end
            end
        end
    end

    # now we write to the file
    write2file(output_path, kernel_d670_all, position_d670_all, boudary_disc_data.NSPEC2D_670, nproc)
end

main()

max_stars_repo_name: ziyixi/seisflow

max_issues_repo_name: ziyixi/seisflow

max_forks_repo_name: ziyixi/seisflow


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule MUMPS_seq_MKL_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("MUMPS_seq_MKL")
JLLWrappers.@generate_main_file("MUMPS_seq_MKL", UUID("8dbac9e0-2c8a-501c-9d6e-c02add176d6e"))
end  # module MUMPS_seq_MKL_jll


text:
module PWF

# using packages
using PowerModels

# include PWF parser file
include("pwf2dict.jl")

# include PowerModels converter files
include("pwf2pm/bus.jl")
include("pwf2pm/branch.jl")
include("pwf2pm/transformer.jl")
include("pwf2pm/load.jl")
include("pwf2pm/gen.jl")
include("pwf2pm/dcline.jl")
include("pwf2pm/shunt.jl")
include("pwf2pm/info.jl")

include("pwf2pm/correct/correct.jl")
include("pwf2pm/correct/organon.jl")
include("pwf2pm/correct/anarede.jl")

include("pwf2pm/pwf2pm.jl")

function parse_file(filename::String; validate::Bool=true, software = ANAREDE, pm::Bool = true, add_control_data::Bool=false)
    pm ? parse_pwf_to_powermodels(filename, validate = validate, software = software, add_control_data = add_control_data) : parse_pwf(filename)
end

function parse_file(io::IO; validate::Bool=true, software = ANAREDE, pm::Bool = true)
    pm ? parse_pwf_to_powermodels(io, validate = validate, software = software, add_control_data = add_control_data) : parse_pwf(filename)
end

export parse_file, ANAREDE, Organon

end # end module


text:
# Copyright (c) 2020 California Institute of Technology (“Caltech”). U.S.
# Government sponsorship acknowledged.
#
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#   • Redistributions of source code must retain the above copyright notice,
#     this list of conditions and the following disclaimer.
#   • Redistributions in binary form must reproduce the above copyright notice,
#     this list of conditions and the following disclaimer in the documentation
#     and/or other materials provided with the distribution.
#   • Neither the name of Caltech nor its operating division, the Jet
#   Propulsion Laboratory, nor the names of its contributors may be used to
#   endorse or promote products derived from this software without specific
#   prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
# THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
# PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
# OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
# WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
# OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
# ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

# Purpose: compute knl and gnl(r) for potential boundary condition

# Notes: Currently we use arbitrary precision numbers to calculate the odd
# cases where Float64 is insufficient. There are several caveats to be aware
# of.
#
# First, we save the coefficients cnl and dnl in double precision. In most
# circumstances that will be sufficient. However, there may be edge cases where
# that needs to be changed. A possible solution is to also use arbitrary
# precision to store cnl and dnl when necessary.
#
# Another solution may be to find good implementations of logbesselj() and
# logbessely() to avoid arbitrary precision.
#
# Second, knl are calculate confidently, and Float64 is sufficient to store
# them. We may use that to do entirely different approaches to calculating
# gnl(r): anything from brute-force arbitrary precision to logbessel
# implementations to an integral representation, series representation, or what
# not.
#
# Third, and perhaps more important is the question of how to test the accuracy
# of the calculated basis function gnl(r). Given confidence in knl, we would
# need to test that the gnl(r) adhere to Bessel's differential equation and
# satisfy the boundary conditions at rmin and rmax.

module GNL

export SphericalBesselGnl

using SpecialFunctions
using ArbNumerics  # Nemo conflicts currently
using Roots
using ..Splines

using Distributed
using SharedArrays

# We don't really need high precision. We need large range with floatmax() >>
# 1e300, and Arb with 1024 bits precision does that for us.
const HighprecisionFloat = ArbFloat{1024}

############## calculate zeros ##################3

function doublesecant_method(fn, a, b; maxevals=1000, xrtol=eps(1.0))
    fa = fn(a)
    fb = fn(b)
    #@show typeof(a) typeof(fa)
    #@show fa fb
    @assert sign(fa) * sign(fb) <= 0
    neval = 0
    while a != b && neval < maxevals
        m = (fb - fa) / (b - a)
        c = b - fb / m
        #if c==b || c==a
        if abs(c - b) <= xrtol * c || abs(c - a) <= xrtol * c
            return c
        end
        fc = fn(c)
        d = b - (fb + fc) / m
        #@show a,c,d,b
        #@show a-c,a-d,b-c,b-d
        #@show typeof(a) typeof(b)
        #@show typeof(c) typeof(b-d)
        #@assert a <= c <= b
        #@show fa,fc,fb
        if sign(fa) * sign(fc) <= 0
            b = c
            fb = fc
        elseif sign(fb) * sign(fc) <= 0
            a = c
            fa = fc
        else
            a = c
            fa = fc
            b = c
            fb = fc
        end

        # also try to move the other side of the interval
        if a <= d <= b
            fd = fn(d)
            if sign(fa) * sign(fd) <= 0
                b = d
                fb = fd
            elseif sign(fb) * sign(fd) <= 0
                a = d
                fa = fd
            else
                a = d
                fa = fd
                b = d
                fb = fd
            end
        end

        neval += 1
    end
    @error "maxevals reached" maxevals
    return typeof(a)(NaN)
end

function calc_next_fn_zero(func, x, δ; maxevals=1000, xrtol=1e-10)
    #@show "==="
    f = func(x)
    #@show x,f,typeof(f)
    xnew = x + δ
    fnew = func(xnew)
    while sign(fnew) == sign(f)
        #@show xnew,fnew
        x = xnew
        f = fnew
        xnew += δ
        fnew = func(xnew)
    end
    #@show x,f
    #@show xnew,fnew
    #@show xnew-2δ,func(xnew-2δ)
    #@show xnew-2δ,func(xnew-2δ)*fnew
    interval = (x, xnew)
    #@show interval func.(interval)
    #@time xnew = find_zero(func, xnew, maxevals=maxevals, xrtol=xrtol)
    #@show xnew,func(xnew)
    #xnew = find_zero(func, interval, maxevals=maxevals, xrtol=xrtol)
    #@show xnew,func(xnew)
    #@show x,xnew
    #xatol = xnew*1e-10
    #xnew = find_zero(func, interval, Roots.Brent(), maxevals=maxevals, xrtol=xrtol, xatol=xatol)
    #@show xnew,func(xnew)
    #xnew = find_zero(func, interval, Roots.FalsePosition(), maxevals=maxevals, xrtol=xrtol)
    #@show xnew,func(xnew)
    #xnew = Roots.secant_method(func, interval, maxevals=maxevals, xrtol=xrtol)
    #@show xnew,func(xnew)
    xnew = doublesecant_method(func, x, xnew, maxevals=maxevals, xrtol=xrtol)
    #@show xnew,func(xnew)
    @assert isfinite(xnew)
    return xnew
end

@doc raw"""
    calc_first_n_zeros(func, nmax; δ=π/20, xmin=0.0)

Calculate the first `nmax` zeros of the function `func`. Assumes that zeros are
spaced more than `δ` apart. First zero could be `xmin`.
"""
function calc_first_n_zeros(func, nmax; δ=π/20, xmin=0.0)
    xn = fill(NaN, nmax)
    for n=1:nmax
        xn[n] = try
            calc_next_fn_zero(func, xmin, δ)
        catch e
            isa(e, SpecialFunctions.AmosException) || rethrow(e)
            T = HighprecisionFloat
            Float64(calc_next_fn_zero(x->T(func(x)), T(xmin), T(δ)))
        end
        xmin = xn[n] + δ
    end
    return xn
end

@doc raw"""
    calc_zeros(func, xmin, xmax; δ=π/20)

Calculate all zeros of the function `func` between `xmin` and `xmax`. Assumes
that zeros are spaced more than `δ` apart.
"""
function calc_zeros(func, xmin, xmax; δ=π/20)
    xn = Float64[]
    while xmin <= xmax
        x = try
            calc_next_fn_zero(func, xmin, δ)
        catch e
            #println(e)
            isa(e, SpecialFunctions.AmosException) || rethrow(e)
            T = HighprecisionFloat
            Float64(calc_next_fn_zero(x->T(func(x)), T(xmin), T(δ)))
        end
        if xmin <= x <= xmax
            push!(xn, x)
        end
        xmin = x + δ
    end
    return xn
end

############## calculate knl for potential boundary conditions ##################

function knl_zero_function_potential(k, l, rmin, rmax)
    #@show k,l,rmin,rmax
    scale = k^3 / (1 + abs(k)^3)
    jlmax = besselj(l-1//2, k*rmax)
    (k*rmin == 0) && return jlmax
    ylmax = bessely(l-1//2, k*rmax) * scale
    jlmin = besselj(l+3//2, k*rmin)
    ylmin = bessely(l+3//2, k*rmin) * scale
    return jlmax * ylmin - jlmin * ylmax
end

@doc raw"""
    calc_knl_potential(nmax, lmax, rmin, rmax)
    calc_knl_potential(kmax, rmin, rmax)

Calculate the `knl` for potential boundary conditions.
"""
function calc_knl_potential(nmax, lmax, rmin, rmax)
    knl = fill(NaN, nmax, lmax+1)
    for l=0:lmax
        δ = π/rmax/4
        xmin = (l + 3//2) / rmax

        func0(k) = knl_zero_function_potential(k, l, rmin, rmax)
        knl[:,l+1] = calc_first_n_zeros(func0, nmax, δ=δ, xmin=xmin)
    end
    @assert all(knl .> 0)
    return knl
end

function calc_knl_potential(kmax, rmin, rmax)
    nmax = ceil(Int64, kmax * rmax / π) + 1
    lmax = ceil(Int64, kmax * rmax)
    knl = fill(NaN, nmax, lmax+1)
    for l=0:lmax
        δ = π/rmax/4
        kmin = (l + 3//2) / rmax

        func0(k) = knl_zero_function_potential(k, l, rmin, rmax)
        kn = calc_zeros(func0, kmin, kmax, δ=δ)

        if nmax < length(kn)
            @error "nmax too small" nmax length(kn) l δ kmin rmin rmax
            throw("nmax too small")
            @assert nmax >= length(kn)
        end
        (l==lmax) && @assert length(kn) == 0
        for i=1:length(kn)
            knl[i,l+1] = kn[i]
        end
    end
    nmax = findfirst(isnan.(knl[:,1])) - 1
    lmax = findfirst(isnan.(knl[1,:])) - 1 - 1
    @assert all(knl[@. isfinite(knl)] .> 0)
    return knl
end

############## calculate gnl(r) for potential boundary conditions ##################

function sphericalbesselj(nu, x::T) where {T<:ArbNumber}
    besselj_nuhalf_x = besselj(nu + 1//2, x)
    if abs(x) ≤ sqrt(eps(real(zero(besselj_nuhalf_x))))
        nu == 0 ? one(besselj_nuhalf_x) : zero(besselj_nuhalf_x)
    else
        √((float(T))(π)/2x) * besselj_nuhalf_x
    end
end
sphericalbessely(nu, x::T) where {T<:ArbNumber} = √((float(T))(π)/2x) * bessely(nu + 1//2, x)

sphericalbesselj(nu, x) = SpecialFunctions.sphericalbesselj(nu, x)
sphericalbessely(nu, x) = SpecialFunctions.sphericalbessely(nu, x)

function calc_sphbes_gnl(q, l, c::Tyl, d::Tyl) where {Tyl}
    jl = sphericalbesselj(l, Tyl(q))
    (d == 0) && return c * jl
    yl = sphericalbessely(l, Tyl(q))
    return c * jl + d * yl
end

bes_yl(l, x) = begin
    (x == 0) && return -typeof(x)(Inf)
    isnan(x) && return typeof(x)(NaN)
    return bessely(l,x)
end

function calc_cnl_dnl(knl, n, l, rmin, rmax)
    #dc1 = - besselj(l-1+1//2, knl*rmax) / bessely(l-1+1//2, knl*rmax)
    #@show dc1
    dc2 = - besselj(l+1+1//2, knl*rmin) / bes_yl(l+1+1//2, knl*rmin)
    #@show dc2
    dc = dc2 #@. √(dc1 * dc2)
    #@show dc
    gnl_rmin = calc_sphbes_gnl(knl*rmin, l, one(dc), dc)
    gnl_rmax = calc_sphbes_gnl(knl*rmax, l, one(dc), dc)
    #@show gnl_rmin gnl_rmax
    num_one = (rmax^3 * gnl_rmax^2 - rmin^3 * gnl_rmin^2) / 2
    #@show num_one
    @assert num_one >= 0
    # Note: the sign doesn't really matter
    cnl = (-1)^((1-floor(Int, 1/(l+1))) * (1-floor(Int,1/n))) / √num_one
    dnl = dc * cnl
    return cnl, dnl
end

function calc_cnl_dnl(knl, rmin, rmax, Tyl=Float64)
    nmax, lmax = size(knl) .- (0,1)
    cnl = fill(Tyl(NaN), size(knl))
    dnl = fill(Tyl(NaN), size(knl))
    n_float64 = 0
    n_arbfloat = 0
    for n=1:nmax, l=0:lmax
        !isfinite(knl[n,l+1]) && continue
        c, d = try
            n_float64 += 1
            c1, d1 = calc_cnl_dnl(knl[n,l+1], n, l, rmin, rmax)
            Tyl(c1), Tyl(d1)
        catch e
            isa(e, SpecialFunctions.AmosException) || rethrow(e)
            T = HighprecisionFloat
            n_arbfloat += 1
            c1, d1 = calc_cnl_dnl(T(knl[n,l+1]), n, l, T(rmin), T(rmax))
            Tyl(c1), Tyl(d1)
        end
        cnl[n,l+1] = c
        dnl[n,l+1] = d
    end
    @show n_float64,n_arbfloat
    return cnl, dnl
end

function gen_gnl_cache(knl, rmin, rmax, sphbesg)
    nmax, lmax = size(knl) .- (0,1)
    #T = typeof(Spline1D(0.1:0.1:1.0, Float64.(sphbesg.(1, 0, 0.1:0.1:1.0))))
    gnl = fill(Spline1D(), size(knl))
    #gnl = SharedArray{T}(size(knl)...)
    for n=1:nmax
        for l=0:lmax
            k = knl[n,l+1]
            isfinite(k) || continue
            nr = 12 * ceil(Int, 2 * k * (rmax - rmin))
            Δr = (rmax - rmin) / nr
            myr = rmin:Δr:rmax+Δr/2
            myg = sphbesg.(n,l,myr)
            gnl[n,l+1] = Spline1D(myr, Float64.(myg))
        end
    end
    return gnl
end

##################### SphericalBesselGnl ###############################33

struct SphericalBesselGnl{Tcache,Tyl}
    nmax::Int
    lmax::Int
    rmin::Float64
    rmax::Float64
    knl::Array{Float64,2}
    cnl::Array{Tyl,2}
    dnl::Array{Tyl,2}
    gnl::Tcache
end

function SphericalBesselGnl(nmax, lmax, rmin, rmax, knl, cnl, dnl)
    sphbesg = SphericalBesselGnl(nmax, lmax, rmin, rmax, knl, cnl, dnl, nothing)
    #return sphbesg  # Don't commit
    gnl = gen_gnl_cache(knl, rmin, rmax, sphbesg)
    return SphericalBesselGnl(nmax, lmax, rmin, rmax, knl, cnl, dnl, gnl)
end

@doc raw"""
    SphericalBesselGnl(nmax, lmax, rmin, rmax)
    SphericalBesselGnl(kmax, rmin, rmax)

Generate `gnl(n,l,r)`. Returns a struct that can be called for calculating
`gnl`. Note that the last argument is `r`, *not* `kr`.
"""
function SphericalBesselGnl(nmax, lmax, rmin, rmax; cache=true)
    (rmin < rmax) || @error "rmin >= rmax" rmin rmax
    knl = calc_knl_potential(nmax, lmax, rmin, rmax)
    cnl, dnl = calc_cnl_dnl(knl, rmin, rmax)
    if !cache
        return SphericalBesselGnl(nmax, lmax, rmin, rmax, knl, cnl, dnl, nothing)
    end
    return SphericalBesselGnl(nmax, lmax, rmin, rmax, knl, cnl, dnl)
end

function SphericalBesselGnl(kmax, rmin, rmax; cache=true)
    (rmin < rmax) || @error "rmin >= rmax" rmin rmax
    knl = calc_knl_potential(kmax, rmin, rmax)
    cnl, dnl = calc_cnl_dnl(knl, rmin, rmax)
    nmax, lmax = size(knl) .- (0,1)
    if !cache
        return SphericalBesselGnl(nmax, lmax, rmin, rmax, knl, cnl, dnl, nothing)
    end
    return SphericalBesselGnl(nmax, lmax, rmin, rmax, knl, cnl, dnl)
end

# evaluate gnl(r)

function evaluate(sphbesg::SphericalBesselGnl{Nothing,Tyl}, n, l, r) where{Tyl}
    k = sphbesg.knl[n,l+1]
    c = sphbesg.cnl[n,l+1]
    d = sphbesg.dnl[n,l+1]
    return calc_sphbes_gnl(k*r, l, c, d)
end

function evaluate(sphbesg::SphericalBesselGnl{Array{Tspl,2},Tyl}, n, l, r) where {Tspl,Tyl}
    return sphbesg.gnl[n,l+1](r)
end

(sphbesg::SphericalBesselGnl)(n, l, r) = evaluate(sphbesg, n, l, r)

end

# vim: set sw=4 et sts=4 :


text:
using BenchmarkTools
globalsuite = BenchmarkGroup()

include("profile_ylm.jl")

# judge(minimum(results), minimum(results1)) |> display


text:
module Drivers

# Julia modules used in the drivers
using Gridap
using Gridap.Helpers
using Gridap.Geometry
using Gridap.Arrays
using Gridap.CellData
using Gridap.MultiField: ConsecutiveMultiFieldStyle
using GridapODEs.ODETools
using GridapODEs.TransientFETools
using GridapFSI.WeakForms
using TimerOutputs
using WriteVTK
using LineSearches: BackTracking, HagerZhang
using ForwardDiff
using Test

# Julia modules extended in the drivers
import GridapODEs.TransientFETools: ∂t

export Problem
export FSIProblem
export PotentialFlowProblem

export execute
export get_problem

abstract type Problem end
struct FSIProblem{Kind} <: Problem end
struct PotentialFlowProblem{Kind} <: Problem end

function get_problem(problemName::String, kwargs)
  ptype = _get_kwarg(:ptype,kwargs,"FSI")
  if ptype == "FSI"
    return FSIProblem{Symbol(problemName)}()
  elseif ptype == "PotentialFlow"
    return PotentialFlowProblem{Symbol(problemName)}()
  else
    @notimplemented("The problem type: $ptype is not implemented")
  end
end

execute(problem::Problem; kwargs...) = @notimplemented("The driver for problem: $problem is not implemented")

function _get_kwarg(kwarg,kwargs)
  try
      return kwargs[kwarg]
  catch
      s = "The key-word argument $(kwarg) is mandatory in the $problem driver"
      error(s)
  end
end

function _get_kwarg(kwarg,kwargs,value)
  try
      return kwargs[kwarg]
  catch
      return value
  end
end

include("FSIDrivers.jl")
include("PotentialFlowDrivers.jl")

end


text:
"""
    normalize_attribute_name(name)::String

For `String` names, this simply verifies that they pass the attribute
name production, but are otherwise untouched.

For `Symbol` names, this converts `snake_case` Symbol objects to their
`kebab-case` equivalent. So that keywords, such as `for` could be used,
we strip leading underscores.
"""
function normalize_attribute_name(name::Symbol)
    name = String(name)
    if '_' in name
       if name[1] == '_'
           name = name[2:end]
       end
       name = replace(name, "_" => "-")
    end
    return normalize_attribute_name(name)
end

function normalize_attribute_name(name::String)
    # Attribute names are unquoted and do not have & escaping;
    # the &, % and \ characters don't seem to be prevented by the
    # specification, but they likely signal a programming error.
    for invalid in "/>='<&%\\\"\t\n\f\r\x20\x00"
        if invalid in name
            throw(DomainError(name, "Invalid character ('$invalid') " *
               "found within an attribute name."))
        end
    end
    if isempty(name)
        throw("Attribute name must not be empty.")
    end
    return name
end

"""
    rewrite_inside_tag(expr)

Attempt to speed up serialization of inside_tag by exploring the
expression tree at macro expansion time.
"""
function rewrite_inside_tag(expr)::Vector{Union{String, Expr}}
    if Meta.isexpr(expr, :tuple)
        args = expr.args
    elseif Meta.isexpr(expr, :call) && expr.args[1] == :Dict
        args = expr.args[2:end]
    elseif Meta.isexpr(expr, :call) && expr.args[1] == :(=>)
        args = [expr]
    else
        return [:(inside_tag($(esc(expr))))]
    end
    parts = []
    for pair in args
        if pair isa Symbol || pair isa String
            (name, value) = (pair, true)
        elseif Meta.isexpr(pair, :(=), 2)
            (name, value) = pair.args
        elseif Meta.isexpr(pair, :call, 3) && pair.args[1] == :(=>)
            (_, name, value) = pair.args
            if name isa String
                nothing
            elseif name isa QuoteNode
                name = name.value
            else
                # unexpected, use dynamic method
                return [:(inside_tag($(esc(expr))))]
            end
        else
            # unexpected, use dynamic method
            return [:(inside_tag($(esc(expr))))]
        end
        attribute = normalize_attribute_name(name)
        push!(parts, :(attribute_pair($attribute, $(esc(value)))))
    end
    return parts
end


text:
#########################################################
# Orbit Types and definitions
#########################################################

@inline Base.eltype(::AbstractOrbit{T}) where T = T

# each of them implements _Orbit Protocol_

include("OrbitT/Orbit1.jl")
include("OrbitT/Orbit2.jl")
include("OrbitT/Orbit3.jl")
include("OrbitT/Orbit4.jl")
include("OrbitT/Orbit5.jl")

const Orbit = Orbit1 # but the naive Vector+Dict is still the fastest...

for OrbT in [Symbol("Orbit", i) for i in 1:5]
    @eval begin
@doc doc"""
    Orbit(gens::Vector{<:GrouElem}, pt[, ^])
Compute the orbit of a point `pt` under the action of group generated by `gens`.

Only the consecutive points are stored in the orbit.
"""
        function $OrbT(gens::AbstractVector{<:GroupElem}, pt::T, op=^) where T
            orb = $OrbT(pt)
            for o in orb
                for g in gens
                    γ = op(o, g)
                    if γ ∉ orb
                        push!(orb,γ)
                    end
                end
            end
            return orb
        end

@doc doc"""
    Orbit(::Type{Transversal}, gens::Vector{<:GrouElem}, pt[, ^])
Compute the orbit of a point `pt` under the action of group generated by `gens`.

Along the points the right transversal is stored as explicit group elements.
"""
        function $OrbT(::Type{Transversal}, gens::Vector{<:GroupElem}, pt::T, op=^) where T
            orb = $OrbT(pt, one(parent(first(gens))))
            for o in orb
                for g in gens
                    γ = op(o, g)
                    if γ ∉ orb
                        push!(orb,(γ, orb[o]*g))
                    end
                end
            end
            return orb
        end
    end
end


text:
function import_has_dots(x)
    i = 2
    while CSTParser.is_dot(x.args[i])
        i +=1
    end
    return i - 2
end

function unwrap_imports(x)
    !all(CSTParser.is_colon(a) || a isa CSTParser.PUNCTUATION || a isa CSTParser.IDENTIFIER || a isa CSTParser.KEYWORD for a in x) && return []
    prefix = Symbol[]
    vars = Vector{Symbol}[]

    i = 2
    if any(CSTParser.is_colon(a) for a in x)
        while !CSTParser.is_colon(x.args[i])
            !(x.args[i] isa CSTParser.PUNCTUATION) && push!(prefix, Symbol(x.args[i].val))
            i += 1
        end
        i += 1
    end

    while i <= length(x.args)
        if CSTParser.is_comma(x.args[i])
            i += 1
        end
        var = copy(prefix)
        while i<= length(x.args) && !CSTParser.is_comma(x.args[i])
            !(x.args[i] isa CSTParser.PUNCTUATION) && push!(var, Symbol(x.args[i].val))
            i += 1
        end
        push!(vars, var)
    end
    return vars
end

function is_pkg_available(pkg::Symbol, S)
    string(pkg) in readdir(Pkg.dir())
end

function is_pkg_loaded(pkg::Symbol, S)
    string(pkg) in keys(loaded_mods)
end

function get_imports(x, S) end
function get_imports(x::CSTParser.EXPR{T}, S) where T <: Union{CSTParser.Using,CSTParser.Import,CSTParser.ImportAll}
    u = T == CSTParser.Using
    vars = unwrap_imports(x)
    for v in vars
        if u && string(v[1]) in keys(S.current_scope.names) && S.current_scope.names[string(v[1])][end].t == :Module
            mx = S.current_scope.names[string(v[1])][end].val.args[3].args
            for a in mx
                if a isa CSTParser.EXPR{CSTParser.Export}
                    for i = 2:length(a.args)
                        if a.args[i] isa CSTParser.IDENTIFIER
                            add_binding(x, CSTParser.str_value(a.args[i]), :Any, S::State, S.loc.offset + x.span)
                        end
                    end
                end
            end
        else
            add_binding(x, string(v[end]), :Any, S::State, S.loc.offset + x.span)
        end
    end
end

function mod_names(m::Module, d = Dict{String,Tuple{Set{Symbol},Set{Symbol}}}())
    ext = names(m)
    int = names(m, true, true)
    d[string(m)] = (Set(ext), Set(int))
    for n in int
        if isdefined(m, n) && getfield(m, n) isa Module && !(string(getfield(m, n)) in keys(d))
            mod_names(getfield(m, n), d)
        end
    end
    d
end


text:
include("units.jl")
include("performance.jl")


text:
# bounds.jl

"""
   bound_remainder(::Type{TaylorModel1}f::Function, polf::Taylor1, polfI::Taylor1, x0::Interval, I::Interval)

Bound the absolute remainder of the polynomial approximation of `f` given
by the Taylor polynomial `polf` around `x0` on the interval `I`. It requires
the interval extension `polfI` of the polynomial that approximates `f` for
the whole interval `I`, in order to compute the Lagrange remainder.

If `polfI[end]` has a definite sign, then it is monotonic in the intervals
[I.lo, x0] and [x0.hi, I.hi], which is exploited; otherwise, it is used
to compute the Lagrange remainder. This corresponds to Prop 2.2.1 in Mioara
Joldes PhD thesis (pp 52).

"""
function bound_remainder(::Type{TaylorModel1}, f::Function, polf::Taylor1, polfI::Taylor1, x0, I::Interval)

    _order = get_order(polf) + 1
    fTIend = polfI[_order]
    if (sup(fTIend) < 0 || inf(fTIend) > 0)
        # Absolute remainder is monotonic
        a = Interval(inf(I))
        b = Interval(sup(I))
        Δlo = f(a) - polf(a-x0)
        # Δlo = f(a) - bound_taylor1(polf, a-x0)
        Δhi = f(b) - polf(b-x0)
        # Δhi = f(b) - bound_taylor1(polf, b-x0)
        Δx0 = f(x0) - polf[0]
        Δ = hull(Δlo, Δx0, Δhi)
    else
        # Lagrange bound
        Δ = fTIend * (I-x0)^_order
    end
    return Δ
end
function bound_remainder(::Type{TaylorModel1}, f::Function, polf::Taylor1{TaylorN{T}}, polfI::Taylor1, x0, I::Interval) where {T}
    # The domain of the TaylorN part is assumed to be
    # the normalized symmetric box
    _order = get_order(polf) + 1
    fTIend = polfI[_order]
    if (sup(fTIend) < 0 || inf(fTIend) > 0)
        # Absolute remainder is monotonic
        a = Interval(inf(I))
        b = Interval(sup(I))
        symIbox = IntervalBox(-1 .. 1, get_numvars())
        Δlo = (f(a) - polf(a-x0))(symIbox)
        # Δlo = f(a) - bound_taylor1(polf, a-x0)
        Δhi = (f(b) - polf(b-x0))(symIbox)
        # Δhi = f(b) - bound_taylor1(polf, b-x0)
        Δx0 = (f(x0) - polf[0])(symIbox)
        Δ = hull(Δlo, Δx0, Δhi)
    else
        # Lagrange bound
        Δ = fTIend * (I-x0)^_order
    end
    return Δ
end

"""
   bound_remainder(::Type{RTaylorModel1}, f::Function, polf::Taylor1, polfI::Taylor1, x0::Interval, I::Interval)

Bound the relative remainder of the polynomial approximation of `f` given
by the Taylor polynomial `polf` around `x0` on the interval `I`. It requires
an the interval extension `polfI` of a polynomial that approximates `f` for
the whole interval `I`, in order to compute the Lagrange remainder.

If `polfI[end]` has a definite sign, then it is monotonic in the interval `I`,
which is exploited; otherwise, the last coefficients bounds the relative
remainder. This corresponds to Prop 2.3.7 in Mioara Joldes' PhD thesis (pp 67).

"""
function bound_remainder(::Type{RTaylorModel1}, f::Function, polf::Taylor1, polfI::Taylor1, x0, I::Interval)

    _order = get_order(polf) + 1
    fTIend = polfI[_order+1]
    a = Interval(inf(I))
    b = Interval(sup(I))
    if (sup(fTIend) < 0 || inf(fTIend) > 0) && isempty(a ∩ x0) && isempty(b ∩ x0)
        # Error is monotonic
        denom_lo = (a-x0)^_order
        Δlo = f(a) - polf(a-x0)
        # Δlo = f(a) - bound_taylor1(polf, a-x0)
        Δlo = Δlo / denom_lo
        denom_hi = (b-x0)^_order
        Δhi = f(b) - polf(b-x0)
        # Δhi = f(b) - bound_taylor1(polf, b-x0)
        Δhi = Δhi / denom_hi
        Δ = hull(Δlo, Δhi)
    else
        # Lagrange coefficient
        Δ = polfI[_order]
    end
    return Δ
end

"""
    bound_taylor1(fT::Taylor1, I::Interval)

Compute a *tight* polynomial bound for the Taylor polynomial `fT`
in the interval `I`.

Note: Algorithm 2.1.1 corresponds to `evaluate(fT, I)` or simply `fT(I).
This function uses the roots of the derivative of `ft` to obtain a
tighter bound.

"""
function bound_taylor1(fT::Taylor1, I::Interval)

    # Check if the fT is monotonous (the derivative has a given sign)
    fTd  = TaylorSeries.derivative(fT)
    range_deriv = fTd(I)
    (sup(range_deriv) ≤ 0 || inf(range_deriv) ≥ 0) && return bound_taylor1(fT, fTd, I)

    # Compute roots of the derivative using the second derivative
    # Fix some sort of relative tolerance for Newton root search
    fTd2 = TaylorSeries.derivative(fTd)
    rootsder = roots(x->fTd(x), x->fTd2(x), I, Newton, 1.0e-5*mag(I))

    # Bound the range of fT using the roots and end points
    num_roots = length(rootsder)
    num_roots == 0 && return fT(I)
    rangepoly = hull( fT(Interval(inf(I))), fT(Interval(sup(I))) )
    @inbounds for ind in 1:num_roots
        rangepoly = hull(rangepoly, fT(rootsder[ind].interval))
    end

    return rangepoly
end

"""
    bound_taylor1(fT::Taylor1, fTd::Taylor1, I::Interval)

Compute a *tight* polynomial bound for the Taylor polynomial `fT`
in the interval `I`, considering whether its derivative `ftd` has
a definite sign.

"""
function bound_taylor1(fT::Taylor1{T}, fTd::Taylor1{T}, I::Interval{T}) where {T}
    #
    I_lo = inf(I)
    I_hi = sup(I)
    if inf(fTd(I)) ≥ 0
        return Interval(fT(I_lo), fT(I_hi))
    elseif sup(fTd(I)) ≤ 0
        return Interval(fT(I_hi), fT(I_lo))
    end
    return fT(I)
end
function bound_taylor1(fT::Taylor1{Interval{T}}, fTd::Taylor1{Interval{T}},
        I::Interval{S}) where {T,S}
    #
    I_lo = inf(I)
    I_hi = sup(I)
    if inf(fTd(I)) ≥ 0
        return hull(fT(I_lo), fT(I_hi))
    elseif sup(fTd(I)) ≤ 0
        return hull(fT(I_hi), fT(I_lo))
    end
    return fT(I)
end

"""
    bound_taylor1(fT::TaylorModel1, I=domain(fT)::Interval)

Compute a *tight* polynomial bound for the Taylor model `fT`
in the interval `I`, considering whether its derivative `ftd` has
a definite sign.

"""
bound_taylor1(fT::TaylorModel1, I=domain(fT)::Interval) = bound_taylor1(polynomial(fT), I)

"""
    linear_dominated_bounder(fT::TaylorModel1, ϵ=1e-3::Float64, max_iter=5::Int)

Compute a tighter polynomial bound for the Taylor model `fT` by the linear
dominated bounder algorithm. The linear dominated algorithm is applied until
the bound of `fT` gets tighter than `ϵ` or the number of steps reachs `max_iter`.
The returned bound corresponds to the improved polynomial bound with the remainder
of the `TaylorModel1` included.
"""
function linear_dominated_bounder(fT::TaylorModel1{T, S}; ϵ=1e-3, max_iter=5) where {T, S}
    d = one(T)
    dom = domain(fT)
    x0 = expansion_point(fT)
    pol = polynomial(fT)
    Pm = deepcopy(pol)
    bound = zero_interval(S)
    hi = sup(pol(dom - x0))

    n_iter = 0
    while d > ϵ && n_iter < max_iter
        x0 = mid(dom - x0)
        c = mid(dom)
        update!(Pm, x0)
        non_linear = nonlinear_polynomial(Pm)
        linear = Pm - non_linear
        if T <: Interval
            Li = mid(linear[1])
        else
            Li = linear[1]
        end
        I1 = linear(dom - c)
        Ih = non_linear(dom - c)
        bound = inf(I1) + Ih
        d = diam(bound)
        n_iter += 1
        dom_lo = inf(dom)
        dom_hi = sup(dom)
        if Li == 0
            break
        elseif Li > 0
            new_hi = min(dom_lo + (d / abs(Li)), dom_hi)
            x0 = dom
            dom = Interval(dom_lo, new_hi)
        else
            new_lo = max(dom_hi - (d / abs(Li)), dom_lo)
            x0 = dom
            dom = Interval(new_lo, dom_hi)
        end
    end

    return Interval(inf(bound), hi) + remainder(fT)
end

"""
    linear_dominated_bounder(fT::TaylorModelN, ϵ=1e-3::Float64, max_iter=5::Int)

Compute a tighter polynomial bound for the Taylor model `fT` by the linear
dominated bounder algorithm. The linear dominated algorithm is applied until
the bound of `fT` gets tighter than `ϵ` or the number of steps reachs `max_iter`.
The returned bound corresponds to the improved polynomial bound with the remainder
of the `TaylorModelN` included.
"""
function linear_dominated_bounder(fT::TaylorModelN{N,T,S}; ϵ=1e-5, max_iter=5) where {N, T, S}
    d = one(T)
    dom = domain(fT)
    x0 = expansion_point(fT)
    pol = polynomial(fT)
    Pm = deepcopy(pol)
    bound = zero_interval(S)
    pol_hi = sup(pol(dom - x0))

    n_iter = 0
    x00 = Array{S}(undef, N)
    new_boxes = fill(bound, N)
    linear_coeffs = Array{Float64}(undef, N)
    while d > ϵ && n_iter < max_iter
        x00 .= mid(dom - x0)
        c = mid(dom)
        update!(Pm, x00)
        linear_part = Pm[1]
        if T <: Interval
            @. linear_coeffs = mid(linear_part.coeffs)
        else
            @. linear_coeffs = linear_part.coeffs
        end
        non_linear = nonlinear_polynomial(Pm)
        linear = Pm - non_linear
        centered_domain = dom .- c
        I1 = linear(centered_domain)
        Ih = non_linear(centered_domain)
        bound = inf(I1) + Ih
        d = diam(bound)
        n_iter += 1
        for (idx, box) in enumerate(dom)
            Li = linear_coeffs[idx]
            box_lo = inf(box)
            box_hi = sup(box)
            if Li == 0
                domi = box
            elseif Li < 0
                lo = max(box_hi - (d / abs(Li)), box_lo)
                domi = Interval(lo, box_hi)
            else
                hi = min(box_lo + (d / abs(Li)), box_hi)
                domi = Interval(box_lo, hi)
            end
            new_boxes[idx] = domi
        end
        x0 = dom
        dom = IntervalBox(new_boxes...)
    end

    return Interval(inf(bound), pol_hi) + remainder(fT)
end

"""
    quadratic_fast_bounder(fT::TaylorModel1)

Compute a *tighter* polynomial bound by the quadratic fast bounder.
The returned bound corresponds to the "improved" polynomial bound
with the remainder of the `TaylorModel1` included. This "improved" bound
can be one of the following two:
    1) An improved bound: if the domain of `fT` has a local minimizer,
       then an improved bound is returned.
    2) Original TaylorModel bound: if the local minimizer criteria is not
       satisfied, then the original bound of `fT` is returned.

This algorithm is a slight modification to the Makino & Berz algorithm.
For this algorithm the linear part is bounded by solving a simple
set of linear equations (compared to the iterative procedure by Makino & Berz).
"""
function quadratic_fast_bounder(fT::TaylorModel1)
    P = polynomial(fT)
    dom = domain(fT)
    x00 = expansion_point(fT)
    bound_tm = fT(dom - x00)
    signbit(P[2]) && return bound_tm

    cent_dom = dom - x00
    x0 = -P[1] / (2 * P[2])
    x = Taylor1(P.order)
    Qx0 = (x - x0) * P[2] * (x - x0)
    bound_qfb = (P - Qx0)(cent_dom)
    hi = sup(P(cent_dom))
    bound_qfb = Interval(inf(bound_qfb), hi) + remainder(fT)

    return bound_qfb ∩ bound_tm
end

"""
    quadratic_fast_bounder(fT::TaylorModelN)

Compute a *tighter* polynomial bound by the quadratic fast bounder.
The returned bound corresponds to the "improved" polynomial bound
with the remainder of the `TaylorModelN` included. This "improved" bound
can be one of the following two:
    1) An improved bound: if the domain of `fT` has a local minimizer,
       then an improved bound is returned.
    2) Original TaylorModel bound: if the local minimizer criteria is not
       satisfied, then the original bound of `fT` is returned.

This algorithm is a slight modification to the Makino & Berz algorithm.
For this algorithm the linear part is bounded by solving a simple
set of linear equations (compared to the iterative procedure by Makino & Berz).
"""
function quadratic_fast_bounder(fT::TaylorModelN)
    @assert get_numvars() == get_numvars(fT)
    P = polynomial(fT)
    dom = domain(fT)
    x0 = expansion_point(fT)
    H = Matrix(TaylorSeries.hessian(P))
    bound_tm = fT(dom - x0)
    if isposdef(H)
        P1 = -P[1].coeffs
        xn = H \ P1
        x = get_variables()#set_variables("x", numvars=length(xn))
        Qxn = 0.5 * (x - xn)' * H * (x - xn)
        bound_qfb = (P - Qxn)(dom - x0)
        hi = sup(P(dom - x0))
        bound_qfb = interval(inf(bound_qfb), hi) + remainder(fT)
        bound = bound_qfb ∩ bound_tm
        return bound
    end

    return bound_tm
end


text:
module ELM

using Reexport
@reexport using DataFrames

export ExtremeLearningMachine
export fit!, predict

# Base code
include("base.jl")

end


text:
include("header.jl")
@testset "highorder" begin
    # info("Test higher order gradients...")
    g1 = grad(sin); @test g1(1)==cos(1)
    g2 = grad(g1);  @test g2(1)==-sin(1)
    g3 = grad(g2);  @test g3(1)==-cos(1)
    g4 = grad(g3);  @test g4(1)==sin(1)
    g5 = grad(g4);  @test g5(1)==cos(1)
    g6 = grad(g5);  @test g6(1)==-sin(1)
    g7 = grad(g6);  @test g7(1)==-cos(1)
    g8 = grad(g7);  @test g8(1)==sin(1)
    g9 = grad(g8);  @test g9(1)==cos(1)

    # PR #75: Tape confusion fix
    @test grad(x -> x*grad(y -> x+y)(x))(5.0) == 1
    @test grad(x -> x*grad(y -> x+y)(1x))(5.0) == 1

    # WIP
    # @test (s->grad((@diff (x->x*(t->grad((@diff (y->x+y)(t)),t))(x))(s)),s))(Param(5)) == 1
    # @test (s->grad((@diff (x->x*(t->grad((@diff (y->x+y)(t)),t))(1x))(s)),s))(Param(5)) == 1
    # @test (s->grad(differentiate(x->x*(t->grad(differentiate(y->x+y,t),t))(x),s),s))(Param(5)) == 1
    # @test (s->grad(differentiate(x->x*(t->grad(differentiate(y->x+y,t),t))(1x),s),s))(Param(5)) == 1

    # Issue #44: third gradient of exp(x*x) gives nothing
    @test exp(1) == grad(exp)(1) == grad(grad(exp))(1)
    f(x) = exp(x*x)
    @test f(1) ≈ grad(f)(1) / 2 ≈ grad(grad(f))(1) / 6

    # Issue #62: bug second derivative tanh
    @test grad(tanh)(1) == 1 - tanh(1)^2
    @test grad(grad(tanh))(1) == -2*tanh(1)*(1-tanh(1)^2)

    # Issue Knet#439: hessians for neural networks
    hess(f,i=1) = grad((x...)->grad(f)(x...)[i])
    _nll(p,y)=(p=exp.(p);p=p./sum(p,dims=1);p=p[1,:];mean(-log.(p)))
    w,b,x,y = randn(2,3),randn(2),randn(3,4),randn(2,4)
    @test size(w) == size(hess(w->_nll(w*x.+b,y))(w))
    @test size(b) == size(hess(b->_nll(w*x.+b,y))(b))
    @test size(x) == size(hess(x->_nll(w*x.+b,y))(x))
end

nothing


text:
# Testing note: only run this file once per Julia session
# otherwise the in-memory or mmap'd files can get multiple
# sets of keys added to them and become too full
using Bloom

# First set up 9 sample Bloom filters using different constructors
# (Ordered as in bloom-filter.jl)

# Raw construction
bf0 = BloomFilter(BitVector(10), 4, 10, 0.01, 20, "")

# First group of constructors: capacity, bits per element, k
bf1 = BloomFilter(1000, 20, 4)
bf2 = BloomFilter(open("/tmp/test_array1.array", "w+"), 1000, 20, 4)
bf3 = BloomFilter("/tmp/test_array2.array", 1000, 20, 4)

# Second group of constructors: capacity, error rate, k
bf4 = BloomFilter(1000, 0.01, 5)
bf5 = BloomFilter(open("/tmp/test_array3.array", "w+"), 1000, 0.01, 5)
bf6 = BloomFilter("/tmp/test_array4.array", 1000, 0.01, 5)

# Third group of constructors: capacity and error rate only,
# computes optimal k from a space efficiency perspective
bf7 = BloomFilter(1000, 0.01)
bf8 = BloomFilter(open("/tmp/test_array5.array", "w+"), 1000, 0.01)
bf9 = BloomFilter("/tmp/test_array6.array", 1000, 0.01)

# Now create a larger in-memory Bloom filter and an mmap-backed one for testing
n = 100000
bfa = BloomFilter(n, 0.001, 5)
bfb = BloomFilter("/tmp/test_array_lg.array", n, 0.001, 5)

# Test with random strings
random_chars = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789"
test_keys = Array(String, n)
for i in 1:n
    temp_str = ""
    for j in 1:8
        temp_str = string(temp_str, random_chars[rand(1:62)])
    end
    test_keys[i] = temp_str
end

println("For insertions:")
@time(
for test_key in test_keys
    add!(bfa, test_key)
end
)

println("For lookups:")
@time(
for test_key in test_keys
    assert(contains(bfa, test_key))
end
)

println("For insertions (mmap-backed):")
@time(
for test_key in test_keys
    add!(bfb, test_key)
end
)

println("For lookups (mmap-backed):")
@time(
for test_key in test_keys
    assert(contains(bfb, test_key))
end
)

# Probabilistic tests – note these may fail every once in a while...
# (but very very rarely)
test_keys_p = Array(String, n)
for i in 1:n
    temp_str = ""
    for j in 1:9  # Longer so not in by definition
        temp_str = string(temp_str, random_chars[rand(1:62)])
    end
    test_keys_p[i] = temp_str
end

false_positives_a = 0
false_positives_b = 0
for test_key_p in test_keys_p
    if contains(bfa, test_key_p)
        false_positives_a += 1
    end
    if contains(bfb, test_key_p)
        false_positives_b += 1
    end
end

# Less than 2x requested false positives... (to avoid failing too often out of pure chance)
@printf "System is %s-bit\n" string(typeof(1))[4:5]
@printf "In-memory Bloom is %.2f%% full\n" (100 * sum(bfa.array) / bfa.n_bits)
@printf "%d false positives in %d tests\n" false_positives_a n
@printf "Error rate for in-memory Bloom filter %.2f%%\n" (bfa.error_rate * 100)
# assert((false_positives_a / n) <= (2.0 * bfa.error_rate))

@printf "Mmap'd Bloom is %.2f%% full\n" (100 * sum(bfb.array) / bfb.n_bits)
@printf "%d false positives in %d tests\n" false_positives_b n
@printf "Error rate for in-memory Bloom filter %.2f%%\n" (bfb.error_rate * 100)
# assert((false_positives_b / n) <= (2.0 * bfb.error_rate))
assert(false_positives_a == false_positives_b)  # Must be true since bfa and bfb should be identical at bit-level

# Test re-opening bfb
bfb = 0
gc()

bfb = BloomFilter(open("/tmp/test_array_lg.array", "r+"), n, 0.001, 5)
println("For lookups after re-opening (mmap-backed):")
@time(
for test_key in test_keys
    assert(contains(bfb, test_key))
end
)

bfb = 0
gc()

bfb = BloomFilter("/tmp/test_array_lg.array", n, 0.001, 5)
println("For lookups after re-opening second time (mmap-backed):")
@time(
for test_key in test_keys
    assert(contains(bfb, test_key))
end
)

## Note: This doesn't work as hash(x::String, seed::Int)
## is only defined for strings in dict.jl
# # Test insertions of non-string types
# test_other_a = 17    # Int
# test_other_b = 15.6  # Float
# test_other_c = "String" #("Tuples", "of", "strings")

# # Test alias'd command
# add!(bfb, test_other_a)
# add!(bfb, test_other_b)
# add!(bfb, test_other_c)

# assert(contains(bfb, test_other_a))
# assert(contains(bfb, test_other_b))
# assert(contains(bfb, test_other_b))

# Clean up mmap-backed temp files (otherwise can end up re-opening them and writing multiple key sets to one file!)
rm("/tmp/test_array1.array")
rm("/tmp/test_array2.array")
rm("/tmp/test_array3.array")
rm("/tmp/test_array4.array")
rm("/tmp/test_array5.array")
rm("/tmp/test_array6.array")
rm("/tmp/test_array_lg.array")


text:
#-------------------------------------------------------------------------------
# File: genetic.jl
# Description: This files contains all functions that are used for genetic.
# Date: November 17, 2019
# Authors: Killian Fretaud.
#-------------------------------------------------------------------------------

"""

"""
function selection(roulette_wheel, population, instance::Instance, population_size::Int)
    parent_1 = rand(1:roulette_wheel[population_size])
    i = 1
    while parent_1 > roulette_wheel[i] && i < population_size
        i += 1
    end
    parent_1 = i
    parent_2 = parent_1
    while parent_2 == parent_1
        parent_2 = rand(1:roulette_wheel[population_size])
        i = 1
        while parent_2 > roulette_wheel[i] && i < population_size
            i += 1
        end
        parent_2 = i
    end

    return parent_1, parent_2
end

"""

"""
function crossover(parent_1, parent_2, population, instance::Instance)
    index = rand(1:instance.nb_jobs)
    start = parent_1[index]

    offspring_1 = zeros(Int, instance.nb_jobs) ; offspring_1[index] = start
    offspring_2 = zeros(Int, instance.nb_jobs) ; offspring_2[index] = parent_2[index]

    while parent_2[index] != start
        i = 1
        while (parent_1[i] != parent_2[index])
            i += 1
        end
        index = i
        offspring_1[index] = parent_1[index]
        offspring_2[index] = parent_2[index]
    end

    for index in 1:instance.nb_jobs
        if offspring_1[index] == 0
            offspring_1[index] = parent_2[index]
            offspring_2[index] = parent_1[index]
        end
    end

    return offspring_1, offspring_2
end

"""

"""
function mutation(genotype, phenotype)
    n = length(genotype)
    indice_1 = rand(1:n) ; indice_2 = indice_1
    while indice_1 == indice_2
        indice_2 = rand(1:n)
    end

    gene_1 = genotype[indice_1] ; gene_2 = genotype[indice_2]
    genotype[indice_1] = gene_2 ; genotype[indice_2] = gene_1

    for k in 1:length(phenotype)
        if phenotype[k] == gene_1
            phenotype[k] = gene_2
        elseif phenotype[k] == gene_2
            phenotype[k] = gene_1
        end
    end

    return genotype, phenotype
end

"""

"""
function evaluation(phenotype, instance::Instance)
    m = instance.nb_machines
    K = length(phenotype)
    n = instance.nb_jobs
    occurence = zeros(Int, n) # sur quelle machine est la tâche
    machines  = zeros(Int, m) # quelle tâche est sur la machine
    C_of_m = zeros(Int, m)
    C_of_task = zeros(Int, m)
    time_robot = 0 ; travel = 0 # last machine visited

    # For all move
    for pos in 1:K
        # A job goes to a machine
        job = phenotype[pos]
        occurence[job] = occurence[job] + 1
        j = occurence[job]

        # Compute the robot's set-up time
        travel_time = 0
        for i in j:travel
            travel_time = travel_time + instance.travel_times[i]
        end
        for i in (travel+1):(j-1)
            travel_time = travel_time + instance.travel_times[i]
        end
        travel = j

        # Is it a machine or the end of our cycle ?
        if j < m+1
            # Process start when :
            #   - machine is available
            #   - job is available and delivered,
            #   - robot has finish is set-up and delivered the job
            debut = max(C_of_m[j], max(C_of_task[job], time_robot + travel_time) + instance.travel_times[j])
            time_robot = debut
            C_of_m[j] = debut + instance.processing_times[job, j]
            C_of_task[job] = C_of_m[j]
        else
            debut = max(C_of_task[job], time_robot + travel_time) + instance.travel_times[j]
            time_robot = debut
        end
    end

    return time_robot
end

"""

"""
# should give cmax without more calcul after its call
function construction_and_evaluation(genotype, instance::Instance)
    #Step 1: Initialization
    n = instance.nb_jobs ; m = instance.nb_machines
    phenotype = [] ; pos = 1 # position a remplir
    occurence = zeros(Int, n) # sur quelle machine est la tâche
    machines  = zeros(Int, m) # quelle tâche est sur la machine

    # evaluation
    C_of_m = zeros(Int, m)
    time_robot = 0 ; travel = 0 # last machine visited
    # end evaluation

    # At first, there is only one job available
    E = [1] # set of indexes

    #Step 2: General procedure
    while pos <= n*(m+1)
        # Next job is chosen randomly
        k = rand(1:length(E))
        indice = E[k] ; job = genotype[indice]
        # We push it
        push!(phenotype, job) ; pos = pos + 1
        occurence[job] = occurence[job] + 1
        j = occurence[job]

        # Was the job on a machine ?
        if j > 1 && j <= m
            machines[j-1] = 0
        end

        # Is next job available ?
        if occurence[job] == 1 && indice < n
            push!(E, indice+1)
        end

            # evaluation
                travel_time = 0
                for i in j:travel
                    travel_time = travel_time + instance.travel_times[i]
                end
                for i in (travel+1):(j-1)
                    travel_time = travel_time + instance.travel_times[i]
                end
                travel = j

                if j <= m
                    temp = 0
                    if j > 1
                        temp = C_of_m[j-1]
                    end
                    debut = max(C_of_m[j], max(temp, time_robot + travel_time) + instance.travel_times[j])
                    time_robot = debut
                    C_of_m[j] = debut + instance.processing_times[job, j]
                else
                    debut = max(C_of_m[j-1], time_robot + travel_time) + instance.travel_times[j]
                    time_robot = debut
                end
            # end evaluation

        # The machine was not available
        while (j <= m && machines[j] != 0)
            # Swap
            Tk = machines[j] ; machines[j] = job
            job = Tk
            # Again we push it
            push!(phenotype, job) ; pos = pos + 1
            occurence[job] = occurence[job] + 1
            j = occurence[job]

                # evaluation
                    travel = j
                    if j <= m
                        debut = max(C_of_m[j], time_robot + instance.travel_times[j])
                        time_robot = debut
                        C_of_m[j] = debut + instance.processing_times[job, j]
                    else
                        time_robot = time_robot + instance.travel_times[j]
                    end
                # end evaluation
        end

        if j == m + 1
            deleteat!(E, 1)
            if machines[m] == job
                machines[m] = 0
            end
        elseif j <= m
            machines[j] = job
        end
    end

    return phenotype, time_robot
end

"""

"""
function evaluate_genotype(genotype, instance::Instance, nb_of_constructions::Int)
    phenotype, score = construction_and_evaluation(genotype, instance)
    for j in 1:nb_of_constructions
        p, s = construction_and_evaluation(genotype, instance)
        if s < score
            phenotype = p ; score = s
        end
    end
    return (genotype, phenotype, score)
end

"""

"""
function first_generation(instance::Instance, population_size::Int, nb_of_constructions::Int)
    population = []
    for i in 1:population_size
        genotype = randperm(instance.nb_jobs)
        push!(population, evaluate_genotype(genotype, instance, nb_of_constructions))
    end
    return population
end

"""
    genetic(instance::Instance,
            probability_of_crossover::Float64 = 0.8,
            probability_of_mutation::Float64 = 0.2,
            population_size::Int = 100, nb_max_of_generation::Int = 100,
            nb_of_constructions::Int = 5)

Applies the genetic algorithm on `instance`. Uses a few parameters. Default
values are from the article. Return the best solution found.
"""
function genetic(instance::Instance,
                 probability_of_crossover::Float64 = 0.8,
                 probability_of_mutation::Float64 = 0.2,
                 population_size::Int = 100, nb_max_of_generation::Int = 100,
                 nb_of_constructions::Int = 5
                 ;
                 verbose = false)
    if verbose
        println("----------------------------  START of GA's execution")
    end

    # inital population
    population = first_generation(instance, population_size, nb_of_constructions)
    sort!(population, by = x -> x[3])

    maximum_score = population[population_size][3]
    score_for_roulette_wheel = zeros(Int, population_size)
    total = 0
    for i in 1:population_size
        score_for_roulette_wheel[i] = maximum_score .- population[i][3]
        total += score_for_roulette_wheel[i]
        score_for_roulette_wheel[i] = total
    end

    # time pass
    generation = 1
    while generation <= nb_max_of_generation
        if verbose
            if (generation-1) % 10 == 0
                println("it. ", generation)
            end
        end

        nb_offsprings = 0
        while nb_offsprings < (population_size * probability_of_crossover)

            p1, p2 = selection(score_for_roulette_wheel, population, instance, population_size)

            o1, o2 = crossover(population[p1][1], population[p2][1], population, instance)
            o1, o1_phenotype, o1_score = evaluate_genotype(o1, instance, nb_of_constructions)
            o2, o2_phenotype, o2_score = evaluate_genotype(o2, instance, nb_of_constructions)

            if rand() <= probability_of_mutation
                o1, o1_phenotype = mutation(o1, o1_phenotype)
                o1_score = evaluation(o1_phenotype, instance)
            end
            if rand() <= probability_of_mutation
                o2, o2_phenotype = mutation(o2, o2_phenotype)
                o2_score = evaluation(o2_phenotype, instance)
            end

            push!(population, (o1, o1_phenotype, o1_score))
            push!(population, (o2, o2_phenotype, o2_score))

            nb_offsprings = nb_offsprings + 2
        end

        sort!(population, by = x -> x[3])
        for i in 1:nb_offsprings
            pop!(population)
        end

        if verbose
            if generation % 10 == 0
                println("\tbest  Cmax: ", population[1][3])
                println("\tworst Cmax: ", population[100][3])
            end
        end

        generation = generation + 1
    end
    if verbose
        println()
        println()
        println("----------------------------  END of GA's execution")
        println("Best solution: ", population[1][1])
        println("\t", population[1][2])
        println("Cmax is ", population[1][3])
    end

    return population[1]
end

"""

"""
function construction(genotype, instance)
    #Step 1: Initialization
    n = instance.nb_jobs ; m = instance.nb_machines
    phenotype = [] ; pos = 1 # position a remplir
    occurence = zeros(Int, n) # sur quelle machine est la tâche
    machines  = zeros(Int, m) # quelle tâche est sur la machine

    # At first, there is only one job available
    E = [1] # set of indexes

    #Step 2: General procedure
    while pos <= n*(m+1)

        # Next job is chosen randomly
        k = rand(1:length(E))
        indice = E[k] ; job = genotype[indice]
        # We push it
        push!(phenotype, job) ; pos = pos + 1
        occurence[job] = occurence[job] + 1
        j = occurence[job]

        # Was the job on a machine ?
        if j > 1 && j <= m
            machines[j-1] = 0
        end

        # Is next job available ?
        if occurence[job] == 1 && indice < n
            push!(E,indice+1)
        end

        # The machine was not available
        while (j <= m && machines[j] != 0)
            # Swap
            Tk = machines[j] ; machines[j] = job
            job = Tk
            # Again we push it
            push!(phenotype, job) ; pos = pos + 1
            occurence[job] = occurence[job] + 1
            j = occurence[job]
        end

        if j == m+1
            deleteat!(E,1)
            if machines[m] == job
                machines[m] = 0
            end
        elseif j <= m
            machines[j] = job
        end
    end

    return phenotype
end

max_stars_repo_name: k-fretaud/Flowshop-robotic-cell

max_issues_repo_name: k-fretaud/Flowshop-robotic-cell

max_forks_repo_name: k-fretaud/Flowshop-robotic-cell


text:
include(joinpath(pwd(), "src/models/biped.jl"))
include(joinpath(pwd(), "src/objectives/velocity.jl"))
include(joinpath(pwd(), "src/objectives/acceleration.jl"))
include(joinpath(pwd(), "src/constraints/contact.jl"))

# Visualize
# - Pkg.add any external deps from visualize.jl
include(joinpath(pwd(), "src/models/visualize.jl"))
vis = Visualizer()
open(vis)

# Configurations
# 1: x pos
# 2: z pos
# 3: torso angle (rel. to downward vertical)
# 4: thigh 1 angle (rel. to downward vertical)
# 5: calf 1 (rel. to thigh 1)
# 6: thigh 2 (rel. to downward vertical)
# 7: calf 2 (rel. to thigh 2)
θ = pi / 10.0
# q1 = initial_configuration(model, θ) # generate initial config from θ
# qT = copy(q1)
# qT[1] += 1.0
# q1[3] -= pi / 30.0
# q1[4] += pi / 20.0
# q1[5] -= pi / 10.0
q1, qT = loop_configurations(model, θ)
qT[1] += 1.0
visualize!(vis, model, [qT])

# Horizon
T = 11

# Time step
tf = 1.0
h = tf / (T - 1)

# Bounds

# control
# u = (τ1..4, λ1..2, β1..4, ψ1..2, η1...4, s1)
# ul <= u <= uu
_uu = Inf * ones(model.m)
_uu[model.idx_u] .= 100.0

_ul = zeros(model.m)
_ul[model.idx_u] .= -100.0
ul, uu = control_bounds(model, T, _ul, _uu)

xl, xu = state_bounds(model, T, x1 = [q1; q1])

# Objective
q_ref = linear_interp(q1, qT, T)
X0 = configuration_to_state(q_ref)

# penalty on slack variable
obj_penalty = PenaltyObjective(1.0e3, model.m)

Qq = 1.0 * Diagonal(ones(model.nq)) # configuration penalty
Q = cat(0.5 * Qq, 0.5 * Qq, dims = (1, 2)) # state penalty
QT = cat(0.5 * Qq, 100.0 * Diagonal(ones(model.nq)), dims = (1, 2))
R = Diagonal([1.0e-1 * ones(model.nu)..., zeros(model.m - model.nu)...]) # no penalty on "contact" controls

# quadratic tracking objective
# Σ (x - xref)' Q (x - x_ref) + (u - u_ref)' R (u - u_ref)
obj_tracking = quadratic_tracking_objective(
    [t < T ? Q : QT for t = 1:T],
    [R for t = 1:T-1],
    [X0[t] for t = 1:T],
    [zeros(model.m) for t = 1:T]
    )
# quadratic velocity penalty
# Σ v' Q v
obj_velocity = velocity_objective(
    [Diagonal(1.0 * ones(model.nq)) for t = 1:T-1],
    model.nq,
    h = h)
# quadratic acceleration penalty
# Σ a' Q a
# obj_acceleration = acceleration_objective(
#     [Diagonal(1.0e-1 * ones(model.nq)) for t = 1:T-1],
#     model.nq,
#     h = h)
obj = MultiObjective([obj_tracking,
                      obj_penalty,
                      obj_velocity])#,
                      # obj_acceleration])

# Constraints
con_contact = contact_constraints(model, T)
# con_pinned = pinned_foot_constraint(model, q1, T)
con = multiple_constraints([con_contact])#, con_pinned])

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
               h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con
               )

# trajectory initialization
U0 = [1.0e-5 * rand(model.m) for t = 1:T-1] # random controls

# Pack trajectories into vector
Z0 = pack(X0, U0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem
@time Z̄ = solve(prob, copy(Z0),
    nlp = :ipopt,
    tol = 1.0e-3, c_tol = 1.0e-3)#, mapl = 5)

check_slack(Z̄, prob)
X̄, Ū = unpack(Z̄, prob)

# using Plots
# plot(hcat(Ū...)[1:model.nu,:]', linetype = :steppost)

# Visualize
visualize!(vis, model, state_to_configuration(X̄), Δt = h)

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
using OED
using jInv.Utils
using jInv.ForwardShare
using Base.Test
using jInv.Mesh

n = 100
nrhs = 5

A = sdiag(1./linspace(1,n,n).^4)
A = A[1:2:end,:]
N = size(A,1)
x = rand(n,nrhs); x[4:7,:] = 1;
D = A*x  + 1e-8* randn(N,nrhs)/1e3
OEDparamA = OEDTypeAParam(A,D,zeros(0,0))

# setup OED type B param
r = 20;
l = 2;
Mimg  = getRegularMesh([-1 1 -1 1.],[10,10])
Mgrid = getRegularMesh([-1.2 2.4 -1.6 1.4],[20,20])
OEDparamB = getOEDTypeBParam(r,l,x,Mimg,Mgrid,noiseLevel=0.)
OEDparamC = getOEDTypeCParam(r,l,x,Mimg,Mgrid,noiseLevel=0.)

L = 1e-1*speye(n)
pForA = UnconstrainedOEDParam(OEDparamA,L,[],[])
pForB = UnconstrainedOEDParam(OEDparamB,L,[],[])
pForC = UnconstrainedOEDParam(OEDparamC,0*L,[],[])

pFors = (pForA,pForB,pForC)

println("\n\t---test OED with unconstrained lower-level problem---")
for k=1:3
	display(pFors[k])
	np  =  getNumberOfDesignParameters(pFors[k].Design)
	p0      = rand(np)

	println("\t\tcheck reconstruction")
	xr, = getData(p0,pFors[k])
	println("\t\t\trel.error: $(norm(xr-x)/norm(x))")

	print("\t\tcheckDerivative...")

	pass, = checkDerivative(p0,pFors[k])
	@test pass
	print("passed!\n")

	print("\t\tadjoint test...")
	(nw,nv) = getSensMatSize(pFors[k])
	w = rand(nw)
	v = rand(nv)
	t1 = dot(w,getSensMatVec(v,p0,pFors[k]))
	t2 = dot(v,getSensTMatVec(w,p0,pFors[k]))
	@test norm(t1-t2)/norm(t1)<1e-10
	print("passed!\n")
end
println("\t---least squares IP: sensitivities OK ---")


text:
module PyJuliaTester

import Python_jll
using PyVenvs: @pyvenv, PyVenvs

@pyvenv TEST_VENV requirements="julia[test]"

function test_cmd(args = ``; rebuild::Union{Bool,Nothing} = nothing)
    py = PyVenvs.python(TEST_VENV)
    jl = Base.julia_cmd().exec[1]
    cmd = `$py -m julia.runtests $args`
    env = copy(ENV)
    # Using `PYJULIA_TEST_RUNTIME` instead of `--julia-runtime=$jl` as
    # `with_rebuilt` ignores `--julia-runtime` ATM.
    env["PYJULIA_TEST_RUNTIME"] = jl
    env["JULIA_PROJECT"] = Base.active_project()
    if rebuild !== nothing
        # https://pyjulia.readthedocs.io/en/latest/testing.html
        env["PYJULIA_TEST_REBUILD"] = rebuild ? "yes" : "no"
    end
    cmd = setenv(cmd, env)
    return cmd
end

"""
    PyJuliaTester.runtests([command_args::Cmd]; rebuild::Bool::Bool)

Run `python -m julia.runtests \$command_args`.

If `rebuild == true` (`false`), set `PYJULIA_TEST_REBUILD=yes` (`no`).
See PyJulia's documentation for the effect.

See also:
* [`PyJuliaTester.help`](@ref)
* https://pyjulia.readthedocs.io/en/latest/testing.html
"""
function runtests(args = ``; kwargs...)
    Python_jll.python() do python
        PyVenvs.init(TEST_VENV; python = python)
        cmd = test_cmd(args; kwargs...)
        cmd_print = `$(cmd.exec)`
        @info "Run: $cmd_print"
        run(cmd)
    end
end

"""
    PyJuliaTester.help()

Show help message for PyJulia's `julia.runtests` API.
"""
function help()
    runtests(`--help`)
    return
end

end


text:
#
# 1. Custom initial condition setup
#    - ClimateMachine.HydrostaticBoussinesq.ocean_init_state!
#

# Setup types and globals we need
using StaticArrays
og=ClimateMachine.HydrostaticBoussinesq.OceanGyre

"""
    ocean_init_state!(::og)

initialize u,v,η with 0 and θ linearly distributed between 9 at z=0 and 1 at z=H

# Arguments
- `p`: OceanGyre problem object, used to dispatch on and obtain ocean height H
- `Q`: state vector
- `A`: auxiliary state vector, not used
- `coords`: the coordidinates
- `t`: time to evaluate at, not used
"""
function ClimateMachine.HydrostaticBoussinesq.ocean_init_state!(p::OceanGyre, Q, A, coords, t)
    @inbounds y = coords[2]
    @inbounds z = coords[3]
    @inbounds H = p.H

    Q.u = @SVector [0, 0]
    Q.η = 0

    # Wavy
    # Q.θ = (5 + 4 * cos(y * π / p.Lʸ)) * (1 + z / H)

    # Constant everywhere
    # Q.θ = 9

    # Big number
    # Q.θ = 50

    # Unstable column
    #
    # Define function with inversion near surface and stable stratification below
    #
    # $ \theta(z) = Ae^{\frac{-(z+L)}{L}} - Ce^{\frac{-D(z+L)}{L}} + B + Ez$
    #
    H=4000.;L=H/10.;A=20.;C=50.;D=2.5;B=8.;E=5.e-4;
    ft(xx,L)=exp(-xx/L);
    th1(zz)=A*ft.(-zz .+ L, L);
    th2(zz)=C*ft(D*(-zz .+ L), L);
    phi1=th1.(z)
    phi2=th2.(z)
    Q.θ = phi1 .- phi2 .+ B .+ E .* z;

    return nothing
end

max_stars_repo_name: christophernhill/climahacks

max_issues_repo_name: christophernhill/climahacks

max_forks_repo_name: christophernhill/climahacks


text:
@testset "Constructors" begin
    #
    # missings
    #

    @test ismissing(Missing())
    @test ismissing(missing)

    #
    # DataVector's
    #

    dv = DataArray([1, 2, 3], falses(3))
    @test isequal(dv.data, [1, 2, 3])
    @test isequal(dv.na, falses(3))
    @test isequal(dv, DataArray([1, 2, 3], [false, false, false]))
    @test isequal(dv, DataArray([1, 2, 3]))

    dv = convert(DataArray, trues(3))
    @test isequal(dv.data, [true, true, true])
    @test isequal(dv.na, falses(3))
    @test isequal(dv, convert(DataArray, trues(3)))

    dv = DataArray([1, 2, 3], falses(3))
    @test isequal(dv, convert(DataArray, 1:3))

    dv = DataArray(Int, 3)
    @test isequal(eltype(dv), Union{Int,Missing})
    @test isequal(dv.na, trues(3))

    dv = convert(DataArray, zeros(3))
    @test isequal(dv, convert(DataArray, zeros(3)))

    dv = convert(DataArray, zeros(Int, 3))
    @test isequal(dv, convert(DataArray, zeros(Int, 3)))

    dv = convert(DataArray, ones(3))
    @test isequal(dv, convert(DataArray, ones(3)))

    dv = convert(DataArray, ones(Int, 3))
    @test isequal(dv, convert(DataArray, ones(Int, 3)))

    dv = convert(DataArray, falses(3))
    @test isequal(dv, convert(DataArray, falses(3)))

    dv = convert(DataArray, trues(3))
    @test isequal(dv, convert(DataArray, trues(3)))

    @test DataArray([1, missing]) isa DataVector{Int}
    @test isequal(DataArray([1, missing]), [1, missing])
    @test DataArray{Int}([1, missing]) isa DataVector{Int}
    @test isequal(DataArray{Int}([1, missing]), [1, missing])
    @test DataArray{Any}([1, missing]) isa DataVector{Any}
    @test isequal(DataArray{Any}([1, missing]), [1, missing])
    @test DataArray{Int, 1}([1, missing]) isa DataVector{Int}
    @test isequal(DataArray{Int, 1}([1, missing]), [1, missing])
    @test DataArray{Any, 1}([1, missing]) isa DataVector{Any}
    @test isequal(DataArray{Any, 1}([1, missing]), [1, missing])

    #
    # PooledDataArray's
    #

    pdv = PooledDataArray([1, 2, 3], falses(3))
    @test all(pdv .== [1, 2, 3])
    @test all(ismissing.(pdv) .== falses(3))

    @test isequal(pdv, PooledDataArray([1, 2, 3], [false, false, false]))
    @test isequal(pdv, PooledDataArray([1, 2, 3]))

    pdv = convert(PooledDataArray, trues(3))
    @test all(pdv .== [true, true, true])
    @test all(ismissing.(pdv) .== falses(3))
    @test isequal(pdv, convert(PooledDataArray, trues(3)))

    pdv = PooledDataArray([1, 2, 3], falses(3))
    @test isequal(pdv, convert(PooledDataArray, 1:3))
    @test isequal(pdv, convert(PooledDataArray, PooledDataArray([1, 2, 3])))

    pdv = PooledDataArray(Int, 3)
    @test isequal(eltype(pdv), Union{Int, Missing})
    @test all(ismissing.(pdv) .== trues(3))

    pdv = convert(PooledDataArray, zeros(3))
    @test isequal(pdv, convert(PooledDataArray, zeros(3)))

    pdv = convert(PooledDataArray, zeros(Int, 3))
    @test isequal(pdv, convert(PooledDataArray, zeros(Int, 3)))

    pdv = convert(PooledDataArray, ones(3))
    @test isequal(pdv, convert(PooledDataArray, ones(3)))

    pdv = convert(PooledDataArray, ones(Int, 3))
    @test isequal(pdv, convert(PooledDataArray, ones(Int, 3)))

    pdv = convert(PooledDataArray, falses(3))
    @test isequal(pdv, convert(PooledDataArray, falses(3)))

    pdv = convert(PooledDataArray, trues(3))
    @test isequal(pdv, convert(PooledDataArray, trues(3)))

    #
    # DataMatrix
    #

    dm = DataArray([1 2; 3 4], falses(2, 2))
    @test isequal(dm.data, [1 2; 3 4])
    @test isequal(dm.na, falses(2, 2))

    @test isequal(dm, DataArray([1 2; 3 4], [false false; false false]))
    @test isequal(dm, DataArray([1 2; 3 4]))

    dm = convert(DataArray, trues(2, 2))
    @test isequal(dm.data, trues(2, 2))
    @test isequal(dm.na, falses(2, 2))

    @test isequal(dm, convert(DataArray, trues(2, 2)))

    dm = DataArray(Int, 2, 2)
    @test isequal(eltype(dm), Union{Int, Missing})
    @test isequal(dm.na, trues(2, 2))

    @test_nowarn convert(DataArray, zeros(2, 2))

    @test_nowarn convert(DataArray, zeros(Int, 2, 2))

    @test_nowarn convert(DataArray, ones(2, 2))
    @test_nowarn convert(DataArray, ones(Int, 2, 2))

    @test_nowarn convert(DataArray, falses(2, 2))
    @test_nowarn convert(DataArray, trues(2, 2))

    @test_nowarn convert(DataArray, eye(3, 2))
    @test_nowarn convert(DataArray, eye(2))
    @test_nowarn convert(DataArray, diagm(Float64[pi, pi]))

    @test DataArray([1 missing]) isa DataMatrix{Int}
    @test isequal(DataArray([1 missing]), [1 missing])
    @test DataArray{Int}([1 missing]) isa DataMatrix{Int}
    @test isequal(DataArray{Int}([1 missing]), [1 missing])
    @test DataArray{Any}([1 missing]) isa DataMatrix{Any}
    @test isequal(DataArray{Any}([1 missing]), [1 missing])
    @test DataArray{Int, 2}([1 missing]) isa DataMatrix{Int}
    @test isequal(DataArray{Int, 2}([1 missing]), [1 missing])
    @test DataArray{Any, 2}([1 missing]) isa DataMatrix{Any}
    @test isequal(DataArray{Any, 2}([1 missing]), [1 missing])
end


text:
const scores_list = [(1, "aeioulnrst"), (2, "dg"), (3, "bcmp"),
    (4, "fhvwy"), (5, "k"), (8, "jx"), (10, "qz")]

const scores =  Dict(c=>score for (score, chars) in scores_list for c in chars)

function score(str)
    sum(map(c -> get(scores, c, 0), collect(lowercase(str))))
end

max_stars_repo_name: koladilip/julia-exercism

max_issues_repo_name: koladilip/julia-exercism

max_forks_repo_name: koladilip/julia-exercism


text:
# Surface flux calculations

"""
    calculate_land_sfc_fluxes(model::ModelSetup, state, aux, t; MO_params)

Calculate surface fluxes using the bulk gradient diffusion theory.
"""
function calculate_land_sfc_fluxes(model::ModelSetup, state, aux, t; MO_params = nothing) # should pass in the coupler state (also move to coupler), so can access states of both models derectly -e.g. callback?

    FT = eltype(state)
    p = model.parameters

    p_sfc = p.pₒ
    θ_a = state.ρθ / state.ρ
    T_a = θ_a
    T_sfc = aux.T_sfc
    θ_sfc = T_sfc

    R_SW = (FT(1) - p.α) * p.τ * p.F_sol * (FT(1) .+ sin(t * 2π / p.τ_d))
    R_LW = p.ϵ * (p.σ * θ_sfc .^ 4 - p.F_a)
    SH = state.ρ * p.cp_d * p.g_a * (T_sfc - T_a) # g_a could be substituded by g_a=f(MO_params) (see Bonan P90), but first need to modify SurfaceFluxes.jl

    #LH   = p.lambda * p.g_w * (q_sat(T_sfc, p_sfc) - q_a)
    F_tot = -(R_SW - R_LW - SH)#- LH

end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
#   Unit examples of the FourierAnalysis Package for julia language
#
#   MIT License
#   Copyright (c) 2019-2022,
#   Marco Congedo, CNRS, Grenobe, France:
#   https://sites.google.com/site/marcocongedo/home

# ? CONTENTS :
#   This example shows how to compute and plot spectra
#   and how to extract inforation from them.

#   ~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~  #
#                                                                             #
#   ~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~  #

using FourierAnalysis, FFTW, LinearAlgebra, Statistics,
      Plots, Plots.PlotMeasures

# add module for reading the two EEG text files to be used ater
push!(LOAD_PATH, @__DIR__)
using IOtxt

# Get EEG file names with complete path (they have extension .txt)
S=getFilesInDir(@__DIR__; ext=(".txt",))

# read the two EEG data files and put them in a Matrix object
X1=readEEG(S[1])
X2=readEEG(S[2])

######## compute the FFT of one epoch of a time series ########
sr, t, f, a = 128, 128, 8, 2
v=sinusoidal(a, f, sr, t, 0.6)
plot(v) # plot the time series
tapering=hann

### The long way using FFTW.jl only:
# get a plan.  The normalization 2/t outputs correct peak amplitude
P=plan_rfft(v)*(2/t)

# get a tapering window
a=taper(tapering, t)

# plot the tapered signal
plot!(v.*a.y)

# do the FFT on tapered data
Ψ=P*(v.*a.y)

# get the amplitude spectrum
Σ=abs.(Ψ)
bar!(Σ[2:end])

# check the amplitude at frequency f and neighbooring bins
pos=f2b(f, sr, t, DC=true)

Σ[pos-1]
Σ[pos]
Σ[pos+1]

# The fast way using FourierAnalysis.jl:
Σ2=spectra(v, sr, t; tapering=tapering, func=√)
bar!(Σ2.y)

# equivalently:
Σ2=spectra(v, sr, t; tapering=taper(hann, t), func=√)
bar!(Σ2.y)

####################################################

# test spectra with odd samples and with DC
sr=65; t=65
f, a = fres(sr, t)*4, 2
v=sinusoidal(a, f, sr, t, 0.6; DC=3)
plot(v) # plot the time series
Σ2=spectra(v, sr, t; tapering=rectangular, DC=true, func=√)
bar!(Σ2.y)

####################################################
### Check amplitude spectra at all Fourier Frequency (rectangular taper) ###
# the amplitudes are in increment of 10 along frequencies
# NB when t is even, correct amplitude for the last frequency is obtained
# only if the sinusoidal has a phase=π/6
sr, t, = 16, 32
t½ = t÷2
V=Matrix{Float64}(undef, t, t½)
for i=1:t½ V[:, i]=sinusoidal(10*i, b2f(i, sr, t), sr, t, π/6) end

# plot the sinusoids in V
Z=deepcopy(V)
for i=2:size(Z, 2), j=1:size(Z, 1) Z[j, i]+=(10.0*i) end
plot(Z, label="")

# long way using FFTW only.jl
P=plan_rfft(V, 1)*(2/t)
Ψ=P*V
Ψ[end, :]*=0.5*2^0.25 # correction for Nyquist frequency
Σ=abs.(Ψ)
bar(Σ[brange(t, DC=true), :], labels="")

# fast way using FourierAnalysis.jl
Σ2=spectra(V, sr, t; tapering=rectangular, func=√, DC=true)
bar(Σ2.y[brange(t, DC=true), :], labels="")

#############################################################################
### Check amplitude spectra on long signals obtained by welch methods
# one sinusoidal is at an exact discrete Fourier Frequency and the other not
# Rectangular window
sr, t, f, a = 128, 128, 10, 0.5
v=sinusoidal(a, f, sr, t*16)+sinusoidal(a, f*3.5+0.5, sr, t*16)+randn(t*16)
Σ=spectra(v, sr, t; tapering=rectangular, func=√)
bar(Σ.y, labels="rectangular")

# harris4 window (default)
Σ2=spectra(v, sr, t; func=√)
bar!(Σ2.y, labels="harris4")

#smooth spectra
Σ3=smooth(blackmanSmoother, Σ2)
bar!(Σ3.y, labels="smoothed")

######## Check amplitude spectrum obtained with dpss (slepian tapers) ########
sr, t = 128, 128
t½ = t÷2
V=Matrix{Float64}(undef, t, t½)
for i=1:t½ V[:, i]=sinusoidal(1, b2f(i, sr, t), sr, t, 0 ) end
# plot half of the the first 7 sinusoids in V
plot(V[1:t½+1, 1:7], color=:gray, legend=false)

Σ=spectra(V, sr, t; tapering=taper(slepian, t; α=2, n=3), func=√)
plot(Σ.y, labels="")
#############################################################

#### Compute Welch amplitude spectra of EEG data ####
t, sr, slide, tapering = 1024, 128, 512, harris4

# gather some attributes to obtain nice spectra plots
spectraArgs=(fmax          = 48,
             left_margin   = 2mm,
             bottom_margin = 2mm,
             xtickfont     = font(10, "Times"),
             ytickfont     = font(10, "Times"))

S=spectra(X1, sr, t; tapering=tapering, slide=sr, func=sqrt)
plot(S; ytitle="Amplitude (\\muV)", spectraArgs...)

#smooth spectra
S2=smooth(blackmanSmoother, S)
plot(S2; ytitle="Amplitude (\\muV)", spectraArgs...)

# extract spectra in alpha range (8Hz to 12Hz) at all electrodes
e=extract(S, (8, 12))

# You can use any combination of integer and real numbers
e=extract(S, (8, 12.5))
e=extract(S, (8.5, 12))

# the following equivalently extract the spectra at 10Hz only (1 bin)
e=extract(S, 10)
e=extract(S, (10))
e=extract(S, (10, 10))

# extract spectra at all frequencies. Equivalent to S.y
e=extract(S, :)

# All these ways to indicate frequency ranges work with the mean function
# as well, e.g.,

# compute average spectra in alpha range for each time-series
bar(mean(S, (8, 12)))
# or
bar(mean(S, (8.0,12.0)))

# this actually extract the power at 10Hz for each time series
# and is equivalent to extract(S, 10)
mean(S, 10)
# thus this computes the average power at 10Hz across time-series
mean(extract(S, 10)) # or mean(mean(S, 10))
# and this computes the mean across time-series and across frequencies in range [8Hz, 12Hz]
mean(extract(S, (8, 12))) # or equivalently mean(mean(S, (8, 12)))

# average spectrum across all frequencies for each time-series
bar(mean(S, :))

# plot average spectra in 2Hz-band-pass regions for all time-series
plot(bands(S, 2))

# plot average spectra in in 2Hz-band-pass regions for series 15-19
plot(bands(S, 2)[:, 15:19])

# plot average spectra in in 2Hz-band-pass regions for series 19
plot(bands(S, 2)[:, 19])

# plot average spectra in alpha range for series 15-19
bar(mean(S, (8,12))[15:19])

# extract spectra in alpha range for series 15-19
e=extract(S, (8, 12))[:, 15:19]

# Compute several spectra altogether
𝐗=[X1, X2]
𝐒=spectra(𝐗, sr, t; tapering=tapering, slide=sr, func=√)
plot(𝐒[1]; ytitle="Amplitude (\\muV)", spectraArgs...)
plot(𝐒[2]; ytitle="Amplitude (\\muV)", spectraArgs...)

# do the same thing using a fast FFTW plan (wait up to 10s for computing the plan)
plan=Planner(plan_exhaustive, 10.0, t, eltype(𝐗[1]))
𝐒_fast=spectra(𝐗, sr, t; planner=plan, slide=sr, func=√)

# compute the average spectra in the alpha range for all time-series and all subjects
mean(𝐒, (8,12))

# compute the average spectra in the alpha range across subjects for all time-series
mean(mean(𝐒, (8,12)))
# compute and plot it
plot(mean(mean(𝐒, (8,12))))

# extract spectra in alpha range for all time-series and all subjects
extract(𝐒, (8, 12))

# if you enter en illegal range, nothing will be done and you will get
# an error in the REPL explaining what is wrong in your argument
extract(𝐒, (0, 12))
mean(𝐒, (0, 128))

# extract 4Hz-band-pass average spectra for all electrodes and all subjects
bands(𝐒, 4)

# Apply smoothing in the above spectra computations
𝐒=spectra(𝐗, sr, t;
    tapering=tapering, smoothing=blackmanSmoother, func=√)
plot(𝐒[1]; ytitle="Amplitude (\\muV)", spectraArgs...)
plot(𝐒[2]; ytitle="Amplitude (\\muV)", spectraArgs...)

# plot the average spectrum across all electrodes for the two files
# using Julia standard mean function
plot(mean(𝐒[1].y[:, i] for i=1:size(𝐒[1].y, 2)))
plot!(mean(𝐒[2].y[:, i] for i=1:size(𝐒[2].y, 2)))

# plot spectra in in 1Hz band-pass regions for all electrodes in 𝐒[1]
plot(bands(𝐒[1], 1))

# use slepian multi-tapering on EEG data
S_sl=spectra(X1, sr, t; tapering=slepians(sr, t, 1.25), func=√)
plot(S_sl; ytitle="Amplitude (\\muV)", spectraArgs...)

# use slepian multi-tapering and get smoothed spectra
S_sl=spectra(X1, sr, t;
    tapering=slepians(sr, t, 1.25), func=√, smoothing=blackmanSmoother)
plot(S_sl; xspace=8, ytitle="Amplitude (\\muV)", spectraArgs...)

# See how the variance of the spectra estimation is lower
# as compared to using a normal tapering window !
S=spectra(X1, sr, t;
    tapering=harris4, func=√, smoothing=blackmanSmoother)
plot(S; ytitle="Amplitude (\\muV)", spectraArgs...)


text:
using Test
using Spandex

@testset "amd 1" begin
    local g = DirectedGraph{Int64}(3)
    g[1, 1] = 6
    g[1, 2] = -4
    g[1, 3] = 7
    g[2, 1] = 4
    g[2, 2] = -3
    g[2, 3] = 0
    g[3, 1] = 1
    g[3, 2] = 2
    g[3, 3] = 5
    local a = from_graph(g, 3, 3)

    local pt = PermuteTable(a)

    @test 3 == length(pt.permuted)
    @test 3 == pt.permuted[1]
    @test 2 == pt.permuted[2]
    @test 1 == pt.permuted[3]

    @test 3 == length(pt.primary)
    @test 3 == pt.primary[1]
    @test 2 == pt.primary[2]
    @test 1 == pt.primary[3]
end

@testset "amd 2" begin
    local g = DirectedGraph{Float64}(5)
    g[1, 1] = 0.44561
    g[1, 2] = 0.92962
    g[4, 2] = 0.87043
    g[5, 2] = 0.47918
    g[1, 3] = 0.10582
    g[4, 3] = 0.060281
    g[1, 4] = 0.87256
    g[2, 4] = 0.41421
    g[5, 4] = 0.42161
    g[3, 5] = 0.58921
    g[5, 5] = 0.15305
    g[4, 6] = 0.63730
    local a = from_graph(g, 5, 6)

    local pt = PermuteTable(a)

    @test 5 == pt.permuted[1] || 6 == pt.permuted[1]
    @test 5 == pt.permuted[2] || 6 == pt.permuted[2]

    local uniques = length(unique(pt.permuted))
    @test 6 == uniques

    uniques = length(unique(pt.primary))
    @test 6 == uniques
end