text:
function make_dict2(args,x)
    nargs = 0
    for i in [1:length(args);]
        if typeof(args[i])==String # if the entry in args is a string, then there's one corresponding scalar entry in x0
            nargs += 1
        else
            nargs += args[i][2]    # otherwise, the entry in args should be a  [varnamestring, nvals] vector,
            # indicating that the next nvals entries in x0 are all a single vector, belonging to variable
            # with name varnamestring.
        end
    end
    if nargs != length(x)
        error("Oy! args and x must indicate the same total number of variables!")
    end

    i=1
    param=Dict()
    for  (iarg,arg) in enumerate(args)
        if typeof(arg)==String
            param[arg]=x[i:i]
            i+=1
        else
            j=arg[2]
            param[arg[1]]=x[i:i+j-1]
            i+=j
        end
    end

    return param
end

function arg_labels(args)
    labels=[]
    for arg in args
        if arg[2]!=0
            if length(arg[2])==1
                append!(labels,[arg[1]])
            else
                for ip in arg[2]
                    append!(labels, [ arg[1]*string(ip) ])
                end
            end
        end
    end

    return labels
end

function make_dict(xargs,x,zargs,z)
    """
    Function to combine paramters to fit (args,x)
    and constants (zargs,z)
    """

    args=[]
    for iarg in 1:length(xargs)
        if xargs[iarg][2]==0
            ii=length(zargs[iarg][2])
            if ii==1
                push!( args,xargs[iarg][1] )
            else
                push!( args,( xargs[iarg][1],ii) )
            end

        elseif zargs[iarg][2]==0
            ii=length(xargs[iarg][2])
            if ii==1
                push!( args,xargs[iarg][1] )
            else
                push!( args,( xargs[iarg][1],ii) )
            end

        else
            push!( args,(xargs[iarg][1],length(zargs[iarg][2])+length(xargs[iarg][2])) )
        end

    end

    y=zeros(typeof(x[1]),length(x)+length(z))
    ix=1
    iz=1
    iy=1
    for iarg in 1:length(xargs)

        if xargs[iarg][2]==0

            ii=length(zargs[iarg][2])
            # for ind in 1:ii
            #     push!(y,z[iz:iz+ind-1])
            # end

            y[iy:iy+ii-1]=z[iz:iz+ii-1]
            iz+=ii
            iy+=ii

        elseif zargs[iarg][2]==0

            ii=length(xargs[iarg][2])

            # for ind in 1:ii
            #     push!(y,z[iz:iz+ind-1])
            # end

            y[iy:iy+ii-1]=x[ix:ix+ii-1]
            ix+=ii
            iy+=ii
        else
            ii=length(xargs[iarg][2])+length(zargs[iarg][2])
            iix=length(xargs[iarg][2])
            iiz=length(zargs[iarg][2])

            y[iy+xargs[iarg][2]-1]=x[ix:ix+iix-1]
            y[iy+zargs[iarg][2]-1]=z[iz:iz+iiz-1]

            ix+=iix
            iy+=ii
            iz+=iiz
        end

    end

    nargs = 0
    for i in [1:length(args);]
        if typeof(args[i])==String # if the entry in args is a string, then there's one corresponding scalar entry in x0
            nargs += 1
        else
            nargs += args[i][2]    # otherwise, the entry in args should be a  [varnamestring, nvals] vector,
            # indicating that the next nvals entries in x0 are all a single vector, belonging to variable
            # with name varnamestring.
        end
    end
    if nargs != length(y)
        error("Oy! args and x must indicate the same total number of variables!")
    end

    i=1
    param=Dict()
    for  (iarg,arg) in enumerate(args)
        if typeof(arg)==String
            param[arg]=y[i:i]
            i+=1
        else
            j=arg[2]
            param[arg[1]]=y[i:i+j-1]
            i+=j
        end
    end

    return param
end

xargs=[ ("c2",[1,2]),("c4",0),("st",[2,3]),("x0",[1,2]),("bias_d",[1]),("hbias_d",[1]),("sigma_a",[1])]
zargs=[("c2",0),("c4",[1]),("st",[1,4]),("x0",0),("bias_d",0),("hbias_d",0),("sigma_a",0)]

x=[1,2,5,6,8,9,10,11,12]
z=[3,4,7]

args=[]
for iarg in 1:length(xargs)
    if xargs[iarg][2]==0
        ii=length(zargs[iarg][2])
        if ii==1
            push!( args,xargs[iarg][1] )
        else
            push!( args,( xargs[iarg][1],ii) )
        end

    elseif zargs[iarg][2]==0
        ii=length(xargs[iarg][2])
        if ii==1
            push!( args,xargs[iarg][1] )
        else
            push!( args,( xargs[iarg][1],ii) )
        end

    else
        push!( args,(xargs[iarg][1],length(zargs[iarg][2])+length(xargs[iarg][2])) )
    end

end

function get_model_f(model,reg; simulations=false)
    global model_f

    if (model=="StLogisticAdaptatioAbDwX0MubReg"*string(reg) )

        model_f=Dict(:c2_f=>c2_const,:c4_f=>c4_const,:c6_f=>c6_zeros,:bias_d_f=>bias_d_const,
        :hbias_d_f=>history_bias_d_const,:st_f=>st_logistic_adaptation_bias,:x0_f=>initial_condition_bias_hbias,
        :ro_f=>read_out_perfect,:lambda=>reg)

    elseif (model=="StLinearPiReg"*string(reg) )
        model_f=Dict(:bias_d_f=>bias_d_const,:hbias_d_f=>history_bias_d_const,
        :st_f=>st_linear,:x0_f=>initial_condition_bias_hbias,
        :ro_f=>read_out_perfect,:lambda=>reg)

    elseif (model=="StLogisticPiReg"*string(reg) )

        model_f=Dict(:bias_d_f=>bias_d_const,:hbias_d_f=>history_bias_d_const,
        :st_f=>st_logistic,:x0_f=>initial_condition_bias_hbias,
        :ro_f=>read_out_perfect,:lambda=>reg)

    elseif (model=="StLogisticAdaptationDwReg"*string(reg))
        model_f=Dict(:c2_f=>c2_const,:c4_f=>c4_const,:c6_f=>c6_zeros,:bias_d_f=>bias_d_const,
        :hbias_d_f=>history_bias_d_const,:st_f=>st_logistic_adaptation,:x0_f=>initial_condition_bias_hbias,
        :ro_f=>read_out_perfect,:lambda=>reg)

    elseif (model=="StLinearAdaptationExpDwReg"*string(reg))
        model_f=Dict(:c2_f=>c2_const,:c4_f=>c4_const,:c6_f=>c6_zeros,:bias_d_f=>bias_d_const,
        :hbias_d_f=>history_bias_d_const,:st_f=>st_linear_adaptation_exp,:x0_f=>initial_condition_bias_hbias,
        :ro_f=>read_out_perfect,:lambda=>reg)

    elseif (model=="StLinearDwReg"*string(reg))
        model_f=Dict(:c2_f=>c2_const,:c4_f=>c4_const,:c6_f=>c6_zeros,:bias_d_f=>bias_d_const,
        :hbias_d_f=>history_bias_d_const,:st_f=>st_linear,:x0_f=>initial_condition_bias_hbias,
        :ro_f=>read_out_perfect,:lambda=>reg)

    elseif (model=="StLogisticDwReg"*string(reg))
        model_f=Dict(:c2_f=>c2_const,:c4_f=>c4_const,:c6_f=>c6_zeros,:bias_d_f=>bias_d_const,
        :hbias_d_f=>history_bias_d_const,:st_f=>st_logistic,:x0_f=>initial_condition_bias_hbias,
        :ro_f=>read_out_perfect,:lambda=>reg)

    elseif (model=="StLogisticAdaptationC4Reg"*string(reg))
        model_f=Dict(:c2_f=>c2_const,:c4_f=>c4_const,:c6_f=>c6_zeros,:bias_d_f=>bias_d_const,
        :hbias_d_f=>history_bias_d_const,:st_f=>st_logistic_adaptation,:x0_f=>initial_condition_bias_hbias,
        :ro_f=>read_out_perfect,:lambda=>reg)

    elseif (model=="StLogisticC4Reg"*string(reg))
        model_f=Dict(:c2_f=>c2_const,:c4_f=>c4_const,:c6_f=>c6_zeros,:bias_d_f=>bias_d_const,
        :hbias_d_f=>history_bias_d_const,:st_f=>st_logistic,:x0_f=>initial_condition_bias_hbias,
        :ro_f=>read_out_perfect,:lambda=>reg)

    else:
        println("not found")
    end
    if simulations
        model_f[:ro_f]=read_out_perfect_sim
    end
    return model_f

end

function get_args_consts(model,reg)
    args=[]
    consts=[]
    lb=[]
    ub=[]
    z=[]
    if (model=="StLogisticAdaptatioDwReg"*string(reg))
        args=[ ("c2",[1]),("c4",[1]),("st",[1,2,3]),("x0",0),("bias_d",0),("hbias_d",0),("sigma_a",0)]
        consts=[("c2",0),("c4",0),("st",0),("x0",[1,2]),("bias_d",[1]),
        ("hbias_d",[1]),("sigma_a",[1])]

        z=1.0*[0.0,0.0,0.0,0.0,1.0]
        lb=[-50., 0.0,  0.0,0.0,-50.0 ]
        ub=[50., 50.,  50., 50.0,  50.0]

    elseif (model=="StLogisticAdaptationC4Reg"*string(reg))
        args=[ ("c2",0),("c4",[1]),("st",[1,2,3]),("x0",0),("bias_d",0),("hbias_d",0),("sigma_a",0)]
        consts=[("c2",[1]),("c4",0),("st",0),("x0",[1,2]),("bias_d",[1]),
        ("hbias_d",[1]),("sigma_a",[1])]

        z=1.0*[0.0,0.0,0.0,0.0,0.0,1.0]
        lb=[0., 0.0,  0.0,-50.0]
        ub=[50., 50.,  50., 50.0]

    elseif (model=="StLogisticC4Reg"*string(lambda))
        args=[ ("c2",0),("c4",[1]),("st",[1,2]),("x0",0),("bias_d",0),("hbias_d",0),("sigma_a",0)]
        consts=[("c2",[1]),("c4",0),("st",0),("x0",[1,2]),("bias_d",[1]),
        ("hbias_d",[1]),("sigma_a",[1])]

        z=1.0*[0.0,0.0,0.0,0.0,0.0,1.0]
        lb=[0., 0.0,  0.0]
        ub=[50., 50.,  50.]

    else
        println("model not found")
    end

    return args,consts,z,lb,ub
end

function create_stimuli(mu,sigmas,N,Nframes)
    m=mu*ones(N)
    m[Int(floor(N/2)):end]=-m[Int(floor(N/2)):end]
    stim=randn(N,Nframes)
    sigmas_vec=[]
    i_sigmas=[]
    for i in 1:Int(N/length(sigmas))
        sigmas_vec=vcat(sigmas_vec,sigmas)
        i_sigmas=vcat(i_sigmas,0:length(sigmas)-1)
    end

    for i in 1:N
        stim[i,:]=m[i]+sigmas_vec[i]*(stim[i,:] -mean(stim[i,:]) )/std(stim[i,:])
    end
    return stim,convert(Array{Float64,1},sigmas_vec)[1:N],convert(Array{Float64,1},i_sigmas)[1:N],sign.(m)

end

max_stars_repo_name: emmafdz/FinalDegreeProject

max_issues_repo_name: emmafdz/FinalDegreeProject

max_forks_repo_name: emmafdz/FinalDegreeProject


text:
export
        ConvexPolygon,
        CPAMemory,
        CollisionCheckResult,

        LineSegment,

        to_oriented_bounding_box!,
        get_oriented_bounding_box,
        is_colliding,
        is_potentially_colliding,
        get_collision_time,
        get_first_collision,
        get_time_and_dist_of_closest_approach,
        is_collision_free,
        get_distance,
        get_edge,
        get_side,

        CollisionCallback

include("minkowski.jl")

ext: jl


text:
import Compat: readstring

# Patches should allow using imported bindings in the body of the patch
@test_throws UndefVarError Minute
@test isdefined(Base.Dates, :Minute)
import Base.Dates: Minute, Hour

let
    myminute(x::Integer) = Minute(x)

    # Patches should work when referencing bindings imported in the file where the patch
    # is created.
    patch = @patch myminute(x::Integer) = Minute(Hour(x))
    apply(patch) do
        @test (@mock myminute(5)) == Minute(300)
    end
end

# Patches should allow using <module>.<binding> syntax in the signature
@test_throws UndefVarError AbstractCmd
@test isdefined(Base, :AbstractCmd)

patch = @patch readstring(cmd::Base.AbstractCmd) = "bar"
apply(patch) do
    @test (@mock readstring(`foo`)) == "bar"
end

# Patches should allow using imported bindings syntax in the signature
import Base: AbstractCmd

patch = @patch readstring(cmd::AbstractCmd) = "bar"
apply(patch) do
    @test (@mock readstring(`foo`)) == "bar"
end


text:
##### structs.jl #####

# In this file you will find all the structures needed to generalize code.
# You will find also some functions to help create the global structures.

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

export BinaryGene, IntegerGene, FloatGene
export Crossover, Selection
export selection, crossover, bin
export GAExternal, ext_init, distributed_ga

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

"""
    BinaryGene(value ::Bool, name ::AbstractString)

Creates a `BinaryGene` structure. This gene represents a `Bool` variable (1 or 0). Useful to make decisions (will a reactor work or not, will a SMB pathway be used or not, etc.). `name` is the variable name for presentation purposes.
"""
mutable struct BinaryGene <: AbstractGene
    value ::Bool
    name  ::AbstractString

    function BinaryGene(value ::Bool, name ::AbstractString)
        return new(value, name)
    end
end

"""
    BinaryGene()

Creates a `BinaryGene` structure with a random `Bool` value.
"""
BinaryGene(name ::AbstractString) = BinaryGene(rand(Bool), name)

"""
    BinaryGene()

Creates a `BinaryGene` structure with a random `Bool` value and a default variable name.
"""
BinaryGene() = BinaryGene(rand(Bool), "bin")

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

"""
    IntegerGene(value ::BitVector, name ::AbstractString)

Creates a `IntegerGene` structure. This gene represents an integer variable as `BitVector`. To convert the `BitVector` in an integer, just look at the `bin` function from this package by typing `?bin` on the command prompt. `name` is a string that represents the name of the variable. It's needed for result presentation purposes.
"""
mutable struct IntegerGene <: AbstractGene
    value  ::BitVector
    lbound ::Real
    ubound ::Real
    name   ::AbstractString
end

"""
    IntegerGene(mutation ::Symbol)

Chooses the mutation type. Does NOT create a structure. This dispatch was created because, when creating a population of integer genes, each gene would create this function, which was unnecessary work. You still need to run one of the other dispatches to create the gene.

Below are the types of mutations supported:

| Symbol | Algorithm          |
|--------|--------------------|
| :FM    | Flip Mutation      |
| :InvM  | Inversion Mutation |
| :InsM  | Insertion Mutation |
| :SwM   | Swap Mutation      |
| :ScrM  | Scramble Mutation  |
| :ShM   | Shifting Mutation  |
"""
function IntegerGene(mutation ::Symbol)
    int_func = mutate(mutation)
    @eval mutate(gene ::IntegerGene) = $int_func(gene.value)
    return nothing
end

"""
    IntegerGene(n ::Int64)

Creates a `IntegerGene` structure in which the `BitVector` is of length `n` with random `Bool` values.
"""
function IntegerGene(n ::Int64, name ::AbstractString;
                     lb ::Real = -Inf, ub ::Real = Inf )
    value  = BitVector(undef, n)
    for i in 1:n
        @inbounds value[i] = rand(Bool)
    end
    return IntegerGene(value, lb, ub, name)
end

"""
    bin(gene ::IntegerGene)

Returns the integer number represented by the `BitVector` of `gene`.
"""
function bin(gene ::IntegerGene)
    bin_val = 0
    for (i,j) in enumerate(gene.value)
        bin_val += j ? 2^(i-1) : 0
    end
    return bin_val
end

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

"""
    FloatGene(value ::Vector{Float64}, range ::Vector{Float64}, m ::Int64)

Creates a `FloatGene` structure. `value` is a vector with the variables to be changed. `range` is a vector with the minimum and maximum values a variable can take. `m` is just a parameter that changes how much in a mutation the variables change, the bigger the value, the bigger the change in each mutation. If the range of a variable is 0.5, then the biggest mutation a variable can suffer in one iteration is 0.5, for instance.
"""
mutable struct FloatGene <: AbstractGene
    value  ::Vector{Float64}
    range  ::Vector{Float64}
    name   ::Vector{AbstractString}
    m      ::Int64
    lbound ::Vector{<:Real}
    ubound ::Vector{<:Real}

    function FloatGene( value  ::Vector{Float64}           ,
                        range  ::Vector{Float64}           ,
                        name   ::Vector{String}            ,
                        m      ::Int64                     ,
                        lbound ::Union{Vector{<:Real},Real},
                        ubound ::Union{Vector{<:Real},Real})
        if length(value) != length(range)
            error("value and range have different lengths")
        end
        if typeof(lbound) <: Real
            lb = Float64[lbound for i in value]
        else
            if length(value) != length(lbound)
                error("value and lbound have different lengths")
            else
                lb = lbound
            end
        end
        if typeof(ubound) <: Real
            ub = Float64[ubound for i in value]
        else
            if length(value) != length(ubound)
                error("value and ubound have different lengths")
            else
                ub = ubound
            end
        end
        return new(value, range, name, m, lb, ub)
    end
end

function FloatGene( value ::Vector{Float64}                  ,
                    range ::Vector{Float64}                  ,
                    name  ::Vector{String}                   ;
                    m     ::Int64                      =  20 ,
                    lb    ::Union{Vector{<:Real},Real} = -Inf,
                    ub    ::Union{Vector{<:Real},Real} =  Inf)

    return FloatGene(value, range, name, m, lb, ub)
end

"""
    FloatGene(value ::Float64, range ::Float64; m ::Int64 = 20)

Creates a `FloatGene` structure. Handy for creating just one real number variable.
"""
function FloatGene( value ::Float64                          ,
                    range ::Float64                          ,
                    name  ::AbstractString                   ;
                    m     ::Int64                      =  20 ,
                    lb    ::Union{Vector{<:Real},Real} = -Inf,
                    ub    ::Union{Vector{<:Real},Real} =  Inf)
    return FloatGene([value], [range], [name], m, lb, ub)
end

"""
    FloatGene(value ::Vector{Float64}, range ::Float64; m ::Int64 = 20)

Creates a `FloatGene` structure. Handy for creating a vector of real numbers with the same range.
"""
function FloatGene( value ::Vector{Float64}                  ,
                    range ::Float64                          ,
                    name  ::Vector{String}                   ;
                    m     ::Int64                      =  20 ,
                    lb    ::Union{Vector{<:Real},Real} = -Inf,
                    ub    ::Union{Vector{<:Real},Real} =  Inf)
    range_vec = Float64[range for i in value]
    return FloatGene(value, range_vec, name, m, lb, ub)
end

"""
    FloatGene(n ::Int64)

Creates a `FloatGene` structure. Creates a vector of length `n` with random variables and random ranges. Used particularly for testing purposes.
"""
function FloatGene(n ::Int64, name ::AbstractString; m ::Int64 = 20)
    value = rand(Float64, n)
    range = rand(Float64, n)
    lb    = rand(Float64, n)
    ub    = rand(Float64, n)
    vec_name = Vector{String}(undef, n)
    for i in 1:n
        vec_name[i] = string(name, i)
    end
    return FloatGene(value, range, vec_name, m, lb, ub)
end

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

function isbound(gene ::FloatGene)
    lb = findmin( gene.value  .- gene.lbound )[1]
    ub = findmin( gene.ubound .- gene.value  )[1]
    if lb < 0.0
        @info("lower boundary violated, retrying...")
        return false
    end
    if ub < 0.0
        @info("upper boundary violated, retrying...")
        return false
    end
    return true
end

function isbound(gene ::IntegerGene)
    return bin(gene) >= gene.lbound && bin(gene) <= gene.ubound
end

isbound(gene ::BinaryGene) = true

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

"""
    Crossover(cross ::Symbol                                    ;
              w     ::Union{Nothing, Vector{Float64}} = nothing ,
              d     ::Union{Nothing, Float64        } = nothing )

Creates a `Crossover` structure. `cross` is a Symbol that represents the type of crossover that would be used. `w` and `d` are not mandatory but need to be set for some types of crossovers. All algorithms will be shown in the table below:

| Symbol | Algorithm                            | Optional Arguments |
|--------|--------------------------------------|--------------------|
| :SPX   | Single Point Crossover               | not needed         |
| :TPX   | Two Point Crossover                  | not needed         |
| :UX    | Uniform Crossover                    | not needed         |
| :DX    | Discrete Crossover                   | not needed         |
| :WMX   | Weighted Mean Crossover              | needs `w`          |
| :IRX   | Intermediate Recombination Crossover | needs `d`          |
| :LRX   | Line Recombination Crossover         | needs `d`          |
| :PMX   | Partially Mapped Crossover           | not needed         |
| :O1X   | Order 1 Crossover                    | not needed         |
| :O2X   | Order 2 Crossover                    | not needed         |
| :CX    | Cycle Crossover                      | not needed         |
| :PX    | Position-based Crossover             | not needed         |
"""
mutable struct Crossover
    cross ::Symbol
    w     ::Union{Nothing, Vector{Float64}}
    d     ::Union{Nothing, Float64        }

    function Crossover(cross ::Symbol                                    ;
                       w     ::Union{Nothing, Vector{Float64}} = nothing ,
                       d     ::Union{Nothing, Float64        } = nothing )
            cross_func = nothing
            if cross == :SPX
                cross_func = singlepoint
            elseif cross == :TPX
                cross_func = twopoint
            elseif cross == :UX
                cross_func = uniform
            elseif cross == :DX
                cross_func = discrete
            elseif cross == :WMX
                if isnothing(w)
                    error("value `w` must be given a value")
                end
                cross_func = waverage(w)
            elseif cross == :IRX
                if isnothing(d)
                    error("value `d` must be given a value")
                end
                cross_func = intermediate(d)
            elseif cross == :LRX
                if isnothing(d)
                    error("value `d` must be given a value")
                end
                cross_func = line(d)
            elseif cross == :PMX
                cross_func = pmx
            elseif cross == :O1X
                cross_func = ox1
            elseif cross == :O2X
                cross_func = ox2
            elseif cross == :CX
                cross_func = cx
            elseif cross == :PX
                cross_func = pos
            end

            @eval begin
                function crossover(v1 ::T, v2 ::T) where {T <: AbstractVector}
                    return $cross_func(v1, v2)
                end
            end
        return new(cross, w, d)
    end
end

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

"""
    Selection( select    ::Symbol                            ;
               sp        ::Union{Nothing, Float64} = nothing ,
               μ         ::Union{Nothing,   Int64} = nothing ,
               groupsize ::Union{Nothing,   Int64} = nothing )

Creates a `Selection` structure. `select` is a symbol that represents the type of selection that will be used. `sp`, `μ` and `groupsize` are optional but need to be set for some types of selections. All algorithms will be shown in the table below:

| Symbol | Algorithm                               | Optional Arguments |
|--------|-----------------------------------------|--------------------|
| :RBS   | Rank-based Selection                    | needs `sp`         |
| :URS   | Uniform-Ranking Selection               | needs `μ`          |
| :RWS   | Roulette Wheel Selection                | not needed         |
| :SUSS  | Stochastic Universal Sampling Selection | not needed         |
| :TrS   | Truncation Selection                    | not needed         |
| :ToS   | Tournament Selection                    | needs `groupsize`  |
"""
mutable struct Selection
    select ::Symbol
    sp     ::Union{Nothing, Float64}
    μ      ::Union{Nothing,   Int64}
    gsize  ::Union{Nothing,   Int64}

    function Selection( select    ::Symbol                            ;
                        sp        ::Union{Nothing, Float64} = nothing ,
                        μ         ::Union{Nothing,   Int64} = nothing ,
                        groupsize ::Union{Nothing,   Int64} = nothing )
        selec_func = nothing
        if select == :RBS
            if isnothing(sp)
                error("need to specify `sp` value")
            end
            selec_func = ranklinear(sp)
        elseif select == :URS
            if isnothing(μ)
                error("need to specify `μ` value")
            end
            selec_func = uniformranking(μ)
        elseif select == :RWS
            selec_func = roulette
        elseif select == :SUSS
            selec_func = sus
        elseif select == :TrS
            selec_func = truncation
        elseif select == :ToS
            if isnothing(groupsize)
                error("need to specify `groupsize` value")
            end
            selec_func = tournament(groupsize)
        else
            error("Unknown parameter " * string(select))
        end

        @eval begin
            function selection(fit ::Vector{<:Real}, N ::Int)
                return $selec_func(fit, N)
            end
        end
        return new(select, sp, μ, groupsize)
    end
end

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

"""
    function GAExternal( program  ::AbstractString                   ,
                         pipename ::AbstractString                   ;
                         nworkers ::Int64          = Sys.CPU_THREADS ,
                         parallel ::Bool           = false           )

Creates communication pipes for the external program `program`. If `parallel` is `true`, then, considering N workers available, N pipes for reading and N pipes for writing will be created. `pipename` is just a handle for the name of the pipes. If `pipename` is `pipe`, then the pipe names will be `pipe_in` and `pipe_out` if `parallel=false` and `pipe_inn` and `pipe_outn` if `parallel=true`, with `n` being one of the N workers. `nworkers` is the number of cores to be used, including the number of cores of a remote computer. `parallel` sets the the external program to run in several workers.
"""
mutable struct GAExternal
    program       ::AbstractString
    pipes_in      ::Union{Vector{<:AbstractString},DArray{String,1,Vector{String}}}
    pipes_out     ::Union{Vector{<:AbstractString},DArray{String,1,Vector{String}}}
    avail_workers ::Vector{Int64}
    parallel      ::Bool

    function GAExternal( program  ::AbstractString                   ,
                         pipename ::AbstractString                   ;
                         nworkers ::Int64          = Sys.CPU_THREADS ,
                         parallel ::Bool           = false           )
        pipes = Dict{String,Vector{String}}()
        pipes["in" ] = Vector{String}(undef, 0)
        pipes["out"] = Vector{String}(undef, 0)
        avail_workers = workers()[1:nworkers]
        if parallel
            # create one pipe for reading and another for writing
            # for each worker
            for i in ["in","out"]
                for p in avail_workers
                    f = string(pipename, "_", i, p)
                    push!(pipes[i], f)
                    remotecall_fetch(rm, p, f, force=true)
                    remotecall_fetch(run, p, `mkfifo $f`)
                end
            end
            pin  = distribute(pipes["in" ]; procs=avail_workers)
            pout = distribute(pipes["out"]; procs=avail_workers)
        else
            # create one pipe for reading and another for writing
            for i in ["in","out"]
                f = string(pipename, "_", i)
                push!(pipes[i], f)
                rm(f, force=true)
                run(`mkfifo $f`)
            end
            pin = pipes["in"]
            pout = pipes["out"]
        end

        #pin  = distribute(pipes["in" ]; procs=avail_workers)
        #pout = distribute(pipes["out"]; procs=avail_workers)

        # activate writing pipes for a big amount of time
        for (i,p) in enumerate(pin)
            id  = 1*nworkers+1 + i
            id1 = 2*nworkers+1 + i
            @spawnat id  run(pipeline(`sleep 100000000`; stdout=p))
            @spawnat id1 run(pipeline(`$program`; stdin=p))
        end

        # Open reading pipes in separate processes.
        # The pipes have to be open before writing to them,
        # otherwise we get SIGPIPE errors when writing.
        function spawn_readpipes(pipe)
            f = open(pipe, "r")
            return nothing
        end
        for (i,p) in enumerate(pout)
            v = 3*nworkers+1 + i
            @spawnat v spawn_readpipes(p)
        end

        # delete all pipes when exiting julia
        function external_atexit()
            for k in keys(pipes)
                for (i,p) in enumerate(pipes[k])
                    id = i+1
                    remotecall_fetch(rm, id, p)
                end
            end
            return nothing
        end
        atexit(external_atexit)

        return new(program, pin, pout, avail_workers, parallel)
    end
end

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

"""
    ext_init(gaext ::GAExternal, codeline ::AbstractString)

Function to help initialize the external program.

When using external programs, you usually have an initial script to be run before performing the genetic algorithm. After creating the `GAExternal` structure, just send it to this function along with the line of code to initialize the optimization process.

## Example

In case of AMPL, if we have a script with the initial model called `init.ampl`, then we could do:

```
ext = GAExternal("ampl", "pipe")
ext_init(ext, "include init.ampl;")
```
"""
function ext_init(gaext ::GAExternal, codeline ::AbstractString)
    function ext(gaext ::GAExternal)
        if gaext.parallel
            pin = gaext.pipes_in[:L]
        else
            pin = gaext.pipes_in
        end
        for p in pin
            open(p, "w") do file
                write(file, codeline)
            end
        end
    end
    if gaext.parallel
        spmd(ext, gaext, pids=gaext.avail_workers)
    else
        ext(gaext)
    end
    return nothing
end

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

"""
    distributed_ga( ;
                    localcpu ::Int64           = Sys.CPU_THREADS ,
                    cluster  ::Vector{<:Tuple} = [()]            ,
                    dir      ::AbstractString  = pwd()           )

Function to help set up the local computer or a cluster for parallel run of the Genetic Algorithm. `localcpu` is the number of cores to be used in your local computer. `cluster` is a vector of machine specifications. To know more about this, type `?addprocs` in the command prompt after importing the `Distributed` package. `dir` is the directory where you want julia to run in each remote computer.
"""
function distributed_ga( ;
                         localcpu ::Int64           = Sys.CPU_THREADS ,
                         cluster  ::Vector{<:Tuple} = [()]            ,
                         dir      ::AbstractString  = pwd()           )

    nworkers = localcpu
    if cluster != [()]
        for i in cluster
            nworkers += i[2]
        end
    end

    @eval using Distributed
    for i in 1:4
        if localcpu != 0
            addprocs(localcpu)
        end
        if cluster != [()]
            addprocs(cluster, dir=dir)
        end
    end

    @eval @everywhere using Distributed
    @eval @everywhere using DistributedArrays
    @eval @everywhere using DistributedArrays.SPMD
    @eval @everywhere using Evolutionary

    return nworkers
end


text:
# # Grover Search
using Yao
using LinearAlgebra
using QuAlgorithmZoo: groveriter, inference_oracle, prob_match_oracle

# ## Target Space and Evidense
num_bit = 12
oracle = matblock(Diagonal((v = ones(1<<num_bit); v[100:101]*=-1; v)))
target_state = zeros(1<<num_bit); target_state[100:101] .= sqrt(0.5)

# now we want to search the subspace with [1,3,5,8,9,11,12] fixed to 1 and [4,6] fixed to 0.
evidense = [1, 3, -4, 5, -6, 8, 9, 11, 12]

# ## Search
# then solve the above problem
it = groveriter(uniform_state(num_bit), oracle)
for (i, psi) in enumerate(it)
    overlap = abs(statevec(psi)'*target_state)
    println("step $(i-1), overlap = $overlap")
end

# ## Inference Example
# we have a state psi0, we know how to prepair it
psi = rand_state(num_bit)

"""
Doing Inference, psi is the initial state, we want to search target space with specific evidense.
e.g. evidense [1, -3, 6] means the [1, 3, 6]-th bits take value [1, 0, 1].
"""
oracle_infer = inference_oracle(evidense)(nqubits(psi))
it = groveriter(psi, oracle_infer)
for (i, psi) in enumerate(it)
    p_target = prob_match_oracle(psi, oracle_infer)
    println("step $(i-1), overlap^2 = $p_target")
end


text:
using Test
using WavePropBase
using WavePropBase.Geometry
using WavePropBase.Integration
using WavePropBase.Interpolation
using WavePropBase.Mesh
using LinearAlgebra
using StaticArrays

@testset "Lagrange elements" begin
    @testset "LagrangeLine" begin
        ## line in 3d
        vtx = ((0.,0.,0.),(1.,1.,1.))
        l   = LagrangeLine(vtx)
        @test domain(l) == ReferenceLine()
        @test length(return_type(l)) == 3
        @test dimension(domain(l)) == 1
        # makes no sense to ask for normal here. Make sure error is thrown.
        @test_throws AssertionError normal(l,SVector(0.5)) == SVector(-1.,1.)/√2
        ## line in 2d
        a = (0.,0.)
        b = (1.,1.)
        l   = LagrangeLine((a,b))
        @test domain(l) == ReferenceLine()
        @test length(return_type(l)) == 2
        @test dimension(domain(l)) == 1
        @test normal(l,SVector(0.5)) == SVector(1.,-1.)/√2
    end
    @testset "LagrangeTriangle" begin
        # triangle in 2d
        vtx = SVector(SVector(0.,0.),SVector(0.,1.),SVector(-1.,0))
        t   = LagrangeTriangle(vtx)
        @test length(return_type(t)) == 2
        @test dimension(domain(t)) == 2
        # triangle in 3d
        vtx = SVector(SVector(0.,0.,0.),SVector(0.,1.,0.),SVector(-1.,0,0.))
        t   = LagrangeTriangle(vtx)
        @test length(return_type(t)) == 3
        @test dimension(domain(t)) == 2
        @test normal(t,SVector(0.1,0.1)) == SVector(0,0,1.)
    end
    @testset "Tetrahedron" begin
        # TODO: add tests
    end
end

max_stars_repo_name: WaveProp/WavePropBase

max_issues_repo_name: WaveProp/WavePropBase

max_forks_repo_name: WaveProp/WavePropBase


text:
#####
##### Parse `dnnl_blocking_desc_t`
#####

# Define a long dispatch chain to unwrap types until we finally get to the final
# blocking description.

layout(x::Ref{Lib.dnnl_memory_desc_t}) = layout(unwrap_ref(x))
function layout(memory_desc::Lib.dnnl_memory_desc_t)
    # Get the blocking desription.
    # We use this in conjunction with the `memory_desc_t` to reconstruct the format.
    blocking_desc = _blocking_desc(memory_desc)
    ndims = memory_desc.ndims

    # Construct the inner tilings
    tiles = TiledArrays.Tile[]
    inner_num_blocks = blocking_desc.inner_nblks
    inner_blks = reverse(blocking_desc.inner_blks[1:inner_num_blocks])
    #inner_idxs = blocking_desc.inner_idxs[1:inner_num_blocks]
    inner_idxs = ndims .- reverse(blocking_desc.inner_idxs[1:inner_num_blocks])
    for i in 1:inner_num_blocks
        # Convert from index-0 to index-1.
        push!(tiles, TiledArrays.Tile(inner_idxs[i], inner_blks[i]))
    end
    # Order the indices.
    _strides = collect(reverse(blocking_desc.strides[1:ndims]))
    outer = sortperm(_strides)
    return (tiles..., outer...)
end

function _blocking_desc(x::Lib.dnnl_memory_desc_t)
    if x.format_kind != Lib.dnnl_blocked
        @show x.format_kind
        error("Wrong Format Kind!")
    end
    return x.format_desc
    #return _blocking_desc(x.format_desc)
end

# TODO: Fix this weird "ANONYMOUS1" thing caused by the Union ...
#_blocking_desc(x::Lib.__JL_Ctag_89) = x.blocking

function generate_linear_indices(layout::TiledArrays.MemoryLayout, size, padded_size)
    indexer = TiledArrays.TiledIndexer{layout}(size, padded_size)
    return vec([TiledArrays.genindex(indexer, Tuple(I)) for I in CartesianIndices(size)])
end


text:
using SemidefiniteModels
using Compat
using Compat.Test
using Compat.LinearAlgebra

include("sdinterface.jl")
using CSDP
sdtest(CSDP.CSDPSolver(), duals=true)


text:
module DataStructures

### Constants

const DEFAULT_MAX_SIZE = 1024

### List related

import Base: append!, insert!, pop!, push!
import Base: eltype, isempty, length, show

export AbstractList, AbstractNode
export ArrayBasedList, LinkedList, LinkNode

export AbstractStack
export LinkedStack

export AbstractQueue
export LinkedQueue, PriorityQueue, KeyValuePair

export AbstractDictionary

export movetostart!, movetoposition!, movetoend!, prev!, next!, currentposition, getvalue
export topvalue

export enqueue!, dequeue!, frontvalue

export remove!, removeany!, find

# TODO put linknode out into own file and put it under abstract includes
include("./list/abstract_list.jl")
include("./list/abstract_queue.jl")
include("./list/abstract_stack.jl")
include("./list/abstract_dictionary.jl")

include("./list/array_based_list.jl")
#include("./list/array_based_queque.jl")
#include("./list/array_based_stack.jl")

include("./list/linked_list.jl")
include("./list/linked_queue.jl")
include("./list/linked_stack.jl")

### Binary trees related

import Base: getkey, parent, peek

export BinaryNode
export BinarySearchNode

export getvalue, setvalue!, setleftchild!, setrightchild!, isleaf
export getkey, setkey!

export AbstractTree
export BinaryHeap, BinarySearchTree, CompleteBinaryTree, MaxHeap, MinHeap

export parent, leftchild, rightchild, leftsibling, rightsibling
export buildheap!

include("./tree/abstract_binarynode.jl")
include("./tree/abstract_tree.jl")

include("./tree/binarysearch_node.jl")
include("./tree/binarysearch_tree.jl")
include("./tree/completebinary_tree.jl")
include("./tree/heap.jl")

include("./list/priorityqueue.jl") # belongs to list section

### Graph related

import Base: first, ==

export AbstractGraph
export AdjacencyList, AdjacencyMatrix, Vertex

export countedges, countvertices, first, next, setedge!, remove!, isedge
export getweight, getmark, setmark!

export breadthfirstsearch!, depthfirstsearch!, topologicalsort!

export dijkstra, jarnik

include("./graph/abstract_graph.jl")

include("./graph/adjacency_list.jl")
include("./graph/adjacency_matrix.jl")

include("./graph/bfs.jl")
include("./graph/dfs.jl")
include("./graph/topological_sort.jl")

include("./graph/utils.jl")
include("./graph/dijkstra.jl")
include("./graph/jarnik.jl")

### Internal sorting related

export insertionsort!, mergesort!, selectionsort!, shellsort!
export quicksort!, heapsort!

include("./sort/swap.jl")
include("./sort/heapsort.jl")
include("./sort/insertion_sort.jl")
include("./sort/mergesort.jl")
include("./sort/selection_sort.jl")
include("./sort/shellsort.jl")
include("./sort/quicksort.jl")

### Search related

export binarysearch

include("./search/binarysearch.jl")

end


text:
using Test, AdmittanceModels
using ProgressMeter: @showprogress

# set to false to create the figures in the paper, though it will take longer
run_fast = true
# to save data so figures in the paper can be created, set to true
write_data = false
# to make plots and print information set to true
show_results = false

if write_data
    using DataFrames, CSV, NPZ
    output_folder = joinpath(@__DIR__, "radiative_loss_results")
    try
        mkdir(output_folder)
    catch SystemError
    end
end

if show_results
    using PlotlyJS
end

#######################################################
# Shared parameters
#######################################################

begin
    ν = 1.2e8 # propagation_speed
    Z0 = 50.0 # characteristic_impedance
    δ = run_fast ? 200e-6 : 50e-6 # discretization length
    cg = 7e-15 # qubit - resonator coupling capacitance
    cj = 100e-15 # qubit capacitance
    lj = 10e-9 # qubit inductance
    x0 = 800e-6 # coupler location on resonator
end

#######################################################
# Transmission line model
#######################################################

begin
    L = .00499171 # resonator length
    F = 2e-3 # length of tline
    y0 = F/2 # coupler location on tline
    cc = 10e-15 # resonator - tline coupling capacitance
    tline_params = (L=L, x0=x0, F=F, y0=y0, cc=cc)
end

function create_tline_components(δ::Real)
    resonator = TransmissionLine(["qr_coupler", "er_coupler_0", "mode_2_max", "short"],
                                 ν, Z0, L, locations=[x0, 2*L/3], δ=δ)
    environment = TransmissionLine(["in", "er_coupler_1", "out"],
                                   ν, Z0, F, locations=[y0], δ=δ)
    er_coupler = SeriesComponent("er_coupler_0", "er_coupler_1", 0, 0, cc)
    qr_coupler = SeriesComponent("qubit", "qr_coupler", 0, 0, cg)
    return [resonator, environment, er_coupler, qr_coupler]
end
create_tline_pso(δ) = short_ports(
    connect(PSOModel.(create_tline_components(δ))), "short")

tline_cascade_comp    = create_tline_components(δ)
tline_cascade_pso     = create_tline_pso(δ)
tline_cascade_bbox(ω) = short_ports(
    connect(Blackbox.(Ref(ω), create_tline_components(δ))), "short")

#######################################################
# Purcell filter model
#######################################################

begin
    L = .00503299# resonator length
    F = .0102522 # length of pfilter
    y0 = 5e-3 # coupler location on pfilter
    cc = 2.6e-15 # resonator - pfilter coupling capacitance
    st = 880e-6 # pfilter stub length
    pfilter_params = (L=L, x0=x0, F=F, y0=y0, cc=cc, st=st)
end

function create_pfilter_components(δ::Real)
    resonator = TransmissionLine(["qr_coupler", "er_coupler_0", "mode_2_max", "short 1"],
                                 ν, Z0, L, locations=[x0, 2*L/3], δ=δ)
    environment = TransmissionLine(["short 2", "in", "er_coupler_1", "out", "short 3"],
                                   ν, Z0, F, locations=[st, y0, F-st], δ=δ)
    er_coupler = SeriesComponent("er_coupler_0", "er_coupler_1", 0, 0, cc)
    qr_coupler = SeriesComponent("qubit", "qr_coupler", 0, 0, cg)
    return [resonator, environment, er_coupler, qr_coupler]
end
create_pfilter_pso(δ) = short_ports(connect(
    PSOModel.(create_pfilter_components(δ))), ["short 1", "short 2", "short 3"])

pfilter_cascade_comp = create_pfilter_components(δ)
pfilter_cascade_pso  = create_pfilter_pso(δ)
pfilter_cascade_bbox(ω) = short_ports(connect(
    Blackbox.(Ref(ω), pfilter_cascade_comp)), ["short 1", "short 2", "short 3"])
#######################################################
# Scattering parameters with qubit
#######################################################

begin
    factor = run_fast ? 2 : 10
    ω = [range(1.5, stop=5.89, length=10 * factor);
         range(5.89, stop=5.92, length=60 * factor); # first mode for both
         range(5.92, stop=17.69, length=10 * factor);
         range(17.69, stop=17.72, length=100 * factor); # 2nd mode pfilter
         range(17.72, stop=17.76, length=10 * factor);
         range(17.76, stop=17.85, length=100 * factor); # 2nd mode tline
         range(17.85, stop=25, length=10 * factor)] * 2π * 1e9
    qubit = ParallelComponent("qubit", 1/lj, 0, cj)

    tline_pso_bbox = let
        model = connect(tline_cascade_pso, PSOModel(qubit))
        model = open_ports_except(model, ["in", "out"])
        Blackbox(ω, model)
    end
    tline_pso_bbox_s = [x[1,2] for x in scattering_matrices(tline_pso_bbox, [Z0, Z0])]
    tline_bbox = let
        model = connect(tline_cascade_bbox(ω), Blackbox(ω, qubit))
        model = open_ports_except(model, ["in", "out"])
    end
    tline_bbox_s = [x[1,2] for x in scattering_matrices(tline_bbox, [Z0, Z0])]
    pfilter_pso_bbox = let
        model = connect(pfilter_cascade_pso, PSOModel(qubit))
        Blackbox(ω, open_ports_except(model, ["in", "out"]))
    end
    pfilter_pso_bbox_s = [x[1,2] for x in scattering_matrices(pfilter_pso_bbox, [Z0, Z0])]
    pfilter_bbox = let
        model = connect(pfilter_cascade_bbox(ω), Blackbox(ω, qubit))
        open_ports_except(model, ["in", "out"])
    end
    pfilter_bbox_s = [x[1,2] for x in scattering_matrices(pfilter_bbox, [Z0, Z0])]
end

if write_data
    CSV.write(joinpath(output_folder, "S12.csv"), DataFrame(
        freqs = ω/(2π),
        tline_S = tline_pso_bbox_s,
        pfilter_S = pfilter_pso_bbox_s))
end

if show_results
    plot([scatter(x=ω/(2π*1e9), y=abs.(tline_pso_bbox_s), name="tline pso"),
          scatter(x=ω/(2π*1e9), y=abs.(tline_bbox_s), name="tline bbox"),
          scatter(x=ω/(2π*1e9), y=abs.(pfilter_pso_bbox_s), name="pfilter pso", line_dash="dash"),
          scatter(x=ω/(2π*1e9), y=abs.(pfilter_bbox_s), name="pfilter bbox", line_dash="dash")],
          Layout(xaxis_title="frequency", yaxis_title="|S21|"))
end

#######################################################
# Find effective capacitance
#######################################################

function terminated(model::PSOModel)
    resistors = ParallelComponent.(["in", "out"], 0, 1/Z0, 0)
    return connect([model; PSOModel.(resistors)])
end
function admittances(model::Blackbox)
    resistors = ParallelComponent.(["in", "out"], 0, 1/Z0, 0)
    model = connect([model; Blackbox.(Ref(model.ω), resistors)])
    bb = canonical_gauge(open_ports_except(model, "qubit"))
    return [x[1,1] for x in admittance_matrices(bb)]
end

begin
    fit_ω = collect(range(0.1e9, stop=2.5e9, length=run_fast ? 100 : 1000)) * 2π
    plot_ω = collect(range(0.01e9, stop=10e9, length=run_fast ? 100 : 1000)) * 2π
    linreg(x, y) = hcat(fill!(similar(x), 1), x) \ y
    tline_pso, Y = terminated(tline_cascade_pso), admittances(tline_cascade_bbox(fit_ω))
    _, tline_slope = linreg(fit_ω, imag.(Y))
    tline_Y = admittances(tline_cascade_bbox(plot_ω))
    Y = admittances(pfilter_cascade_bbox(fit_ω))
    _, pfilter_slope = linreg(fit_ω, imag.(Y))
    pfilter_pso, pfilter_Y = terminated(tline_cascade_pso), admittances(tline_cascade_bbox(plot_ω))
    tline_eff_capacitance = tline_slope + cj
    pfilter_eff_capacitance = pfilter_slope + cj
end

if show_results
    println("tline_slope: $tline_slope, pfilter_slope: $pfilter_slope, cg: $cg")
    println("tline_eff_capacitance: $tline_eff_capacitance, pfilter_eff_capacitance: $pfilter_eff_capacitance")
    density(m) = count(!iszero, m)/(size(m,1) * size(m,2))
    println("tline matrix densities: $(map(density, get_Y(tline_pso)))")
    println("pfilter matrix densities: $(map(density, get_Y(tline_pso)))")
end

if write_data
    CSV.write(joinpath(output_folder, "Y.csv"), DataFrame(
        freqs = plot_ω/(2π),
        tline_Y = tline_Y,
        pfilter_Y = pfilter_Y))
end

if show_results
    plot([scatter(x=plot_ω/(2π), y=imag.(tline_Y), name="imag(tline_Y)"),
          scatter(x=plot_ω/(2π), y=imag.(pfilter_Y), name="imag(pfilter_Y)", line_dash="dash"),
          scatter(x=plot_ω/(2π), y=cg * plot_ω, name="cg * ω")],
          Layout(xaxis_title="frequency", yaxis_title="imag(Y)", yaxis_range=[-.001, .002]))
end

#######################################################
# Eigenmodes
#######################################################

function get_bare_qubit_freqs(num_points)
    f0 = ν/(4 * tline_params.L)
    f1 = ν/(2 * (tline_params.L - tline_params.x0))
    f2 = 3 * f0
    _eps = f0/10
    freqs_0 = range(1.5e9, stop=f0-_eps, length=num_points)
    freqs_1 = range(f0-_eps, stop=f0+_eps, length=num_points)
    freqs_2 = range(f0+_eps, stop=f1-_eps, length=num_points)
    freqs_3 = range(f1-_eps, stop=f1+_eps, length=num_points)
    freqs_4 = range(f1+_eps, stop=f2-_eps, length=num_points)
    freqs_5 = range(f2-_eps, stop=f2+_eps, length=num_points)
    freqs_6 = range(f2+_eps, stop=25e9, length=num_points)
    return [freqs_0; freqs_1; freqs_2; freqs_3; freqs_4; freqs_5; freqs_6]
end

function eigenmodes(lj::Real, pso::PSOModel)
    qubit = ParallelComponent("qubit", 1/lj, 0, cj)
    pso_full = connect(pso, PSOModel(qubit))
    eigenvalues, eigenvectors = lossy_modes_dense(pso_full, min_freq=1e9, max_freq=25e9)
    port_inds = ports_to_indices(pso_full, "qubit", "qr_coupler", "mode_2_max")
    mode_inds = AdmittanceModels.match_vectors(get_P(pso_full)[:, port_inds], eigenvectors)
    return eigenvalues[mode_inds]
end

function eigenmodes(casc::PSOModel, eff_capacitance::Real, ω::Vector{<:Real})
    pso, Y = terminated(casc), admittances(Blackbox(ω, casc))
    T1_Y = eff_capacitance ./ real.(Y)
    ljs = (1 ./ (ω .^2 * eff_capacitance))
    if show_results
        println("dimension: $(size(pso.K,1)), num_ljs: $(length(ljs))")
    end
    λ_lists = @showprogress [eigenmodes(lj, pso) for lj in ljs]
    λs = cat(λ_lists..., dims=2)
    return λs, T1_Y
end

ω = 2π * get_bare_qubit_freqs(run_fast ? 10 : 120)
tline_λs, tline_T1_Y = eigenmodes(tline_cascade_pso, tline_eff_capacitance, ω)

if show_results
    plot([scatter(x=ω/(2π*1e9), y=imag.(tline_λs[i,:])/(2π*1e9), name="$i",
        mode="lines+markers") for i in 1:size(tline_λs,1)],
        Layout(xaxis_title="bare qubit freq", yaxis_title="mode freq", title="tline freqs"))
end

if show_results
    traces = [scatter(x=ω/(2π*1e9), y=-.5 ./ real.(tline_λs[i,:]), name="$i",
               mode="lines+markers") for i in 1:size(tline_λs,1)]
    push!(traces, scatter(x=ω/(2π*1e9), y=tline_T1_Y, name="C/Re[Y]", mode="lines+markers"))
    plot(traces, Layout(xaxis_title="bare qubit freq", yaxis_title="T1", yaxis_type="log", title="tline T1s"))
end

if write_data
    df = DataFrame(transpose(tline_λs))
    df[:freqs] = ω/(2π)
    df[:T1_Y] = tline_T1_Y
    CSV.write(joinpath(output_folder, "tline_eigenmodes.csv"), df)
end

pfilter_λs, pfilter_T1_Y = eigenmodes(pfilter_cascade_pso, pfilter_eff_capacitance, ω)

if show_results
    plot([scatter(x=ω/(2π*1e9), y=imag.(pfilter_λs[i,:])/(2π*1e9), name="$i",
        mode="lines+markers") for i in 1:size(pfilter_λs,1)],
        Layout(xaxis_title="bare qubit freq", yaxis_title="mode freq", title="pfilter freqs"))
end

if show_results
    traces = [scatter(x=ω/(2π*1e9), y=-.5 ./ real.(pfilter_λs[i,:]), name="$i",
        mode="lines+markers") for i in 1:size(pfilter_λs,1)]
        push!(traces, scatter(x=ω/(2π*1e9), y=pfilter_T1_Y, name="C/Re[Y]", mode="lines+markers"))
        plot(traces, Layout(xaxis_title="bare qubit freq", yaxis_title="T1", yaxis_type="log", title="pfilter T1s"))
end

if write_data
    df = DataFrame(transpose(pfilter_λs))
    df[:freqs] = ω/(2π)
    df[:T1_Y] = pfilter_T1_Y
    CSV.write(joinpath(output_folder, "pfilter_eigenmodes.csv"), df)
end

#######################################################
# Convergence
#######################################################

freqs_func(vects, mode_ind, δ_ind) = imag.(vects[δ_ind][mode_ind,:])/(2π*1e9)
T1_func(vects, mode_ind, δ_ind) = -.5 ./ real.(vects[δ_ind][mode_ind,:])
function diff(vects, func, mode_ind, δ_ind_bad, δ_ind_good)
    f_bad = func(vects, mode_ind, δ_ind_bad)
    f_good = func(vects, mode_ind, δ_ind_good)
    return abs.((f_bad - f_good) ./ f_good) * 100
end

if !run_fast
    δs = [60, 50, 40] * 1e-6
    tline_λ_vects = [eigenmodes(create_tline_pso(δ), tline_eff_capacitance, ω)[1] for δ in δs]

    if write_data
        npzwrite(joinpath(output_folder, "tline_eigenvalues.npy"), cat(tline_λ_vects..., dims=3))
    end

    if show_results
        plot([scatter(x=ω/(2π*1e9), y=diff(tline_λ_vects, freqs_func, i, 1, 3), name="$i", mode="markers") for i in 1:3],
            Layout(xaxis_title="bare qubit freq", yaxis_title="% error", yaxis_type="log", title="tline freqs"))
    end

    if show_results
        plot([scatter(x=ω/(2π*1e9), y=diff(tline_λ_vects, T1_func, i, 2, 3), name="$i", mode="markers") for i in 1:3],
            Layout(xaxis_title="bare qubit freq", yaxis_title="% error", yaxis_type="log", title="tline freqs"))
    end

    pfilter_λ_vects = [eigenmodes(create_pfilter_pso(δ), pfilter_eff_capacitance, ω)[1] for δ in δs]

    if write_data
        npzwrite(joinpath(output_folder, "pfilter_eigenvalues.npy"), cat(pfilter_λ_vects..., dims=3))
    end

    if show_results
        plot([scatter(x=ω/(2π*1e9), y=diff(pfilter_λ_vects, freqs_func, i, 2, 3), name="$i", mode="markers") for i in 1:3],
            Layout(xaxis_title="bare qubit freq", yaxis_title="% error", yaxis_type="log", title="tline freqs"))
    end

    if show_results
        plot([scatter(x=ω/(2π*1e9), y=diff(pfilter_λ_vects, T1_func, i, 2, 3), name="$i", mode="markers") for i in 1:3],
            Layout(xaxis_title="bare qubit freq", yaxis_title="% error", yaxis_type="log", title="tline freqs"))
    end
end

#######################################################
# Declare victory
#######################################################

@testset "radiative loss" begin
    @test true
end


text:
"""
`rpc(r,g)`: Regional purchase coefficient

```math
\\rho_{r,g}^{cfs}
=
\\begin{cases}
\\dfrac{\\bar{d}_{r,g}}
       {\\bar{d}_{r,g} - \\bar{mn}_{r,g}}  & r\\neq uti, \\bar{d}_{r,g} \\neq \\bar{mn}_{r,g}
\\\\
0.0                                        & r\\neq uti, \\bar{d}_{r,g} = \\bar{mn}_{r,g}
\\\\
0.9                                        & r = uti
\\end{cases}
```
"""
function share_rpc!(d::Dict, set::Dict)
    println("  Calculating rpc(r,g), regional purchase coefficient")

    _set_ng!(d, set)
    _share_mrt0!(d)
    _share_d0!(d, set)
    _share_xn0!(d, set)
    _share_mn0!(d, set)

    df = dropnan(d[:d0] / (d[:d0] + d[:mn0]))

    df = edit_with(df, Drop(:g, "uti", "=="))
    df_uti = fill_with((r = set[:r], g = "uti"), 0.9)

    d[:rpc] = sort([df; df_uti])
end

"""
`ng`: Sectors not included in the CFS.
"""
function _set_ng!(d::Dict, set::Dict)
    set[:ng] = setdiff(set[:g], unique(d[:cfs][:,:g]))
end

"""
`d0(r,g)`: Local supply-demand (CFS), trade that remains within the same region.

```math
\\bar{d}_{r,ng\\ni g} = \\left\\{{cfs}\\left(orig,dest,g\\right)
\\;\\vert\\; orig=dest, \\, g \\right\\}
```
Calling [`SLiDE._avg_ng`](@ref) returns ``\\bar{d}_{r,ng\\in g}``.
"""
function _share_d0!(d::Dict, set::Dict)
    println("\td0(r,g), local supply-demand")
    df = copy(d[:cfs])

    df = df[df[:,:orig] .== df[:,:dest],:]
    df = edit_with(df, [Rename(:orig, :r), Deselect([:dest], "==")])

    d[:d0] = _avg_ng(df, set)
    return d[:d0]
end

"""
`mrt0(orig,dest,g)`: Interstate trade (CFS)
```math
\\bar{mrt}_{orig,dest,ng\\ni g} = \\left\\{{cfs}\\left(orig,dest,g\\right)
\\;\\vert\\; orig\\neq dest, \\, g \\right\\}
```
"""
function _share_mrt0!(d::Dict)
    df = copy(d[:cfs])

    d[:mrt0] = df[df[:,:orig] .!= df[:,:dest],:]
    return d[:mrt0]
end

"""
`mn0(r,g)`: National demand (CFS)
```math
\\bar{mn}_{r,ng\\ni g} = \\sum_{orig} \\bar{mrt}_{orig,dest,ng\\ni g}
```
Calling [`SLiDE._avg_ng`](@ref) returns ``\\bar{mn}_{r,ng\\in g}``.
"""
function _share_mn0!(d::Dict, set::Dict)
    println("\tmn0(r,g), national demand")
    df = copy(d[:mrt0])

    df = edit_with(combine_over(df, :orig), Rename(:dest, :r))
d[:mn0] = _avg_ng(df, set)
    return d[:mn0]
end

"""
`xn0(r,g)`: National exports (CFS)
```math
\\bar{xn}_{r,ng\\ni g} = \\sum_{dest} \\bar{mrt}_{orig,dest,ng\\ni g}
```
Calling [`SLiDE._avg_ng`](@ref) returns ``\\bar{xn}_{r,ng\\in g}``.
"""
function _share_xn0!(d::Dict, set::Dict)
    println("\txn0(r,g), national exports")
    df = copy(d[:mrt0])

    df = edit_with(combine_over(df, :dest), Rename(:orig, :r))
    d[:xn0] = _avg_ng(df, set)
    return d[:xn0]
end

"""
```math
\\bar{x}_{r,ng\\in g} = \\dfrac{\\sum_g \\bar{x}_{r,g}}
    {\\text{length}(ng)}
```
"""
function _avg_ng(df::DataFrame, set::Dict)
    df_ng = crossjoin(
        DataFrame(g=set[:ng]),
        combine_over(copy(df), :g; fun=Statistics.sum),
    )

    df_ng[!,:value] ./= length(setdiff(set[:g], set[:ng]))

    return vcat(df, df_ng)
end

max_stars_repo_name: NREL/SLiDE

max_issues_repo_name: NREL/SLiDE

max_forks_repo_name: NREL/SLiDE


text:
using Base.Test

include("leap.jl")

@testset "Year not divisible by 4: common year" begin
  @test !is_leap_year(2015)
end

@testset "Year divisible by 4, not divisible by 100: leap year" begin
  @test is_leap_year(2016)
end

@testset "Year divisible by 100, not divisible by 400: common year" begin
  @test !is_leap_year(2100)
end

@testset "Year divisible by 400: leap year" begin
  @test is_leap_year(2000)
end

max_stars_repo_name: dceresoli/exercism

max_issues_repo_name: dceresoli/exercism

max_forks_repo_name: dceresoli/exercism


text:
module EigenvalueSolver

# A Julia package for solving systems of polynomial equations using eigenvalue methods.
# Matías R. Bender and Simon Telen - 14-05-2021

using DynamicPolynomials
using LinearAlgebra
using GenericSchur
using GenericSVD
using Statistics
using SmithNormalForm  # to add: pkg> add https://github.com/wildart/SmithNormalForm.jl.git

import TaylorSeries

include("PolytopeUtils.jl")

# Main function.
# Solves the system f = 0 in the variables x using the Sylvester map defined by A₀, E, D.
# INPUT:
# f     vector of polynomials of length s
# x     array of variables
# A₀    array representing the exponents of the monomials used to construct generalized multiplication matrices
# E     vector of arrays of length s+1, represents the monomial multipliers of f₀, f₁, …, fs (rows are exponent vectors)
# D     array representing a monomial basis of the codomain of the Sylvester map (rows are exponent vectors)
# OUTPUT:
# a set of candidate solutions to f(x) = 0.
# ------------------ It is assumed that (f,A₀,E,D) is an admissible tuple
# ------------------ and that all arrays in A₀, E, D have the same number of columns (n)
# ------------------ If the cokernel N of the Sylvester map was already computed, it is possible to pass N as an input
function solve_EV(
    f,
    x,
    A₀,
    E,
    D;
    rankTol = 1e3 * eps(), #tolerance for rank determination and nullspace computation
    complex = true,
    check_criterion = true,
    EVtol = 1e-4, #tolerance for extracting eigenvalues
    N = nothing,
    compress = true, #compress Sylv before computing cokernel
    schurfac = false, #compute eigenvalues via Schur factorization
    filter = false, #filter the solutions at the end
    filtertol = 1e-6,
    verbose::Bool=false, #print info
)
    x = variables(f)
    n = length(x) # number of variables
    s = length(f) # number of equations, s ≥ n
    #-------------------------------------------------------------------------
    Σ = exptomon(D, x)
    if N == nothing
        verbose && println("constructing Sylvester matrix...")
        σ = [exptomon(e, x) for e ∈ E[2:end]]
        @time Sylv = getRes(f, Σ, σ, x; complex = complex)
        verbose && println("       Sylv is a matrix of size $(size(Sylv))")
        #-------------------------------------------------------------------------
        verbose && println("computing cokernel...")
        if compress && size(Sylv, 2) > size(Sylv, 1)
            verbose && println("compressing Sylv...")
            @time Sylv = Sylv * randn(size(Sylv, 2), size(Sylv, 1))
        end
        @time svdobj = svd(transpose(Sylv), full = true)
        firstzerosingval = findfirst(sv -> sv < rankTol * svdobj.S[1], svdobj.S)
        if !isnothing(firstzerosingval)
            rk = firstzerosingval - 1
            verbose && println("       rank = $(rk)")
            verbose && println("       relative size of last nz singular value = $(svdobj.S[rk]/svdobj.S[1])")
            verbose && println("       gap = $(svdobj.S[rk]/svdobj.S[rk+1])")
        else
            rk = length(svdobj.S)
            verbose && println("       relative size of last nz singular value = $(svdobj.S[rk]/svdobj.S[1])")
        end
        N = transpose(svdobj.V[:, rk+1:end])
    end
    verbose && println("N has size $(size(N))")
    #-------------------------------------------------------------------------
    Σ₀ = exptomon(E[1], x)
    Σα0 = exptomon(A₀, x)
    if check_criterion
        verbose && println("Checking criterion...")
        chck = checkCriterion(Σα0, Σ₀, Σ, N, x; verbose = verbose)
    end
    #-------------------------------------------------------------------------
    if !check_criterion || chck
        verbose && println("finding a basis using QR-P...")
        @time M, κ = get_multmtces_QR_toric(N, Σ, Σ₀, Σα0, x; complex = complex)
        #-------------------------------------------------------------------------
        verbose && println("extracting relevant eigenvalues...")
        if schurfac
            @time monsol =
                simulDiag_schur_toric(M; clustering = true, tol = 1e-4, complex = complex)
        else
            @time monsol = extractSolutions_fast(M, Σα0; tol = EVtol)
        end
        #-------------------------------------------------------------------------
        verbose && println("inverting monomial map given by A₀...")
        A₀2 = convert(Array{Int64,2}, A₀')
        zeroind = findfirst(ℓ -> norm(A₀2[:, ℓ]) == 0, 1:size(A₀2, 2))
        @time logsol, sol = solvebinom_SNF_array(A₀2, [ms / ms[zeroind] for ms ∈ monsol])
        verbose && println("found $(length(sol)) solutions")
        if filter
            verbose && println("filtering...")
            res = get_residual(f, sol, x)
            sol = sol[findall(r -> r < filtertol, res)]
        end
        return sol
    else
        println("Criterion was not satisfied.")
    end
end

# Returns a vector of monomials in x corresponding to the exponent tuples in the rows of M.
function exptomon(M, x)
    mons = fill(prod(x), size(M, 1))
    mtx = x' .^ M
    mons = map(k -> prod(mtx[k, :]), 1:size(mtx, 1)) # this gives a vector of monomials
end

# Converts an array of monomials into an array of exponents.
function montoexp(Σ)
    exparray = exponents.(Σ)
    M = fill(0, length(Σ), length(exparray[1]))
    for i = 1:length(Σ)
        M[i, :] = exparray[i]
    end
    return M
end

# Computes the Sylvester map in the monomial bases σ[i] for R_{A_i} and Σ for R_D.
# x should contain the variables of the system.
function getRes(f, Σ, σ, x; complex = false)
    if complex
        res = fill(0.0 + 0.0 * im, length(Σ), sum(map(i -> length(i), σ)))
    else
        res = fill(0.0, length(Σ), sum(map(i -> length(i), σ)))
    end
    n = length(x)
    s = length(f)
    mapping = Dict(Σ .=> 1:length(Σ))
    for i = 1:s
        J = sum(map(k -> length(k), σ[1:i-1]))
        for j = 1:length(σ[i])
            pol = σ[i][j] * f[i]
            res[map(k -> mapping[k], monomials(pol)), J+j] = coefficients(pol)
        end
    end
    return res
end

# Check if N_{f₀} is of rank rank(N) for generic elements f₀ in span(Σα0).
# Here N is the previously computed cokernel of a Sylvester map with support D.
# The monomials in Σ₀ should correspond to E₀ and those in Σ to D.
function checkCriterion(Σα0, Σ₀, Σ, N, x; rankTol = 1e3 * eps(), verbose::Bool=false)
    f₀ = randn(ComplexF64, length(Σα0))' * Σα0
    M0 = getRes([f₀], Σ, [Σ₀], x; complex = true)
    r = rank(N * M0; rtol = rankTol)
    if r == size(N, 1)
        verbose && println("********* criterion satisfied ********* γ = ", r)
        return true
    else
        verbose && println("!!!!!!!!! criterion violated !!!!!!!!!")
        return false
    end
end

# This computes the multiplication matrices via a basis selection which uses QR with optimal column pivoting
function get_multmtces_QR_toric(N, Σ_α_α₀, Σ_α, Σ_α₀, t; complex = false)
    δ = size(N, 1)
    mapping = Dict(Σ_α_α₀ .=> 1:length(Σ_α_α₀))
    n = length(t)
    n_α₀ = length(Σ_α₀)
    if complex
        Nᵢ = fill(fill(0.0 + 0.0 * im, δ, length(Σ_α)), n_α₀)
        M = fill(fill(0.0 + 0.0 * im, δ, δ), n_α₀) # This will contain the multiplication operators corresponding to the variables
    else
        Nᵢ = fill(fill(0.0, δ, length(Σ_α)), n_α₀)
        M = fill(fill(0.0, δ, δ), n_α₀) # This will contain the multiplication operators corresponding to the variables
    end

    for i = 1:n_α₀
        indsᵢ = map(k -> mapping[k*Σ_α₀[i]], Σ_α)
        Nᵢ[i] = N[:, indsᵢ]
    end
    if complex
        Nh₀ = sum(randn(ComplexF64, n_α₀) .* Nᵢ)
    else
        Nh₀ = sum(randn(n_α₀) .* Nᵢ) # Dehomogenization with respect to a generic linear form.
    end
    QRobj = qr(Nh₀, Val(true))
    pivots = QRobj.p
    Nb = triu(QRobj.R[1:δ, 1:δ])
    κ = cond(Nb)
    if κ > 1e12
        println("warning: this might not be zero-dimensional in X")
    end
    for i = 1:n_α₀
        M[i] = Nb \ (QRobj.Q' * Nᵢ[i][:, pivots[1:δ]])
    end
    return M, κ
end

# This computes an upper triangularized, column pivoted version of N via a pivoted QR factorization of
# the map N: R_D → C^γ whose columns are indexed by the monomials in Σ_α_α₀ ~ D.
# the monomials in the subspace W (largest subspace s.t. W⁺ ⊂ V) are in Σ_α ~ E₀.
function get_QRbasis_general(N, Σ_α_α₀, Σ_α)
    mapping = Dict(Σ_α_α₀ .=> 1:length(Σ_α_α₀))
    δ = size(N, 1)
    indsW = map(k -> mapping[k], Σ_α)
    indsNotW = findall(i -> i ∉ indsW, 1:length(Σ_α_α₀))
    QRobj = qr(N[:, indsW], Val(true))
    pivots = QRobj.p
    N = triu(QRobj.Q' * N[:, vcat(indsW[pivots], indsNotW)])
    Σ = Σ_α_α₀[vcat(indsW[pivots], indsNotW)]
    return N, Σ
end

# simultaneous diagonalization of the matrices in M via Schur factorization of a random linear combination.
# returns an array with the solutions (`third mode` of the CPD)
function simulDiag_schur_toric(M; clustering = false, tol = 1e-4, complex = false)
    solMatrix =
        simulDiag_schur_mtx(M; clustering = clustering, tol = tol, complex = complex)
    δ = size(solMatrix, 1)
    monsol = map(i -> solMatrix[i, :], 1:δ)
    return monsol
end

# This is a function for clustering the eigenvalues.
function cluster(v; tol = 1e-4)
    clusters = fill(zero(typeof(v[1])), 0)
    inds = fill([0], 0)
    currentcluster = [v[1]]
    currentinds = [1]
    remaininginds = 2:length(v)
    lastind = 1
    while length(remaininginds) > 0
        nextind = remaininginds[argmin(abs.(v[remaininginds] .- v[lastind]))]
        if abs(v[nextind] - v[lastind]) < tol * (abs(v[lastind]) + 1)
            push!(currentcluster, v[nextind])
            push!(currentinds, nextind)
        else
            push!(clusters, mean(currentcluster))
            push!(inds, currentinds)
            currentcluster = [v[nextind]]
            currentinds = [nextind]
        end
        remaininginds = setdiff(remaininginds, nextind)
        lastind = nextind
    end
    push!(clusters, mean(currentcluster))
    push!(inds, currentinds)
    return clusters, inds
end

# Match the points in sol to the points in refsol by permuting them as defined by 'usedinds'
# and compute the relative distance between the matched points.
function getForwardError(sol, refsol)
    errs = fill(0.0, 0)
    usedinds = fill(1, 0)
    δ = length(refsol)
    for s ∈ sol
        remaininginds = setdiff(1:δ, usedinds)
        if length(remaininginds) > 1
            nrms = norm.([(ref - s)/norm(ref) for ref ∈ refsol[remaininginds]])
            amin = argmin(nrms)
            matchind = remaininginds[amin]
            push!(usedinds, matchind)
            push!(errs, nrms[amin])
        else
            push!(usedinds, remaininginds[1])
            push!(errs, norm(refsol[remaininginds[1]] - s)/norm(refsol[remaininginds[1]]))
        end
    end
    return errs, usedinds
end

# Extracts the useful eigenvalues from an array M of matrices M_{x^α} for all
# monomials in x^α ∈ Σα0.
function extractSolutions_fast(M, Σα0; tol = 1e-5, clustertol = 1e-7)
    Mg = sum([randn(ComplexF64) * MM for MM ∈ M])
    Mgprime = sum([randn(ComplexF64) * MM for MM ∈ M])
    eigenobj = eigen(transpose(Mg))
    eigenvals = eigenobj.values
    eigenvecs = eigenobj.vectors
    clusters, inds = cluster(eigenvals; tol = clustertol)
    monsol = []
    γ = size(Mg, 1)
    k = 1
    for i = 1:length(clusters)
        λ = clusters[i]
        Vλ = eigenvecs[:, inds[i]]
        if size(Vλ, 2) > 0
            if size(Vλ, 2) == 1
                aux = [transpose(Vλ) * Mgprime; transpose(Vλ)]
                svdobj = svd(aux)
                ratio = svdobj.S[2] / svdobj.S[1]
                if ratio < tol
                    eval = [mean(aux[1, :] ./ aux[2, :])]
                    espace = [[1]]
                else
                    eval = []
                    espace = []
                end
            else
                Vλ, = svd(Vλ)
                eval, espace = solveREP(transpose(Vλ) * Mgprime, transpose(Vλ); rtol = tol)
            end
        end
        if length(eval) >= 1
            Vλnew = hcat(espace...)
            Vλ = Vλ * Vλnew'
            m = size(Vλ, 2)
            if m == 1 && !isnan(eval[1])
                monsol = push!(monsol, [(transpose(Vλ)*MM*Vλ)[1] for MM ∈ M])
            elseif m > 1 && !isnan(eval[1])

                O = randn(ComplexF64, γ, m)
                aux = fill(0.0 + 0.0 * im, length(M))
                for i = 1:length(M)
                    aux[i] = tr((transpose(Vλ) * M[i] * O / (transpose(Vλ) * O))) / m
                end
                monsol = push!(monsol, aux)
            end
        end
        k += 1
    end
    return monsol
end

# notation of `Numerical Root Finding via Cox Rings`
# solves the binomial systems {xᵃ = Λ[i]ⱼ} for i = 1...length(Λ) and a is the j-th column of A.
function solvebinom_SNF_array(A, Λ)
    F = smith(A)
    k = size(A, 1)
    δ = length(Λ)
    n_α₀ = size(A, 2)
    invfactors = F.SNF
    r = findfirst(k -> k == 0, invfactors)  #rank
    if isnothing(r)
        r = minimum(size(A))
    else
        r = r - 1
    end
    U = inv(F.S)
    V = inv(F.T)
    Λmat = fill(0.0 + 0.0 * im, length(Λ), n_α₀)
    for i = 1:δ
        Λmat[i, :] = Λ[i]
    end
    w = log.(Complex.(Λmat)) * V[:, 1:r]
    w = w * diagm(1 ./ invfactors[1:r])
    logY = hcat(w, zeros(δ, k - r))
    logX = logY * U
    loghomsol = map(k -> logX[k, :], 1:δ)
    homsol = map(k -> exp.(loghomsol[k]), 1:δ)
    return loghomsol, homsol
end

# Compute the cokernel from a larger Sylvester map, supported on (oldsupport ∪ newsupport), from
# the cokernel 'oldcokernel' of a Sylvester map supported on oldsupport.
# The monomial multipliers of the equations in f to be added are given by an array `newshifts' of monomial vectors.
function updateCokernel(
    f,
    x,
    oldcokernel,
    oldsupport,
    newsupport,
    newshifts;
    complex = false,
    rankTol = 1e3 * eps(),
)
    newsylvprime = getRes(f, newsupport, newshifts, x; complex = complex)
    oldδ = size(oldcokernel, 1)
    growth = length(newsupport) - length(oldsupport)
    mapping = Dict(newsupport .=> 1:length(newsupport))
    oldinds = [mapping[ind] for ind ∈ oldsupport]
    newinds = setdiff(1:length(newsupport), oldinds)
    id = Matrix(I, growth, growth)
    if complex
        prodmtx = fill(0.0 + 0.0im, oldδ + growth, length(newsupport))
    else
        prodmtx = fill(0.0, oldδ + growth, length(newsupport))
    end
    prodmtx[1:oldδ, oldinds] = oldcokernel
    prodmtx[oldδ+1:end, newinds] = id
    L = transpose(nullspace(transpose(prodmtx * newsylvprime), rtol = rankTol))
    return L * prodmtx
end

# Compute an admissible tuple for a dense system of polynomials f by incrementing the degree
# until our criterion is satisfied.
function findAdmissibleTuple_dense(f, x, DMAX; complex = false, rankTol = 1e3 * eps(), verbose::Bool=false, trySemiregular::Bool=false)
    s = length(f)
    ds = fill(0, s)
    for i = 1:s
        ds[i] = maximum(degree.(terms(f[i])))
    end
    if trySemiregular
        varTaylor = TaylorSeries.Taylor1(Integer,sum(ds))
        HScandidate = prod(vv -> (1-varTaylor^vv),ds)/(1-varTaylor)^length(x)
        dmax = max(findfirst(vv -> vv <= 0,HScandidate.coeffs)-1,maximum(ds))
    else
        dmax = maximum(ds)
    end
    verbose && println("degree = ", dmax)
    oldsupport = getMonomials(x, dmax)
    oldshifts = [getMonomials(x, dmax - ds[i]) for i = 1:s]
    Σα0 = getMonomials(x, 1)
    Σ₀ = getMonomials(x, dmax - 1)
    oldSylv = getRes(f, oldsupport, oldshifts, x; complex = complex)
    oldcokernel = transpose(nullspace(transpose(oldSylv), rtol = rankTol))
    chck = checkCriterion(Σα0, Σ₀, oldsupport, oldcokernel, x; verbose = verbose)
    while dmax < DMAX && !chck
        dmax = dmax + 1
        verbose && println("degree = ", dmax)
        newsupport = getMonomials(x, dmax)
        newshifts = [getMonomials(x, dmax - ds[i]; homogeneous = true) for i = 1:s]
        newcokernel = updateCokernel(
            f,
            x,
            oldcokernel,
            oldsupport,
            newsupport,
            newshifts;
            complex = complex,
            rankTol = rankTol,
        )
        testshifts = [getMonomials(x, dmax - ds[i]; homogeneous = false) for i = 1:s]
        oldsupport = newsupport
        oldcokernel = newcokernel
        Σ₀ = getMonomials(x, dmax - 1)
        chck = checkCriterion(Σα0, Σ₀, oldsupport, oldcokernel, x; verbose = verbose)
    end
    topdeg = dmax
    A₀ = montoexp(getMonomials(x, 1))
    extdeg = topdeg .- ds
    extdeg = [topdeg - 1; extdeg]
    E = [montoexp(getMonomials(x, extdeg[i])) for i = 1:length(f)+1]
    D = montoexp(getMonomials(x, topdeg))
    return oldcokernel, oldsupport, chck, dmax, (A₀, E, D)
end

# Compute an admissible tuple for a sparse unmixed system of polynomials f by incrementing the degree
# until our criterion is satisfied.
function findAdmissibleTuple_unmixed(
    f,
    x,
    A,
    α,
    DMAX;
    complex = false,
    rankTol = 1e3 * eps(),
    verbose::Bool=false
)
    s = length(f)
    dmax = maximum(α)
    A₀pol = sum([prod(x .^ (A[i, :])) for i = 1:size(A, 1)])
    A₀polytope = newtonPolytope(A₀pol, x)
    vtcs = vertexRepresentation(A₀polytope)
    A₀pol = sum([prod(x .^ (vtcs[i, :])) for i = 1:size(vtcs, 1)])
    A₀polytope = newtonPolytope(A₀pol, x)
    Aold = newtonPolytope(A₀pol^dmax, x)
    oldlatticepoints = getLatticePoints(Aold)
    oldsupport = exptomon(oldlatticepoints, x)
    oldpolytopes = [newtonPolytope(A₀pol^(dmax - α[i]), x) for i = 1:s]
    oldshiftlatticepoints = [getLatticePoints(P) for P ∈ oldpolytopes]
    oldshifts = [exptomon(oslp, x) for oslp ∈ oldshiftlatticepoints]
    Σα0 = exptomon(getLatticePoints(A₀polytope), x)
    Σ₀ = exptomon(getLatticePoints(newtonPolytope(A₀pol^(dmax - 1), x)), x)
    oldSylv = getRes(f, oldsupport, oldshifts, x; complex = complex)
    oldcokernel = transpose(nullspace(transpose(oldSylv), rtol = rankTol))
    chck = checkCriterion(Σα0, Σ₀, oldsupport, oldcokernel, x; verbose = verbose)
    while dmax < DMAX && !chck
        dmax = dmax + 1
        verbose && println("degree = ", dmax)
        newsupport = exptomon(getLatticePoints(newtonPolytope(A₀pol^(dmax), x)), x)
        newshifts = [
            setdiff(
                exptomon(getLatticePoints(newtonPolytope(A₀pol^(dmax - α[i]), x)), x),
                exptomon(getLatticePoints(newtonPolytope(A₀pol^(dmax - α[i] - 1), x)), x),
            ) for i = 1:s
        ]
        newcokernel = updateCokernel(
            f,
            x,
            oldcokernel,
            oldsupport,
            newsupport,
            newshifts;
            complex = complex,
        )
        oldsupport = newsupport
        oldcokernel = newcokernel
        Σ₀ = exptomon(getLatticePoints(newtonPolytope(A₀pol^(dmax - 1), x)), x)
        chck = checkCriterion(Σα0, Σ₀, oldsupport, oldcokernel, x; rankTol = rankTol, verbose = verbose)
    end
    topdeg = dmax
    extdeg = topdeg .- α
    extdeg = [topdeg - 1; extdeg]
    Eholes = [extdeg[i] * A for i = 1:length(extdeg)]
    Dholes = topdeg * A
    Epols = [sum([prod(x .^ (AA[i, :])) for i = 1:size(AA, 1)]) for AA ∈ Eholes]
    Dpol = sum([prod(x .^ (Dholes[i, :])) for i = 1:size(Dholes, 1)])
    Epolytopes = [newtonPolytope(ff, x) for ff ∈ Epols]
    Dpolytope = newtonPolytope(Dpol, x)
    E = [getLatticePoints(P) for P ∈ Epolytopes]
    D = getLatticePoints(Dpolytope)
    A₀pol = sum([prod(x .^ (A[i, :])) for i = 1:size(A, 1)])
    A₀polytope = newtonPolytope(A₀pol, x)
    A₀ = getLatticePoints(A₀polytope)
    return oldcokernel, oldsupport, chck, dmax, (A₀, E, D)
end

# Compute the kernel of a dense Sylvester map  supported in degree DMAX by using the 'degree by degree' procedure.
function getKernel_dense_DBD(f, x, DMAX; complex = false, rankTol = 1e3 * eps(), verbose::Bool=false)
    verbose && println("Computing cokernel degree by degree...")
    s = length(f)
    ds = fill(0, s)
    for i = 1:s
        ds[i] = maximum(degree.(terms(f[i])))
    end
    dmax = maximum(ds)
    oldsupport = getMonomials(x, dmax)
    oldshifts = [getMonomials(x, dmax - ds[i]) for i = 1:s]
    Σα0 = getMonomials(x, 1)
    Σ₀ = getMonomials(x, dmax - 1)
    oldSylv = getRes(f, oldsupport, oldshifts, x; complex = complex)
    oldcokernel = transpose(nullspace(transpose(oldSylv), rtol = rankTol))
    while dmax < DMAX
        dmax = dmax + 1
        verbose && println("degree = ", dmax)
        newsupport = getMonomials(x, dmax)
        newshifts = [getMonomials(x, dmax - ds[i]; homogeneous = true) for i = 1:s]
        newcokernel = updateCokernel(
            f,
            x,
            oldcokernel,
            oldsupport,
            newsupport,
            newshifts;
            complex = complex,
        )
        testshifts = [getMonomials(x, dmax - ds[i]; homogeneous = false) for i = 1:s]
        oldsupport = newsupport
        oldcokernel = newcokernel
        Σ₀ = getMonomials(x, dmax - 1)
    end
    return oldcokernel, oldsupport, dmax
end

# Compute an admissible tuple for a (square,) dense system of equations.
# This is based on the Macaulay bound for the regularity.
function get_AT_CI_dense(f, x)
    s = length(f)
    n = length(x)
    ds = fill(0, s)
    for i = 1:s
        ds[i] = maximum(degree.(terms(f[i])))
    end
    dreg = sum(ds) - n + 1
    Dmons = getMonomials(x, dreg)
    Emons = [getMonomials(x, dreg - ds[i]) for i = 1:s]
    Emons = [[getMonomials(x, dreg - 1)]; Emons]
    A₀mons = getMonomials(x, 1)
    D = montoexp(Dmons)
    E = [montoexp(EE) for EE ∈ Emons]
    A₀ = montoexp(A₀mons)
    return A₀, E, D
end

# Compute an admissible tuple for a (square,) multi-graded dense system of equations.
# This is based on the Macaulay bound for the regularity.
function get_AT_CI_multi_dense(f, vargroups, ds)
    s = size(ds, 1)
    dreg = sum(ds, dims = 1) - length.(vargroups)' .+ 1
    Dmons = getMultiMonomials(vargroups, dreg)
    Emons = [getMultiMonomials(vargroups, dreg - ds[i, :]') for i = 1:s]
    Emons = [[getMultiMonomials(vargroups, dreg - ones(Int, length(dreg))')]; Emons]
    A₀mons = getMultiMonomials(vargroups, ones(Int, length(dreg))')
    D = montoexp(Dmons)
    E = [montoexp(EE) for EE ∈ Emons]
    A₀ = montoexp(A₀mons)
    return A₀, E, D
end

# Compute an admissible tuple for a (square,) multi-unmixed system of equations.
# This is based on the Macaulay bound for the regularity.
function get_AT_CI_multi_unmixed(f, vargroups, sups, ds)
    s = size(ds, 1)
    codegs = [getCodegree(sups[i]) for i = 1:length(vargroups)]
    dreg = sum(ds, dims = 1) - codegs' .+ 2
    Dmons = getMultiUnmixedMonomials(vargroups, sups, dreg)
    Emons = [getMultiUnmixedMonomials(vargroups, sups, dreg - ds[i, :]') for i = 1:s]
    Emons = [[getMultiUnmixedMonomials(vargroups, sups, dreg - ones(Int, length(dreg))')]; Emons]
    A₀mons = getMultiUnmixedMonomials(vargroups, sups, ones(Int, length(dreg))')
    D = montoexp(Dmons)
    E = [montoexp(EE) for EE ∈ Emons]
    A₀ = montoexp(A₀mons)
    return A₀, E, D
end

# Solve a square, dense system of equations.
# This uses the DBD method to compute the cokernel if DBD = true.
function solve_CI_dense(f, x; DBD = true, complex = false, check_criterion = false, verbose::Bool=false)
    if DBD
        n = length(x)
        s = n
        ds = fill(0, s)
        for i = 1:s
            ds[i] = maximum(degree.(terms(f[i])))
        end
        dreg = sum(ds) - n + 1
        cokernel, support, topdeg = getKernel_dense_DBD(f, x, dreg; complex = true, verbose = verbose)
        A₀ = montoexp(getMonomials(x, 1))
        extdeg = topdeg .- ds
        extdeg = [topdeg - 1; extdeg]
        E = [montoexp(getMonomials(x, extdeg[i])) for i = 1:length(f)+1]
        D = montoexp(getMonomials(x, topdeg))
        solve_EV(
            f,
            x,
            A₀,
            E,
            D;
            complex = complex,
            N = cokernel,
            check_criterion = check_criterion,
            schurfac = true,
            verbose = verbose,
        )
    else
        (A₀, E, D) = get_AT_CI_dense(f, x)
        solve_EV(
            f,
            x,
            A₀,
            E,
            D;
            complex = complex,
            check_criterion = check_criterion,
            schurfac = true,
            verbose = verbose,
        )
    end
end

# Solve a square, multi-graded dense system of equations.
function solve_CI_multi_dense(f, vargroups, ds; complex = false, check_criterion = false, verbose::Bool=false)
    (A₀, E, D) = get_AT_CI_multi_dense(f, vargroups,ds)
    solve_EV(
        f,
        vcat(vargroups...),
        A₀,
        E,
        D;
        complex = complex,
        check_criterion = check_criterion,
        schurfac = true,
        verbose = verbose,
    )
end

# Solve a square, multi-unmixed system of equations.
function solve_CI_multi_unmixed(f, vargroups, sups, ds; AT = nothing, complex = false, check_criterion = false, verbose::Bool=false)
    if isnothing(AT)
        (A₀, E, D) = get_AT_CI_multi_unmixed(f, vargroups, sups, ds)
    else
        A₀ = AT[1]
        E = AT[2]
        D = AT[3]
    end
    sol = solve_EV(
        f,
        vcat(vargroups...),
        A₀,
        E,
        D;
        complex = complex,
        check_criterion = check_criterion,
        schurfac = true,
        verbose = verbose
    )
    return sol, A₀, E, D
end

# Generate a random dense system of equations of degrees specified by ds.
function getRandomSystem_dense(x, ds; complex = false)
    n = length(x)
    s = length(ds)
    M = [getMonomials(x, d) for d ∈ ds]
    if complex
        f = [randn(ComplexF64, length(MM))' * MM for MM ∈ M]
    else
        f = [randn(length(MM))' * MM for MM ∈ M]
    end
    return f[:]
end

# Generate a random multi-graded dense system of equations. The r variable groups
# are given by vargroups, and the multidegrees in the (s x r) array ds.
function getRandomSystem_multi_dense(vargroups, ds; complex = false)
    M = [getMultiMonomials(vargroups, ds[i, :]) for i = 1:size(ds, 1)]
    if complex
        f = [randn(ComplexF64, length(MM))' * MM for MM ∈ M]
    else
        f = [randn(length(MM))' * MM for MM ∈ M]
    end
    return f[:]
end

# Generate a random multi-unmixed system of equations. The r variable groups
# are given by vargroups, the corresponding polytopes are defined by the supports in sups and the multidegrees in the (s x r) array ds.
function getRandomSystem_multi_unmixed(vargroups, sups, ds; complex = false)
    M = [getMultiUnmixedMonomials(vargroups,sups,ds[i, :]) for i = 1:size(ds, 1)]
    if complex
        f = [randn(ComplexF64, length(MM))' * MM for MM ∈ M]
    else
        f = [randn(length(MM))' * MM for MM ∈ M]
    end
    return f[:]
end

# Compute all monomials of multidegree ds in the groups of variables in vargroups.
function getMultiMonomials(vargroups, ds)
    aux = 1
    for i = 1:length(vargroups)
        aux = aux * (1 + sum(vargroups[i]))^(ds[i])
    end
    monomials(aux)
end

# Compute all monomials of multidegree ds in the groups of variables in vargroups with respect to the supports in sups.
function getMultiUnmixedMonomials(vargroups, sups, ds)
    aux = 1
    for i = 1:length(vargroups)
        aux = aux * (sum(exptomon(sups[i],vargroups[i])))^(ds[i])
    end
    P = newtonPolytope(aux,vcat(vargroups...))
    latpts = getLatticePoints(P)
    exptomon(latpts,vcat(vargroups...))
end

# Compute an admissible tuple for a square, sparse polynomial system which is generic with respect to its support,
# in the sense of "Toric Eigenvalue Methods for Solving Sparse Polynomial Systems".
function get_AT_CI_mixed(f, x)
    NP = [newtonPolytope(ff, x) for ff ∈ f]
    NP0 = newtonPolytope(1 + sum(x), x)
    NPS = [NP0; NP]
    Minksum, vtxmtx = sumAll2(NPS)
    D = getLatticePoints(Minksum)
    inds = [setdiff(1:length(NPS), i) for i = 1:length(NPS)]
    NPi = [sumAll2(NPS[ii])[1] for ii ∈ inds]
    E = [getLatticePoints(Δ) for Δ ∈ NPi]
    A₀ = getLatticePoints(NP0)
    return A₀, E, D
end

# Solve a square, sparse polynomial system which is generic with respect to its support,
# in the sense of "Toric Eigenvalue Methods for Solving Sparse Polynomial Systems".
function solve_CI_mixed(f, x; AT = nothing, complex = false, verbose::Bool=false)
    if isnothing(AT)
        (A₀, E, D) = get_AT_CI_mixed(f, x)
    else
        A₀ = AT[1]
        E = AT[2]
        D = AT[3]
    end
    sol = solve_EV(f, x, A₀, E, D; complex = complex, verbose = verbose)
    return sol, A₀, E, D
end

# Compute an admissible tuple for a square, sparse, unmixed polynomial system which is generic with respect to its support,
# in the sense of "Toric Eigenvalue Methods for Solving Sparse Polynomial Systems".
function get_AT_CI_unmixed(x, A, α)
    α = [1; α]
    r = getCodegree(A)
    inds = [setdiff(1:length(α), i) for i = 1:length(α)]
    Eholes = [(sum(α[ii]) - r + 1) * A for ii ∈ inds]
    Dholes = (sum(α) - r + 1) * A
    Epols = [sum([prod(x .^ (AA[i, :])) for i = 1:size(AA, 1)]) for AA ∈ Eholes]
    Dpol = sum([prod(x .^ (Dholes[i, :])) for i = 1:size(Dholes, 1)])
    Epolytopes = [newtonPolytope(ff, x) for ff ∈ Epols]
    Dpolytope = newtonPolytope(Dpol, x)
    E = [getLatticePoints(P) for P ∈ Epolytopes]
    D = getLatticePoints(Dpolytope)
    A₀pol = sum([prod(x .^ (A[i, :])) for i = 1:size(A, 1)])
    A₀polytope = newtonPolytope(A₀pol, x)
    A₀ = getLatticePoints(A₀polytope)
    return A₀, E, D
end

# Generate a random sparse unmixed system such that fᵢ has support in the convex hull of αᵢA.
function getRandomSystem_unmixed(x, A, α; complex = false)
    n = length(x)
    s = length(α)
    monsholes = [exptomon(a * A, x) for a ∈ α]
    pols = [sum(mons) for mons ∈ monsholes]
    exps = [getLatticePoints(newtonPolytope(ff, x)) for ff ∈ pols]
    M = [exptomon(e, x) for e ∈ exps]
    if complex
        f = [randn(ComplexF64, length(MM))' * MM for MM ∈ M]
    else
        f = [randn(length(MM))' * MM for MM ∈ M]
    end
    return f[:]
end

# Solve a square, sparse, unmixed polynomial system which is generic with respect to its support,
# in the sense of "Toric Eigenvalue Methods for Solving Sparse Polynomial Systems".
# !!!!!!!!!!!!!!!!!!! MAKE SURE A CONTAINS 0
function solve_CI_unmixed(f, x, A, α; AT = nothing, verbose::Bool=false)
    if isnothing(AT)
        (A₀, E, D) = get_AT_CI_unmixed(x, A, α)
    else
        A₀ = AT[1]
        E = AT[2]
        D = AT[3]
    end
    sol = solve_EV(f, x, A₀, E, D; verbose = verbose)
    return sol, A₀, E, D
end

# Solve an overdetermined system of dense equations.
function solve_OD_dense(f, x; maxdeg = 100, complex = false, verbose::Bool=false, trySemiregular::Bool=false)
    verbose && println("Looking for an admissible tuple")
    cokernel, support, chck, topdeg, AT =
        findAdmissibleTuple_dense(f, x, maxdeg; complex = true, verbose = verbose, trySemiregular = trySemiregular)
    (A₀, E, D) = AT
    sol = solve_EV(f, x, A₀, E, D; N = cokernel, complex = complex, check_criterion = false, verbose = verbose)
    return sol
end

# Solve an overdetermined system of sparse unmixed equations.
function solve_OD_unmixed(f, x, A, α; maxdeg = 100, complex = false, verbose::Bool=false)
    verbose && println("Looking for an admissible tuple")
    cokernel, support, chck, topdeg, AT =
        findAdmissibleTuple_unmixed(f, x, A, α, maxdeg; complex = true, verbose = verbose)
    (A₀, E, D) = AT
    sol = solve_EV(f, x, A₀, E, D; N = cokernel, complex = complex, check_criterion = false, verbose = verbose)
    return sol
end

# solve vA = λvB for m x n matrices A, B with m <= n.
function solveREP(A, B; rtol = 1e-10)
    (m, n) = size(A) # should be the same as size(B)
    if m == 1
        svdobj = svd([A; B])
        if svdobj.S[2] < rtol * svdobj.S[1]
            return [mean(A ./ B)], [[1]]
        else
            return [], []
        end
    else
        O = randn(n, m)
        qrobj = qr(O)
        O = qrobj.Q[:, 1:m]
        eigenobj = eigen(transpose(A * O), transpose(B * O))
        Λ = eigenobj.values
        eigenvals = fill(0.0 + 0.0 * im, 0)
        eigenspaces = fill(fill(0.0 + 0.0 * im, m, 1), 0)
        for λ ∈ Λ
            ns = nullspace(transpose(A - λ * B), rtol = rtol)
            if size(ns, 2) > 0
                push!(eigenvals, λ)
                push!(eigenspaces, ns')
            end
        end
        return eigenvals, eigenspaces
    end
end

# Return all monomials of degree ≤ d in x (equality if homogeneous = true).
function getMonomials(x, d; homogeneous = false)
    if homogeneous
        return monomials(sum(x)^d)
    else
        return monomials((1 + sum(x))^d)
    end
end

# Return a basis for the space of vanishing polynomials in x supported on mons and vanishing on the points in pts.
function getVanishingPolynomials(pts, mons, x; augm_prec = false)
    if augm_prec
        pts = convert.(Array{Complex{BigFloat},1}, pts)
    end
    exps = montoexp(mons)
    Vdm = [[prod(pt .^ (exps[i, :])) for i = 1:size(exps, 1)] for pt ∈ pts]
    Vdm = [v / norm(v) for v ∈ Vdm]
    Vdm = hcat(Vdm...)
    N = nullspace(transpose(Vdm))
    N = convert.(Complex{Float64}, N)
    κ = cond(Vdm)
    return transpose(N) * mons, κ
end

# Return the degrees of the polynomials in f.
function getDegrees(f)
    d = fill(0, s)
    for i = 1:s
        d[i] = maximum(degree.(terms(f[i])))
    end
    return d
end

# This computes the relative backward errors of the points in sol for the system f in the variables x.
# It uses the definition in Appendix C of my thesis.
function get_residual(f, sol, x)
    n = length(x)
    residuals = zeros(length(sol))
    for k = 1:length(sol)
        ressol = 0
        for i = 1:length(f)
            R = f[i]
            T = terms(R)
            l = length(T)
            Rabs = sum(abs.([T[s](x => abs.(sol[k])) for s = 1:l]))
            R = R(x => sol[k]) / (1 + Rabs)
            ressol = ressol + abs(R)
        end
        residuals[k] = ressol / length(f)
    end
    return residuals
end

# simultaneous diagonalization of the matrices in M via Schur factorization of a random linear combination.
# returns a matrix with the eigenvalues (`third mode` of the CPD)
function simulDiag_schur_mtx(M; clustering = false, tol = 1e-4, complex = false)
    n = length(M)
    δ = size(M[1], 1)
    Mₕ = sum(randn(n) .* M)
    F = schur(Mₕ)
    if complex == false
        F = triangularize(F)
    end
    solMatrix = fill(0.0 + 0.0 * im, δ, n)
    if clustering
        oF = ordschur(F, fill(true, δ))
        v = diag(oF.Schur)
        clusvec = getClusters(v; tol = tol)
        push!(clusvec, δ + 1)
        Q = oF.Z
        for i = 1:n
            Uᵢ = Q' * M[i] * Q
            eᵢ = fill(zero(ComplexF64), δ)
            for j = 1:length(clusvec)-1
                i₀ = clusvec[j]
                i₁ = clusvec[j+1] - 1
                ev = tr(Uᵢ[i₀:i₁, i₀:i₁]) / (i₁ - i₀ + 1)
                eᵢ[i₀:i₁] .= ev
            end
            solMatrix[:, i] = eᵢ
        end
    else
        Q = F.Z
        for i = 1:n
            Uᵢ = Q' * M[i] * Q
            if i == 1 && norm(Uᵢ - triu(Uᵢ)) > 1e-10 * norm(Uᵢ)
                println("there may be singular solutions, try the option clustering = true")
            end
            solMatrix[:, i] = diag(Uᵢ)
        end
    end
    return solMatrix
end

# returns a vector clusvec with the indices of the start of each cluster in v.
# v is a vector of complex number which is assumed to be ordered, e.g. coming from an ordered Schur factorization
function getClusters(v; tol = 1e-4)
    clus = [v[1]]
    clusvec = [1]
    k = 2
    while k <= length(v)
        center = mean(clus)
        if abs(v[k] - center) < tol * (abs(center) + 1)
            push!(clus, v[k])
        else
            push!(clusvec, k)
            clus = [v[k]]
        end
        k = k + 1
    end
    return clusvec
end

end

max_stars_repo_name: simontelen/JuliaEigenvalueSolver

max_issues_repo_name: simontelen/JuliaEigenvalueSolver

max_forks_repo_name: simontelen/JuliaEigenvalueSolver


text:
@testset "Embedding operators" begin

@testset "RandomBound" begin
    @testset "does nothing if within bounds" begin
        @test apply!(RandomBound([(0.0, 1.0)]), [0.0], [0.0]) == [0.0]

        @test apply!(RandomBound([(0.0, 1.0), (10.0, 15.0)]),
                     [0.0, 11.4], [0.1, 12.3] ) == [0.0, 11.4]
    end

    @testset "bounds if lower than min bound" begin
        @test apply!(RandomBound([(0.0, 1.0)]), [-0.1], [0.0]) == [0.0]

        res = apply!(RandomBound([(0.0, 1.0)]), [-0.1], [0.5])
        @test res[1] >= 0.0
        @test res[1] <= 0.5

        res = apply!(RandomBound([(-11.0, 1.0), (0.0, 1.0)]),
                                                            [-11.1, 0.5], [-10.8, 0.5])
        @test (-10.8 >= res[1] >= -11.0)
        @test res[2] == 0.5

        res = apply!(RandomBound([(30.0, 60.0), (-102.0, -1.0)]),
                     [50.4, -103.1], [49.6, -101.4])
        @test res[1] == 50.4
        @test (-101.4 >= res[2] >= -102.0)
    end

    @testset "bounds if higher than max bound" begin
        @test apply!(RandomBound([(0.0, 1.0)]), [1.1], [1.0]) == [1.0]

        res = apply!(RandomBound([(-10.0, 96.0)]), [97.0], [95.0])
        @test (95.0 <= res[1] <= 96.0)
    end
end

end