text:
module AMQPTestCoverage

using AMQPClient, Test, Random

const JULIA_HOME = Sys.BINDIR

const EXCG_DIRECT = "ExcgDirect"
const EXCG_FANOUT = "ExcgFanout"
const QUEUE1 = "queue1"
const ROUTE1 = "key1"
const invalid_auth_params = Dict{String,Any}("MECHANISM"=>"AMQPLAIN", "LOGIN"=>randstring(10), "PASSWORD"=>randstring(10))

function runtests(;virtualhost="/", host="localhost", port=AMQPClient.AMQP_DEFAULT_PORT, auth_params=AMQPClient.DEFAULT_AUTH_PARAMS, amqps=nothing, keepalive=true, heartbeat=true)
    verify_spec()
    test_types()
    @test default_exchange_name("direct") == "amq.direct"
    @test default_exchange_name() == ""
    @test AMQPClient.method_name(AMQPClient.TAMQPMethodPayload(:Basic, :Ack, (1, false))) == "Basic.Ack"

    # test failure on invalid auth_params
    @test_throws AMQPClient.AMQPClientException connection(;virtualhost=virtualhost, host=host, port=port, amqps=amqps, auth_params=invalid_auth_params)

    conn_ref = nothing

    # open a connection
    @info("opening connection")
    connection(;virtualhost=virtualhost, host=host, port=port, amqps=amqps, auth_params=auth_params, send_queue_size=512, keepalive=keepalive, heartbeat=heartbeat) do conn
        # Issue #51
        @test isa(conn, AMQPClient.Connection)

        @test conn.sendq.sz_max == 512

        # open a channel
        @info("opening channel")
        channel(conn, AMQPClient.UNUSED_CHANNEL, true) do chan1
            @test chan1.id == 1
            @test conn.sendq.sz_max == 512

            # test default exchange names
            @test default_exchange_name() == ""
            @test default_exchange_name(EXCHANGE_TYPE_DIRECT) == "amq.direct"

            # create exchanges
            @info("creating exchanges")
            @test exchange_declare(chan1, EXCG_DIRECT, EXCHANGE_TYPE_DIRECT; arguments=Dict{String,Any}("Hello"=>"World", "Foo"=>"bar"))
            @test exchange_declare(chan1, EXCG_FANOUT, EXCHANGE_TYPE_FANOUT)
            # redeclaring the exchange with same attributes should be fine
            @test exchange_declare(chan1, EXCG_FANOUT, EXCHANGE_TYPE_FANOUT)
            # redeclaring an existing exchange with different attributes should fail
            @test_throws AMQPClient.AMQPClientException exchange_declare(chan1, EXCG_FANOUT, EXCHANGE_TYPE_DIRECT)
        end

        chan_ref = nothing
        # must reconnect as channel gets closed after a channel exception
        channel(conn, AMQPClient.UNUSED_CHANNEL, true) do chan1
            @test chan1.id == 1

            # create and bind queues
            @info("creating queues")
            success, queue_name, message_count, consumer_count = queue_declare(chan1, QUEUE1)
            @test success
            @test message_count == 0
            @test consumer_count == 0

            @test queue_bind(chan1, QUEUE1, EXCG_DIRECT, ROUTE1)

            # rabbitmq 3.6.5 does not support qos
            # basic_qos(chan1, 1024*10, 10, false)

            M = Message(Vector{UInt8}("hello world"), content_type="text/plain", delivery_mode=PERSISTENT)

            @info("testing basic publish and get")
            # publish 10 messages
            for idx in 1:10
                basic_publish(chan1, M; exchange=EXCG_DIRECT, routing_key=ROUTE1)
                flush(chan1)
                @test !isready(chan1.conn.sendq)
            end

            # basic get 10 messages
            for idx in 1:10
                result = basic_get(chan1, QUEUE1, false)
                @test result !== nothing
                rcvd_msg = result
                basic_ack(chan1, rcvd_msg.delivery_tag)
                @test rcvd_msg.remaining == (10-idx)
                @test rcvd_msg.exchange == EXCG_DIRECT
                @test rcvd_msg.redelivered == false
                @test rcvd_msg.routing_key == ROUTE1
                @test rcvd_msg.data == M.data
                @test :content_type in keys(rcvd_msg.properties)
                @test convert(String, rcvd_msg.properties[:content_type]) == "text/plain"
            end

            # basic get returns null if no more messages
            @test basic_get(chan1, QUEUE1, false) === nothing

            ## test reject and requeue
            basic_publish(chan1, M; exchange=EXCG_DIRECT, routing_key=ROUTE1)

            result = basic_get(chan1, QUEUE1, false)
            @test result !== nothing
            rcvd_msg = result
            @test rcvd_msg.redelivered == false

            basic_reject(chan1, rcvd_msg.delivery_tag; requeue=true)

            result = basic_get(chan1, QUEUE1, false)
            @test result !== nothing
            rcvd_msg = result
            @test rcvd_msg.redelivered == true

            basic_ack(chan1, rcvd_msg.delivery_tag)

            @info("testing basic consumer")
            # start a consumer task
            global msg_count = 0
            consumer_fn = (rcvd_msg) -> begin
                @test rcvd_msg.exchange == EXCG_DIRECT
                @test rcvd_msg.redelivered == false
                @test rcvd_msg.routing_key == ROUTE1
                global msg_count
                msg_count += 1
                if msg_count <= 10
                    @test rcvd_msg.data == M.data
                else
                    @test rcvd_msg.data == UInt8[]
                end
                println("received msg $(msg_count): $(String(rcvd_msg.data))")
                basic_ack(chan1, rcvd_msg.delivery_tag)
            end
            success, consumer_tag = basic_consume(chan1, QUEUE1, consumer_fn)
            @test success

            # publish 10 messages
            for idx in 1:10
                basic_publish(chan1, M; exchange=EXCG_DIRECT, routing_key=ROUTE1)
            end

            # wait for a reasonable time to receive all messages
            for idx in 1:10
                (msg_count == 10) && break
                sleep(1)
            end
            @test msg_count == 10

            @info("testing empty messages")
            # Test sending and receiving empty message
            M_empty = Message(Vector{UInt8}(), content_type="text/plain", delivery_mode=PERSISTENT)
            basic_publish(chan1, M_empty; exchange=EXCG_DIRECT, routing_key=ROUTE1)

            M_no_ct = Message(Vector{UInt8}(), delivery_mode=PERSISTENT)
            basic_publish(chan1, M_no_ct; exchange=EXCG_DIRECT, routing_key=ROUTE1)

            println("Waiting")
            # wait for a reasonable time to receive last two messages
            for idx in 1:5
                (msg_count == 12) && break
                sleep(1)
            end
            println("Waited")
            @test msg_count == 12

            # cancel the consumer task
            @test basic_cancel(chan1, consumer_tag)

            # test transactions
            @info("testing tx")
            @test tx_select(chan1)
            @test tx_commit(chan1)
            @test tx_rollback(chan1)

            # test heartbeats
            if 120 >= conn.heartbeat > 0
                c = conn
                @info("testing heartbeats (waiting $(3*c.heartbeat) secs)...")
                ts1 = c.heartbeat_time_server
                tc1 = c.heartbeat_time_client
                sleeptime = c.heartbeat/2
                for idx in 1:6
                    (c.heartbeat_time_server > ts1) && (c.heartbeat_time_client > tc1) && break
                    sleep(sleeptime)
                end
                @test c.heartbeat_time_server > ts1
                @test c.heartbeat_time_client > tc1
            elseif conn.heartbeat == 0
                @info("heartbeat disabled")
            else
                @info("not testing heartbeats (wait too long at $(3*conn.heartbeat) secs)")
            end

            @info("closing down")
            success, message_count = queue_purge(chan1, QUEUE1)
            @test success
            @test message_count == 0

            @test queue_unbind(chan1, QUEUE1, EXCG_DIRECT, ROUTE1)
            success, message_count = queue_delete(chan1, QUEUE1)
            @test success
            @test message_count == 0

            # delete exchanges
            @test exchange_delete(chan1, EXCG_DIRECT; nowait=true)
            @test exchange_delete(chan1, EXCG_FANOUT)

            chan_ref = chan1 # to do additional tests on a closed channel
        end

        close(chan_ref) # closing a closed channel should not be an issue
        AMQPClient.wait_for_state(chan_ref, AMQPClient.CONN_STATE_CLOSED)
        @test !isopen(chan_ref)

        conn_ref = conn  # to do additional tests on a closed connection
    end

    # closing a closed connection should not be an issue
    close(conn_ref)
    AMQPClient.wait_for_state(conn_ref, AMQPClient.CONN_STATE_CLOSED)
    @test !isopen(conn_ref)

    @info("done")
    nothing
end

function verify_spec()
    ALLCLASSES = (:Connection, :Basic, :Channel, :Confirm, :Exchange, :Queue, :Tx)
    for n in ALLCLASSES
        @test n in keys(AMQPClient.CLASSNAME_MAP)
    end
    for (n,v) in keys(AMQPClient.CLASSMETHODNAME_MAP)
        @test n in ALLCLASSES
    end
end

function test_types()
    d = Dict{String,Any}(
        "bool"      => 0x1,
        "int"       => 10,
        "uint"      => 0x1,
        "float"     => rand(),
        "shortstr"  => AMQPClient.TAMQPShortStr(randstring(10)),
        "longstr"   => AMQPClient.TAMQPLongStr(randstring(1024)))
    ft = AMQPClient.TAMQPFieldTable(d)
    iob = IOBuffer()
    show(iob, ft)
    @test length(take!(iob)) > 0

    fields = [Pair{Symbol,AMQPClient.TAMQPField}(:bit,          AMQPClient.TAMQPBit(0x1)),
              Pair{Symbol,AMQPClient.TAMQPField}(:shortstr,     AMQPClient.TAMQPShortStr(randstring(10))),
              Pair{Symbol,AMQPClient.TAMQPField}(:longstr,      AMQPClient.TAMQPLongStr(randstring(1024))),
              Pair{Symbol,AMQPClient.TAMQPField}(:fieldtable,   ft)]
    show(iob, fields)
    @test length(take!(iob)) > 0

    mpayload = AMQPClient.TAMQPMethodPayload(:Channel, :Open, ("",))
    show(iob, mpayload)
    @test length(take!(iob)) > 0

    mfprop = AMQPClient.TAMQPFrameProperties(AMQPClient.TAMQPChannel(0), AMQPClient.TAMQPPayloadSize(100))
    show(iob, mfprop)
    @test length(take!(iob)) > 0

    mframe = AMQPClient.TAMQPMethodFrame(mfprop, mpayload)
    show(iob, mframe)
    @test length(take!(iob)) > 0

    fields = AMQPClient.TAMQPFieldValue[
        AMQPClient.TAMQPFieldValue(true),
        AMQPClient.TAMQPFieldValue(1.1),
        AMQPClient.TAMQPFieldValue(1),
        AMQPClient.TAMQPFieldValue("hello world"),
        AMQPClient.TAMQPFieldValue(Dict{String,Int}("one"=>1, "two"=>2)),
    ]

    fieldarray = AMQPClient.TAMQPFieldArray(fields)
    simplified_fields = AMQPClient.simplify(fieldarray)
    @test simplified_fields == Any[
        0x01,
        1.1,
        1,
         "hello world",
         Dict{String, Any}("two" => 2, "one" => 1)
    ]

    iob = PipeBuffer()
    write(iob, hton(AMQPClient.TAMQPLongUInt(10)))
    write(iob, UInt8[1,2,3,4,5,6,7,8,9,0])
    barr = read(iob, AMQPClient.TAMQPByteArray)
    @test barr.len == 10
    @test barr.data == UInt8[1,2,3,4,5,6,7,8,9,0]
end

end # module AMQPTestCoverage


text:
function _make_reference_gui()
    ref_grid=Grid()

    ref_list1=ListStore(Int32)
    for i=1:0
        push!(ref_list1,(i,))
    end
    ref_tv1=TreeView(TreeModel(ref_list1))
    ref_r1=CellRendererText()
    ref_c1=TreeViewColumn("New Reference:",ref_r1, Dict([("text",0)]))

    ref_tv1_s=Gtk.GAccessor.selection(ref_tv1)

    push!(ref_tv1,ref_c1)

    ref_scroll1=ScrolledWindow()
    Gtk.GAccessor.min_content_height(ref_scroll1,350)
    Gtk.GAccessor.min_content_width(ref_scroll1,175)
    push!(ref_scroll1,ref_tv1)

    ref_list2=ListStore(Int32)
    for i=1:0
        push!(ref_list2,(i,))
    end
    ref_tv2=TreeView(TreeModel(ref_list2))
    ref_r2=CellRendererText()
    ref_c2=TreeViewColumn("Apply Reference To:",ref_r2, Dict([("text",0)]))

    ref_tv2_s=Gtk.GAccessor.selection(ref_tv2)
    Gtk.GAccessor.mode(ref_tv2_s,Gtk.GConstants.GtkSelectionMode.MULTIPLE)

    push!(ref_tv2,ref_c2)

    ref_scroll2=ScrolledWindow()
    Gtk.GAccessor.min_content_height(ref_scroll2,350)
    Gtk.GAccessor.min_content_width(ref_scroll2,175)
    push!(ref_scroll2,ref_tv2)

    ref_button2=Button("Select All/None")
    ref_button3=Button("Apply")

    ref_grid[1,1]=ref_scroll1
    ref_grid[2,1]=Canvas(50,350)
    ref_grid[3,1]=ref_scroll2
    ref_grid[3,2]=ref_button2
    ref_grid[2,3]=Canvas(50,50)
    ref_grid[3,4]=ref_button3
    ref_win=Window(ref_grid)
    setproperty!(ref_win, :title, "Reference Channel Select")

    Gtk.showall(ref_win)
    visible(ref_win,false)

    ref_widgets = Reference_Widgets(ref_win,ref_tv1,ref_tv2,ref_list1,ref_list2,ref_button2,ref_button3)
end

function add_reference_callbacks(w,handles,r,fpga)

    signal_connect(ref_b2_cb, w.select_button, "clicked",Void,(),false,(handles,))
    signal_connect(ref_b3_cb, w.apply_button, "clicked",Void,(),false,(handles,r))

    signal_connect(w.win, :delete_event) do widget, event
        visible(w.win, false)
        true
    end
end

function update_ref(rhd::RHD2000,han::Gui_Handles)

    selmodel_l=Gtk.GAccessor.selection(han.ref_widgets.tv1)
    selmodel_r=Gtk.GAccessor.selection(han.ref_widgets.tv2)

    for i=1:size(rhd.v,2)

        myref=rhd.refs[i]

        if myref!=0

            select!(selmodel_l,Gtk.iter_from_index(han.ref_widgets.list1,myref))
            select!(selmodel_r,Gtk.iter_from_index(han.ref_widgets.list2,i))
        end
    end
    nothing
end

#=
Callbacks for reference channel selection
=#

function ref_cb(widget::Ptr,user_data::Tuple{Gui_Handles,RHD2000})

    han, rhd = user_data

    visible(han.ref_widgets.win,true)
    nothing
end

function ref_b1_cb(widget::Ptr, user_data::Tuple{Gui_Handles})

    han, = user_data
    selmodel=Gtk.GAccessor.selection(han.ref_widgets.tv1)
    selectall!(selmodel)

    nothing
end

function ref_b2_cb(widget::Ptr, user_data::Tuple{Gui_Handles})

    han, = user_data
    selmodel=Gtk.GAccessor.selection(han.ref_widgets.tv2)
    selectall!(selmodel)

    nothing
end

function ref_b3_cb(widget::Ptr, user_data::Tuple{Gui_Handles,RHD2000})

    han, rhd = user_data

    #Left
    selmodel=Gtk.GAccessor.selection(han.ref_widgets.tv1)
    myref=parse(Int64,Gtk.get_string_from_iter(TreeModel(han.ref_widgets.list1), Gtk.selected(selmodel)))+1

    #Right
    for i=1:size(rhd.v,2)
        if is_selected(han.ref_widgets.list2,han.ref_widgets.tv2,i-1)
            if i==myref
                rhd.refs[i]=0
            else
                rhd.refs[i]=myref
            end
        end
    end

    f=open(string(rhd.save.backup,"ref.bin"),"w")
    write(f,rhd.refs)
    close(f)

    nothing
end

ext: jl


text:
test_prior_failure_env = {
    "cycle_length" => 12,
    "coverage"=> [3,4,5,5,5,5,5,5,4,4,2,2]
}

test_prior_failure_employees = {
    "Jaeger"=>{
        "max"=>9,
        "min"=>3,
        "availability"=>[1,1,1,1,1,1,1,1,1,1,1,1]
    },
    "Damian"=>{
        "max"=>8,
        "min"=>3,
        "availability"=>[1,1,1,1,1,1,1,1,1,1,1,1]
    },
    "Patrick"=>{
        "max"=>8,
        "min"=>3,
        "availability"=>[1,1,1,1,1,1,1,1,1,1,1,1]
    },
    "Rufus"=>{
        "max"=>8,
        "min"=>3,
        "availability"=>[1,1,1,1,1,1,1,1,1,1,1,1]
    },
    "Shaun"=>{
        "max"=>8,
        "min"=>3,
        "availability"=>[1,1,1,1,1,1,1,1,1,1,1,1]
    },
    "Douglas"=>{
        "max"=>8,
        "min"=>3,
        "availability"=>[1,1,1,1,1,1,1,1,1,1,1,1]
    },
    "Leonardo"=>{
        "max"=>8,
        "min"=>3,
        #"availability"=>[1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,0.0]
        "availability"=>[1,1,1,1,1,1,1,1,1,1,1,0],
    },
}

max_stars_repo_name: shdann/autoshift

max_issues_repo_name: shdann/autoshift

max_forks_repo_name: shdann/autoshift


text:
function readmpdata(fname,data)
open(fname,"r") do io
read!(io,data)
end
end

function writempdata(fname,data)
open(fname,"w") do io
write(io, data)
end
end

function Fname(c,T)
if T==Float64
   fname=string("F3264c",string(c))
else
   fname=string("F16c",string(c))
end
return fname
end

max_stars_repo_name: ctkelley/MPResults


text:
using ForwardDiff
using LinearAlgebra: mul!
using Parameters: @with_kw, @unpack
using Setfield: set
using StaticArrays: SVector, SMatrix, Size

using QuadGK: gauss
using ..PolyUtils: lagrange, dlagrange

import ..Continuations: get_prob_cache, get_u0, residual!, residual_jacobian!,
    isindomain
using ..BifurcationsBase: MutableState, TimeKind

abstract type
    Codim1LCProblem{skind, tkind} <: BifurcationProblem{skind, tkind}
end

struct LimitCycleProblem{
        skind <: MutableState,  # <: StateKind,
        tkind <: TimeKind,
        XS <: AbstractMatrix,
        X,
        R <: Real,
        P,
        L,
        } <: Codim1LCProblem{skind, tkind}
    statekind::skind
    timekind::tkind
    xs0::XS
    l0::R  # period
    t0::R  # parameter  (TODO: rename parameter to something else)
    period_bound::Tuple{R, R}
    t_domain::Tuple{R, R}
    phase_space::Tuple{X, X}
    diameter_bound::Tuple{R, R}
    time_offset::R
    de_prob::P
    param_axis::L
    # TODO: Maybe move `num_mesh` and `degree` to solver options:
    num_mesh::Int
    degree::Int
end

function LimitCycleProblem(;
        xs0,
        l0,
        t0,
        period_bound = (0.0, Inf),
        t_domain = (typemin(typeof(t0)), typemax(typeof(t0))),
        phase_space = (typemin(eltype(xs0)),
                       typemax(eltype(xs0))),
        diameter_bound = (1e-5, Inf),
        time_offset = 0.0,
        de_prob,
        param_axis,
        num_mesh,
        degree,
        )
    return LimitCycleProblem(
        MutableState(),
        timekind(de_prob),
        xs0,
        l0,
        t0,
        period_bound,
        t_domain,
        phase_space,
        diameter_bound,
        time_offset,
        de_prob,
        param_axis,
        num_mesh,
        degree,
    )
end

dim_state(prob::LimitCycleProblem) = size(prob.xs0, 1)
degree(prob::LimitCycleProblem) = prob.degree
num_mesh(prob::LimitCycleProblem) = prob.num_mesh

# ----------------------------------------------------------------------- cache

struct LimitCycleCache{P, VRef, MV, ML, VV} <: AbstractProblemCache{P}
    prob::P
    reference::VRef
    dv::MV
    lagrange_polynomial_vals::ML
    lagrange_polynomial_driv::ML
    gauss_quadrature_weight::VV
end

# Indices:
#   i  ∈  1:m    :  collocation/evaluation/Gauss point/node index; ζᵢ
#   k  ∈  1:m+1  :  "sample" node index; x(τₖ)
# Arrays:       Index:    Size:
#   lpv, lpd    [i, k]    (m, m + 1)
#   weight      [i]       (m,)
# Abbreviations:
#   lpv     :  lagrange_polynomial_vals
#   lpd     :  lagrange_polynomial_driv
#   weight  :  gauss_quadrature_weight

function LimitCycleCache(prob)
    n = dim_state(prob)
    m = degree(prob)
    reference = reshape(copy(prob.xs0), length(prob.xs0))
    dv = similar(prob.xs0, (n, m))

    # Compute τ: equi-distant points (within a mesh)
    interval = 1 / num_mesh(prob)
    dτ = interval / m
    τ = range(0, step=dτ, length=m + 1)
    @assert τ[end] ≈ interval

    ζ_, gauss_quadrature_weight_ = gauss(m)
    ζ = (ζ_ .+ 1) ./ 2 .* interval  # Gauss nodes
    gauss_quadrature_weight = gauss_quadrature_weight_ ./ 2 .* interval
    lagrange_polynomial_vals = [lagrange(k, ζᵢ, τ) for ζᵢ in ζ, k in 1:m+1]
    lagrange_polynomial_driv = [dlagrange(k, ζᵢ, τ) for ζᵢ in ζ, k in 1:m+1]

    return LimitCycleCache(
        prob,
        reference,
        SMatrix{n, m}(dv),
        SMatrix{m, m + 1}(lagrange_polynomial_vals),
        SMatrix{m, m + 1}(lagrange_polynomial_driv),
        SVector{m}(gauss_quadrature_weight),
    )
end
# TODO: dv is only used to hold the dimension of the dynamical system.
# Remove it.

Base.@pure dim_state(cache::LimitCycleCache) = Size(cache.dv)[1]
Base.@pure degree(cache::LimitCycleCache) = Size(cache.dv)[2]
num_mesh(cache::LimitCycleCache) = cache.prob.num_mesh

function set_reference!(cache::LimitCycleCache, u)
    copyto!(cache.reference, @view u[1:length(cache.reference)])
end

# ------------------------------------------------------ continuation interface

get_prob_cache(prob::LimitCycleProblem) = LimitCycleCache(prob)

get_u0(prob::LimitCycleProblem) = vcat(
    view(prob.xs0, :),
    prob.l0,
    prob.t0,
)

u_idx_period(cache) = length(cache.prob.xs0) + 1
u_idx_param(cache::LimitCycleCache) = u_idx_period(cache) + 1
H_idx_phase_condition(cache) = length(cache.prob.xs0) + 1

isindomain(u, cache::LimitCycleCache) = isindomain_lc(u, cache)

function isindomain_lc(u, wrapper)
    cache = as(wrapper, LimitCycleCache)
    tmin, tmax = wrapper.prob.t_domain
    all(tmin .<= u[u_idx_param(wrapper)] .<= tmax) || return false

    period_min, period_max = cache.prob.period_bound
    (period_min <= u[u_idx_period(wrapper)] <= period_max) || return false

    xmin, xmax = cache.prob.phase_space
    n = dim_state(cache)
    for i in 1:degree(cache) * num_mesh(cache)
        x = @view u[n * (i - 1) + 1:n * i]
        all(xmin .<= x .<= xmax) || return false
    end

    xs = reshape((@view u[1:length(cache.prob.xs0)]), size(cache.prob.xs0))
    in_diameter_bound(xs, cache.prob.diameter_bound) || return false

    return true
end

function in_diameter_bound(xs::AbstractMatrix, diameter_bound)
    diameter_min, diameter_max = diameter_bound
    if diameter_max == typemax(diameter_max) && diameter_min == 0.0
        return true
    elseif diameter_max == typemax(diameter_max)
        return any(d -> d >= diameter_min, iter_diameter_candidates(xs))
    elseif diameter_min == 0.0
        return all(d -> d <= diameter_max, iter_diameter_candidates(xs))
    else
        d = maximum(iter_diameter_candidates(xs))
        return (diameter_min <= d <= diameter_max)
    end
end

iter_diameter_candidates(xs::AbstractMatrix) =
    (norm(@views xs[:, i] .- xs[:, j])
     for i in 1:size(xs, 2)
     for j in i + 1:size(xs, 2))

function diameter(u::AbstractVector, prob::LimitCycleProblem)
    xs = reshape((@view u[1:length(prob.xs0)]), size(prob.xs0))
    return maximum(iter_diameter_candidates(xs))
end

# ------------------------------------------------------------------- workspace

@inline function diff_range(cache, j, i)
    n = dim_state(cache)
    m = degree(cache)
    offset = n * m * (j - 1) + n * (i - 1)
    return offset + 1:offset + n
end

# TODO: get_samples --- better name than "samples"?
@generated function get_samples(
        cache::LimitCycleCache{P, VRef, MV, ML, VV},
        u, j,
        ) where {P, VRef, MV, ML, VV}
    n = Size(MV)[1]
    m = Size(MV)[2]

    idxs_rest = [:($(n * m) * (j - 1) + $i) for i in  1:(n * (m + 1))]
    idxs_last = [
        idxs_rest[1:end - n]
        1:n
    ]
    Mat = SMatrix{n, m + 1, eltype(u)}

    quote
        if j == num_mesh(cache)
            @inbounds return $Mat($([:(u[$i]) for i in idxs_last]...))
        else
            @inbounds return $Mat($([:(u[$i]) for i in idxs_rest]...))
        end
    end
end
# TODO: Stop hard-coding output type here. Assuming that m is always
# small (for StaticArrays to be effective) is probably fine but
# assuming n to be small is not OK.  But I need to optimize a lot for
# Bifurcations.jl to be useful for non-small n anyway.

# Indices:
#   i  ∈  1:m    :  collocation/evaluation/Gauss point/node index; ζᵢ
#   k  ∈  1:m+1  :  "sample" node index; x(τₖ)
#   p  ∈  1:n    :  phase space dimension
# Arrays:       Index:    Size:
#   x, dx       [p, i]    (n, m)
#   lpv, lpd    [i, k]    (m, m + 1)
#   xτ          [p, k]    (n, m + 1)

@inline function collocation(cache, u, j)
    lpv = cache.lagrange_polynomial_vals  # ℓᵀ; ℓᵢₖ = [ℓᵀ]ₖᵢ = ℓₖ(ζᵢ)
    lpd = cache.lagrange_polynomial_driv
    xτ = get_samples(cache, u, j) # xₚ(τₖ)
    x = xτ * lpv'                    # x(ζ) = x(τ) * ℓ
    dx = xτ * lpd'                   # x'(ζ) = x'(τ) * ℓ
    return x, dx
end
# TODO: store ℓ (not ℓᵀ) in lpv? (Why not?)

@inline function reference_diff(cache, j)
    dv = cache.dv
    lpd = cache.lagrange_polynomial_driv
    vτ = get_samples(cache, cache.reference, j)
    dv = vτ * lpd'  # x'(ζ) = x'(τ) * ℓ
    return dv
end

"""
    turns(num_mesh::Integer, degree::Integer, j::Integer)
    turns(cache, j::Integer)

Return an iterable of scaled times ([turn]) for `j`-th mesh.

[turn]: https://en.wikipedia.org/wiki/Turn_(geometry)
"""
turns(cache, j) = turns(num_mesh(cache), degree(cache), j)

function turns(num_mesh::Integer, degree::Integer, j::Integer)
    interval = 1 / num_mesh
    dt = interval / degree
    return range(interval * (j - 1), step=dt, length=degree)
end

# ------------------------------------------------------------------- residual!

function residual!(H, u, cache::LimitCycleCache)
    prob = cache.prob
    q = set(prob.de_prob.p, prob.param_axis, u[u_idx_param(cache)])
    return residual_lc!(H, u, q, cache)
end

function residual_lc!(H, u, q, cache::LimitCycleCache)
    prob = cache.prob

    l = u[u_idx_period(cache)]
    weight = cache.gauss_quadrature_weight

    phase_condition = zero(eltype(H))
    for j in 1:num_mesh(cache)
        x, dx = collocation(cache, u, j)
        dv = reference_diff(cache, j)
        @inbounds for (i, t) in enumerate(turns(cache, j))
            r = diff_range(cache, j, i)
            f = @view H[r]
            prob_time = t * l + prob.time_offset
            prob.de_prob.f(f, x[:, i], q, prob_time)

            @. H[r] = dx[:, i] - l * f
            phase_condition += weight[i] * (x[:, i] ⋅ dv[:, i])
        end
    end
    H[H_idx_phase_condition(cache)] = phase_condition

    return H
end
# Note: phase_condition above is integrated at collocation points
# (Gauss nodes) ζᵢ, not the Lagrange nodes τₖ, as mentioned in
# Kuznetsov (1998).

# ---------------------------------------------------------- residual_jacobian!

function residual_jacobian!(H, J, u, cache::LimitCycleCache)
    # TODO: implement residual_jacobian! manually
    ForwardDiff.jacobian!(
        J,
        (y, x) -> residual!(y, x, cache),
        H,  # y
        u,  # x
    )
    return (H, J)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApimachineryPkgApisMetaV1APIGroupList <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    groups::Any # spec type: Union{ Nothing, Vector{IoK8sApimachineryPkgApisMetaV1APIGroup} } # spec name: groups
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind

    function IoK8sApimachineryPkgApisMetaV1APIGroupList(;apiVersion=nothing, groups=nothing, kind=nothing)
        o = new()
        validate_property(IoK8sApimachineryPkgApisMetaV1APIGroupList, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApimachineryPkgApisMetaV1APIGroupList, Symbol("groups"), groups)
        setfield!(o, Symbol("groups"), groups)
        validate_property(IoK8sApimachineryPkgApisMetaV1APIGroupList, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        o
    end
end # type IoK8sApimachineryPkgApisMetaV1APIGroupList

const _property_map_IoK8sApimachineryPkgApisMetaV1APIGroupList = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("groups")=>Symbol("groups"), Symbol("kind")=>Symbol("kind"))
const _property_types_IoK8sApimachineryPkgApisMetaV1APIGroupList = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("groups")=>"Vector{IoK8sApimachineryPkgApisMetaV1APIGroup}", Symbol("kind")=>"String")
Base.propertynames(::Type{ IoK8sApimachineryPkgApisMetaV1APIGroupList }) = collect(keys(_property_map_IoK8sApimachineryPkgApisMetaV1APIGroupList))
Swagger.property_type(::Type{ IoK8sApimachineryPkgApisMetaV1APIGroupList }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApimachineryPkgApisMetaV1APIGroupList[name]))}
Swagger.field_name(::Type{ IoK8sApimachineryPkgApisMetaV1APIGroupList }, property_name::Symbol) =  _property_map_IoK8sApimachineryPkgApisMetaV1APIGroupList[property_name]

function check_required(o::IoK8sApimachineryPkgApisMetaV1APIGroupList)
    (getproperty(o, Symbol("groups")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApimachineryPkgApisMetaV1APIGroupList }, name::Symbol, val)
end


text:
-- prints 36
x = 1 + 3
y = x + 5
z = x * y
print(z)

max_stars_repo_name: JCGrant/JLang

max_issues_repo_name: JCGrant/JLang

max_forks_repo_name: JCGrant/JLang


text:
# This should probably be its own package

_cufunc(f,x) = f
_cufunc(f,x,xs...) = _cufunc(f, xs...)

using MacroTools

isbcastop(x) = isexpr(x, :call) && x.args[1] in :[.*,./,.+,.-].args

broadcast_inputs(ex) =
  ex isa Symbol ? [ex] :
  @capture(ex, f_.(args__)) ? vcat(broadcast_inputs.(args)...) :
  isbcastop(ex) ? vcat(broadcast_inputs.(ex.args[2:end])...) :
  []

cudata(x) = x

macro fix(ex)
  fs = []
  ex = MacroTools.prewalk(ex) do ex
    @capture(ex, f_.(args__)) || return ex
    # May not work in cases like `x .+ log(1.0)` but w/e
    xs = broadcast_inputs(ex)
    f_ = gensym()
    push!(fs, :($f_ = $_cufunc($f, $cudata.(($(xs...),))...)))
    :($f_.($(args...)))
  end
  :($(fs...); $ex) |> esc
end


text:
using LightGraphs
using MetaGraphs

"""
    disconnect_param!(obj::AbstractCompositeComponentDef, comp_def::AbstractComponentDef, param_name::Symbol)

Remove any parameter connections for a given parameter `param_name` in a given component
`comp_def` which must be a direct subcomponent of composite `obj`.
"""
function disconnect_param!(obj::AbstractCompositeComponentDef, comp_def::AbstractComponentDef, param_name::Symbol)
    # If the path isn't set yet, we look for a comp in the eventual location
    path = @or(comp_def.comp_path, ComponentPath(obj, comp_def.name))

    # @info "disconnect_param!($(obj.comp_path), $path, :$param_name)"

    if is_descendant(obj, comp_def) === nothing
        error("Cannot disconnect a component ($path) that is not within the given composite ($(obj.comp_path))")
    end

    filter!(x -> !(x.dst_comp_path == path && x.dst_par_name == param_name), obj.internal_param_conns)

    if obj isa ModelDef
        filter!(x -> !(x.comp_path == path && x.param_name == param_name), obj.external_param_conns)
    end
    dirty!(obj)
end

"""
    disconnect_param!(obj::AbstractCompositeComponentDef, comp_name::Symbol, param_name::Symbol)

Remove any parameter connections for a given parameter `param_name` in a given component
`comp_def` which must be a direct subcomponent of composite `obj`.
"""
function disconnect_param!(obj::AbstractCompositeComponentDef, comp_name::Symbol, param_name::Symbol)
    comp = compdef(obj, comp_name)
    comp === nothing && error("Did not find $comp_name in composite $(printable(obj.comp_path))")
    disconnect_param!(obj, comp, param_name)
end

# Default string, string unit check function
verify_units(unit1::AbstractString, unit2::AbstractString) = (unit1 == unit2)

function _check_labels(obj::AbstractCompositeComponentDef,
                       comp_def::AbstractComponentDef, param_name::Symbol, ext_param::ArrayModelParameter)
    param_def = parameter(comp_def, param_name)

    t1 = eltype(ext_param.values)
    t2 = eltype(param_def.datatype)
    if !(t1 <: Union{Missing, t2})
        error("Mismatched datatype of parameter connection: Component: $(comp_def.comp_id) ($t1), Parameter: $param_name ($t2)")
    end

    comp_dims  = dim_names(param_def)
    param_dims = dim_names(ext_param)

    if ! isempty(param_dims) && size(param_dims) != size(comp_dims)
        d1 = size(comp_dims)
        d2 = size(param_dims)
        error("Mismatched dimensions of parameter connection: Component: $(comp_def.comp_id) ($d1), Parameter: $param_name ($d2)")
    end

    # Don't check sizes for ConnectorComps since they won't match.
    if nameof(comp_def) in (:ConnectorCompVector, :ConnectorCompMatrix)
        return nothing
    end

    # index_values = indexvalues(obj)

    for (i, dim) in enumerate(comp_dims)
        if isa(dim, Symbol)
            param_length = size(ext_param.values)[i]
            comp_length = dim_count(obj, dim)
            if param_length != comp_length
                error("Mismatched data size for a parameter connection: dimension :$dim in $(comp_def.comp_id) has $comp_length elements; external parameter :$param_name has $param_length elements.")
            end
        end
    end
end

"""
    connect_param!(obj::AbstractCompositeComponentDef, comp_name::Symbol, param_name::Symbol, ext_param_name::Symbol;
                   check_labels::Bool=true)

Connect a parameter `param_name` in the component `comp_name` of composite `obj` to
the external parameter `ext_param_name`.
"""
function connect_param!(obj::AbstractCompositeComponentDef, comp_name::Symbol,
                        param_name::Symbol, ext_param_name::Symbol;
                        check_labels::Bool=true)
    comp_def = compdef(obj, comp_name)
    connect_param!(obj, comp_def, param_name, ext_param_name, check_labels=check_labels)
end

function connect_param!(obj::AbstractCompositeComponentDef, comp_def::AbstractComponentDef,
                        param_name::Symbol, ext_param_name::Symbol; check_labels::Bool=true)
    ext_param = external_param(obj, ext_param_name)

    if ext_param isa ArrayModelParameter && check_labels
        _check_labels(obj, comp_def, param_name, ext_param)
    end

    disconnect_param!(obj, comp_def, param_name)    # calls dirty!()

    comp_path = @or(comp_def.comp_path, ComponentPath(obj.comp_path, comp_def.name))
    conn = ExternalParameterConnection(comp_path, param_name, ext_param_name)
    add_external_param_conn!(obj, conn)

    return nothing
end

"""
    _connect_param!(obj::AbstractCompositeComponentDef,
        dst_comp_path::ComponentPath, dst_par_name::Symbol,
        src_comp_path::ComponentPath, src_var_name::Symbol,
        backup::Union{Nothing, Array}=nothing;
        ignoreunits::Bool=false, backup_offset::Int=0)

Bind the parameter `dst_par_name` of one component `dst_comp_path` of composite `obj` to a
variable `src_var_name` in another component `src_comp_path` of the same model using
`backup` to provide default values and the `ignoreunits` flag to indicate the need to
check match units between the two.  The `backup_offset` argument, which is only valid
when `backup` data has been set, indicates that the backup data should be used for
a specified number of timesteps after the source component begins. ie. the value would be
`1` if the destination componentm parameter should only use the source component
data for the second timestep and beyond.
"""
function _connect_param!(obj::AbstractCompositeComponentDef,
                        dst_comp_path::ComponentPath, dst_par_name::Symbol,
                        src_comp_path::ComponentPath, src_var_name::Symbol,
                        backup::Union{Nothing, Array}=nothing;
                        ignoreunits::Bool=false, backup_offset::Union{Nothing, Int}=nothing)

    dst_comp_def = compdef(obj, dst_comp_path)
    src_comp_def = compdef(obj, src_comp_path)

    # remove any existing connections for this dst parameter
    disconnect_param!(obj, dst_comp_def, dst_par_name)  # calls dirty!()

    # @info "dst_comp_def: $dst_comp_def"
    # @info "src_comp_def: $src_comp_def"

    if backup !== nothing
        # If value is a NamedArray, we can check if the labels match
        if isa(backup, NamedArray)
            dims = dimnames(backup)
            check_parameter_dimensions(obj, backup, dims, dst_par_name)
        else
            dims = nothing
        end

        # Check that the backup data is the right size
        if size(backup) != datum_size(obj, dst_comp_def, dst_par_name)
            error("Cannot connect parameter; the provided backup data is the wrong size. ",
                  "Expected size $(datum_size(obj, dst_comp_def, dst_par_name)) but got $(size(backup)).")
        end

        # convert number type and, if it's a NamedArray, convert to Array
        backup = convert(Array{Union{Missing, number_type(obj)}}, backup)

        dst_param = parameter(dst_comp_def, dst_par_name)
        dst_dims  = dim_names(dst_param)
        dim_count = length(dst_dims)

        ti = get_time_index_position(dst_param)

        if ti === nothing # not time dimension
            values = backup
        else # handle time dimension

            # get first and last of the ModelDef, NOT the ComponentDef
            first = first_period(obj)
            last = last_period(obj)

            T = eltype(backup)

            if isuniform(obj)
                stepsize = step_size(obj)
                values = TimestepArray{FixedTimestep{first, stepsize, last}, T, dim_count, ti}(backup)
            else
                times = time_labels(obj)
                values = TimestepArray{VariableTimestep{(times...,)}, T, dim_count, ti}(backup)
            end

        end

        set_external_array_param!(obj, dst_par_name, values, dst_dims)
        backup_param_name = dst_par_name

    else
        # cannot use backup_offset keyword argument if there is no backup
        if backup_offset !== nothing
            error("Cannot set `backup_offset` keyword argument if `backup` data is not explicitly provided")
        end

        # If backup not provided, make sure the source component covers the span of the destination component
        src_first, src_last = first_and_last(src_comp_def)
        dst_first, dst_last = first_and_last(dst_comp_def)
        if (dst_first !== nothing && src_first !== nothing && dst_first < src_first) ||
            (dst_last  !== nothing && src_last  !== nothing && dst_last  > src_last)
            src_first = printable(src_first)
            src_last  = printable(src_last)
            dst_first = printable(dst_first)
            dst_last  = printable(dst_last)
            error("""Cannot connect parameter: $src_comp_path runs only from $src_first to $src_last,
whereas $dst_comp_path runs from $dst_first to $dst_last. Backup data must be provided for missing years.
Try calling:
    `connect_param!(m, comp_name, par_name, comp_name, var_name, backup_data)`""")
        end

        backup_param_name = nothing
    end

    # Check the units, if provided
    if ! ignoreunits && ! verify_units(variable_unit(src_comp_def, src_var_name),
                                       parameter_unit(dst_comp_def, dst_par_name))
        error("Units of $src_comp_path:$src_var_name do not match $dst_comp_path:$dst_par_name.")
    end

    conn = InternalParameterConnection(src_comp_path, src_var_name, dst_comp_path, dst_par_name,
                                       ignoreunits, backup_param_name, backup_offset=backup_offset)
    add_internal_param_conn!(obj, conn)

    return nothing
end

function connect_param!(obj::AbstractCompositeComponentDef,
                        dst_comp_name::Symbol, dst_par_name::Symbol,
                        src_comp_name::Symbol, src_var_name::Symbol,
                        backup::Union{Nothing, Array}=nothing; ignoreunits::Bool=false,
                        backup_offset::Union{Nothing, Int} = nothing)
    _connect_param!(obj, ComponentPath(obj, dst_comp_name), dst_par_name,
                        ComponentPath(obj, src_comp_name), src_var_name,
                        backup; ignoreunits=ignoreunits, backup_offset=backup_offset)
end

"""
    connect_param!(obj::AbstractCompositeComponentDef,
        dst::Pair{Symbol, Symbol}, src::Pair{Symbol, Symbol},
        backup::Union{Nothing, Array}=nothing;
        ignoreunits::Bool=false, backup_offset::Union{Nothing, Int} = nothing)

Bind the parameter `dst[2]` of one component `dst[1]` of composite `obj`
to a variable `src[2]` in another component `src[1]` of the same composite
using `backup` to provide default values and the `ignoreunits` flag to indicate the need
to check match units between the two.  The `backup_offset` argument, which is only valid
when `backup` data has been set, indicates that the backup data should be used for
a specified number of timesteps after the source component begins. ie. the value would be
`1` if the destination componentm parameter should only use the source component
data for the second timestep and beyond.
"""
function connect_param!(obj::AbstractCompositeComponentDef,
                        dst::Pair{Symbol, Symbol}, src::Pair{Symbol, Symbol},
                        backup::Union{Nothing, Array}=nothing; ignoreunits::Bool=false,
                        backup_offset::Union{Nothing, Int} = nothing)
    connect_param!(obj, dst[1], dst[2], src[1], src[2], backup; ignoreunits=ignoreunits, backup_offset=backup_offset)
end

"""
    split_datum_path(obj::AbstractCompositeComponentDef, s::AbstractString)

Split a string of the form "/path/to/component:datum_name" into the component path,
`ComponentPath(:path, :to, :component)` and name `:datum_name`.
"""
function split_datum_path(obj::AbstractCompositeComponentDef, s::AbstractString)
    elts = split(s, ":")
    length(elts) != 2 && error("Cannot split datum path '$s' into ComponentPath and datum name")
    return (ComponentPath(obj, elts[1]), Symbol(elts[2]))
end

"""
    connection_refs(obj::AbstractCompositeComponentDef)

Return a vector of UnnamedReference's to parameters from subcomponents that are either found in
internal connections or that have been already imported in a parameter definition.
"""
function connection_refs(obj::AbstractCompositeComponentDef)
    refs = UnnamedReference[]

    for conn in obj.internal_param_conns
        push!(refs, UnnamedReference(conn.dst_comp_path.names[end], conn.dst_par_name))
    end

    for item in values(obj.namespace)
        if item isa CompositeParameterDef
            for ref in item.refs
                push!(refs, ref)
            end
        end
    end

    return refs
end

"""
    connection_refs(obj::ModelDef)

Return a vector of UnnamedReference's to parameters from subcomponents that are either found in
internal connections or that have been already connected to external parameter values.
"""
function connection_refs(obj::ModelDef)
    refs = UnnamedReference[]

    for conn in obj.internal_param_conns
        push!(refs, UnnamedReference(conn.dst_comp_path.names[end], conn.dst_par_name))
    end

    for conn in obj.external_param_conns
        push!(refs, UnnamedReference(conn.comp_path.names[end], conn.param_name))
    end

    return refs
end

"""
    unconnected_params(obj::AbstractCompositeComponentDef)

Return a list of UnnamedReference's to parameters that have not been connected
to a value.
"""
function unconnected_params(obj::AbstractCompositeComponentDef)
    return setdiff(subcomp_params(obj), connection_refs(obj))
end

"""
    set_leftover_params!(m::Model, parameters::Dict)

Set all of the parameters in model `m` that don't have a value and are not connected
to some other component to a value from a dictionary `parameters`. This method assumes
the dictionary keys are strings that match the names of unset parameters in the model.
"""
function set_leftover_params!(md::ModelDef, parameters::Dict{T, Any}) where T
    for param_ref in unconnected_params(md)
        param_name = param_ref.datum_name
        comp_name = param_ref.comp_name
        comp_def = find_comp(md, comp_name)
        param_def = comp_def[param_name]

        # Only set the unconnected parameter if it doesn't have a default
        if param_def.default === nothing
            # check whether we need to create the external parameter
            if external_param(md, param_name, missing_ok=true) === nothing
                if haskey(parameters, string(param_name))
                    value = parameters[string(param_name)]
                    param_dims = parameter_dimensions(md, comp_name, param_name)

                    set_external_param!(md, param_name, value; param_dims = param_dims)
                else
                    error("Cannot set parameter :$param_name, not found in provided dictionary and no default value detected.")
                end
            end
            connect_param!(md, comp_name, param_name, param_name)
        end
    end
    nothing
end

# Find internal param conns to a given destination component
function internal_param_conns(obj::AbstractCompositeComponentDef, dst_comp_path::ComponentPath)
    return filter(x->x.dst_comp_path == dst_comp_path, internal_param_conns(obj))
end

function internal_param_conns(obj::AbstractCompositeComponentDef, comp_name::Symbol)
    return internal_param_conns(obj, ComponentPath(obj.comp_path, comp_name))
end

function add_internal_param_conn!(obj::AbstractCompositeComponentDef, conn::InternalParameterConnection)
    push!(obj.internal_param_conns, conn)
    dirty!(obj)
end

#
# These should all take ModelDef instead of AbstractCompositeComponentDef as 1st argument
#

# Find external param conns for a given comp
function external_param_conns(obj::ModelDef, comp_path::ComponentPath)
    return filter(x -> x.comp_path == comp_path, external_param_conns(obj))
end

function external_param_conns(obj::ModelDef, comp_name::Symbol)
    return external_param_conns(obj, ComponentPath(obj.comp_path, comp_name))
end

function external_param(obj::ModelDef, name::Symbol; missing_ok=false)
    haskey(obj.external_params, name) && return obj.external_params[name]

    missing_ok && return nothing

    error("$name not found in external parameter list")
end

function add_external_param_conn!(obj::ModelDef, conn::ExternalParameterConnection)
    push!(obj.external_param_conns, conn)
    dirty!(obj)
end

function set_external_param!(obj::ModelDef, name::Symbol, value::ModelParameter)
    # if haskey(obj.external_params, name)
    #     @warn "Redefining external param :$name in $(obj.comp_path) from $(obj.external_params[name]) to $value"
    # end
    obj.external_params[name] = value
    dirty!(obj)
    return value
end

function set_external_param!(obj::ModelDef, name::Symbol, value::Number;
                             param_dims::Union{Nothing,Array{Symbol}} = nothing)
    set_external_scalar_param!(obj, name, value)
end

function set_external_param!(obj::ModelDef, name::Symbol,
                             value::Union{AbstractArray, AbstractRange, Tuple};
                             param_dims::Union{Nothing,Array{Symbol}} = nothing)
    ti = get_time_index_position(param_dims)
    if ti != nothing
        value = convert(Array{number_type(obj)}, value)
        num_dims = length(param_dims)
        values = get_timestep_array(obj, eltype(value), num_dims, ti, value)
    else
        values = value
    end

    set_external_array_param!(obj, name, values, param_dims)
end

"""
    set_external_array_param!(obj::ModelDef,
                              name::Symbol, value::TimestepVector, dims)

Add a one dimensional time-indexed array parameter indicated by `name` and
`value` to the composite `obj`.  In this case `dims` must be `[:time]`.
"""
function set_external_array_param!(obj::ModelDef,
                                   name::Symbol, value::TimestepVector, dims)
    param = ArrayModelParameter(value, [:time])  # must be :time
    set_external_param!(obj, name, param)
end

"""
    set_external_array_param!(obj::ModelDef,
                              name::Symbol, value::TimestepMatrix, dims)

Add a multi-dimensional time-indexed array parameter `name` with value
`value` to the composite `obj`.  In this case `dims` must be `[:time]`.
"""
function set_external_array_param!(obj::ModelDef,
                                   name::Symbol, value::TimestepArray, dims)
    param = ArrayModelParameter(value, dims === nothing ? Vector{Symbol}() : dims)
    set_external_param!(obj, name, param)
end

"""
    set_external_array_param!(obj::ModelDef,
                              name::Symbol, value::AbstractArray, dims)

Add an array type parameter `name` with value `value` and `dims` dimensions to the composite `obj`.
"""
function set_external_array_param!(obj::ModelDef,
                                   name::Symbol, value::AbstractArray, dims)
    param = ArrayModelParameter(value, dims === nothing ? Vector{Symbol}() : dims)
    set_external_param!(obj, name, param)
end

"""
    set_external_scalar_param!(obj::ModelDef, name::Symbol, value::Any)

Add a scalar type parameter `name` with the value `value` to the composite `obj`.
"""
function set_external_scalar_param!(obj::ModelDef, name::Symbol, value::Any)
    param = ScalarModelParameter(value)
    set_external_param!(obj, name, param)
end

"""
    update_param!(obj::AbstractCompositeComponentDef, name::Symbol, value; update_timesteps = nothing)

Update the `value` of an external model parameter in composite `obj`, referenced
by `name`. The update_timesteps keyword argument is deprecated, we keep it here
just to provide warnings.
"""
function update_param!(obj::AbstractCompositeComponentDef, name::Symbol, value; update_timesteps = nothing)
    !isnothing(update_timesteps) ? @warn("Use of the `update_timesteps` keyword argument is no longer supported or needed, time labels will be adjusted automatically if necessary.") : nothing
    _update_param!(obj::AbstractCompositeComponentDef, name, value)
end

function update_param!(mi::ModelInstance, name::Symbol, value)
    param = mi.md.external_params[name]

    if param isa ScalarModelParameter
        param.value = value
    elseif param.values isa TimestepArray
        copyto!(param.values.data, value)
    else
        copyto!(param.values, value)
    end

    return nothing
end

function _update_param!(obj::AbstractCompositeComponentDef,
                        name::Symbol, value)
    param = external_param(obj, name, missing_ok=true)
    if param === nothing
        error("Cannot update parameter; $name not found in composite's external parameters.")
    end

    if param isa ScalarModelParameter
        _update_scalar_param!(param, name, value)
    else
        _update_array_param!(obj, name, value)
    end

    dirty!(obj)
end

function _update_scalar_param!(param::ScalarModelParameter, name, value)
    if ! (value isa typeof(param.value))
        try
            value = convert(typeof(param.value), value)
        catch e
            error("Cannot update parameter $name; expected type $(typeof(param.value)) but got $(typeof(value)).")
        end
    end
    param.value = value
    nothing
end

function _update_array_param!(obj::AbstractCompositeComponentDef, name, value)

    # Get original parameter
    param = external_param(obj, name)

    # Check type of provided parameter
    if !(typeof(value) <: AbstractArray)
        error("Cannot update array parameter $name with a value of type $(typeof(value)).")

    elseif !(eltype(value) <: eltype(param.values))
        try
            value = convert(Array{eltype(param.values)}, value)
        catch e
            error("Cannot update parameter $name; expected array of type $(eltype(param.values)) but got $(eltype(value)).")
        end
    end

    # Check if the parameter dimensions match the model dimensions.  Note that we
    # previously checked if parameter dimensions matched the dimensions of the
    # parameter they were to replace, but given dimensions of a model can be changed,
    # we now choose to enforce that the new dimensions match the current model state,
    # whatever that is.

    expected_size = ([length(dim_keys(obj, d)) for d in dim_names(param)]...,)
    size(value) != expected_size ? error("Cannot update parameter $name; expected array of size $expected_size but got array of size $(size(value)).") : nothing

    # check if updating timestep labels is necessary
    if param.values isa TimestepArray
        time_label_change = time_labels(param.values) != dim_keys(obj, :time)
        N = ndims(value)
        if time_label_change
            T = eltype(value)
            ti = get_time_index_position(param)
            new_timestep_array = get_timestep_array(obj, T, N, ti, value)
            set_external_param!(obj, name, ArrayModelParameter(new_timestep_array, dim_names(param)))
        else
            copyto!(param.values.data, value)
        end
    else
        copyto!(param.values, value)
    end

    dirty!(obj)
    nothing
end

"""
    update_params!(obj::AbstractCompositeComponentDef, parameters::Dict{T, Any}; update_timesteps = nothing) where T

For each (k, v) in the provided `parameters` dictionary, `update_param!`
is called to update the external parameter by name k to value v. Each key k must be a symbol or convert to a
symbol matching the name of an external parameter that already exists in the
component definition.
"""
function update_params!(obj::AbstractCompositeComponentDef, parameters::Dict; update_timesteps = nothing)
    !isnothing(update_timesteps) ? @warn("Use of the `update_timesteps` keyword argument is no longer supported or needed, time labels will be adjusted automatically if necessary.") : nothing
    parameters = Dict(Symbol(k) => v for (k, v) in parameters)
    for (param_name, value) in parameters
        _update_param!(obj, param_name, value)
    end
    nothing
end

function add_connector_comps!(obj::AbstractCompositeComponentDef)
    conns = internal_param_conns(obj)
    i = 1 # counter to track the number of connector comps added

    for comp_def in compdefs(obj)
        comp_name = nameof(comp_def)
        comp_path = comp_def.comp_path

        # first need to see if we need to add any connector components for this component
        internal_conns  = filter(x -> x.dst_comp_path == comp_path, conns)
        need_conn_comps = filter(x -> x.backup !== nothing, internal_conns)

        # isempty(need_conn_comps) || @info "Need connectors comps: $need_conn_comps"

        for conn in need_conn_comps
            add_backup!(obj, conn.backup)

            num_dims = length(size(external_param(obj, conn.backup).values))

            if ! (num_dims in (1, 2))
                error("Connector components for parameters with > 2 dimensions are not implemented.")
            end

            # Fetch the definition of the appropriate connector commponent
            conn_comp_def = (num_dims == 1 ? Mimi.ConnectorCompVector : Mimi.ConnectorCompMatrix)
            conn_comp_name = connector_comp_name(i) # generate a new name
            i += 1 # increment connector comp counter

            # Add the connector component before the user-defined component that
            # required it, and for now let the first and last of the component
            # be free and thus be set to the same as the model
            conn_comp = add_comp!(obj, conn_comp_def, conn_comp_name, before=comp_name)
            conn_path = conn_comp.comp_path

            # add a connection between src_component and the ConnectorComp
            add_internal_param_conn!(obj, InternalParameterConnection(conn.src_comp_path, conn.src_var_name,
                                                                      conn_path, :input1,
                                                                      conn.ignoreunits))

            # add a connection between ConnectorComp and dst_component
            add_internal_param_conn!(obj, InternalParameterConnection(conn_path, :output,
                                                                      conn.dst_comp_path, conn.dst_par_name,
                                                                      conn.ignoreunits))

            # add a connection between ConnectorComp and the external backup data
            add_external_param_conn!(obj, ExternalParameterConnection(conn_path, :input2, conn.backup))

            # set the first and last parameters for WITHIN the component which
            # decide when backup is used and when connection is used
            src_comp_def = compdef(obj, conn.src_comp_path)

            param_last = last_period(obj, src_comp_def)
            param_first = first_period(obj, src_comp_def)
            conn.backup_offset !== nothing ? param_first = param_first + conn.backup_offset : nothing

            set_param!(obj, conn_comp_name, :first, Symbol(conn_comp_name, "_", :first), param_first)
            set_param!(obj, conn_comp_name, :last, Symbol(conn_comp_name, "_", :last), param_last)
        end
    end

    return nothing
end

"""
    _pad_parameters!(obj::ModelDef)

Take each external parameter of the Model Definition `obj` and `update_param!`
with new data values that are altered to match a new time dimension by (1) trimming
the values down if the time dimension has been shortened and (2) padding with missings
as necessary.
"""
function _pad_parameters!(obj::ModelDef)

    model_times = time_labels(obj)

    for (name, param) in obj.external_params
        if (param isa ArrayModelParameter) && (:time in param.dim_names)

           param_times = _get_param_times(param)
           padded_data = _get_padded_data(param, param_times, model_times)
           update_param!(obj, name, padded_data)

        end
    end
end

"""
    _get_padded_data(param::ArrayModelParameter, param_times::Vector, model_times::Vector)

Obtain the new data values for the Array Model Paramter `param` with current
time labels `param_times` such that they are altered to match a new time dimension
with keys `model_times` by (1) trimming the values down if the time dimension has
been shortened and (2) padding with missings as necessary.
"""
function _get_padded_data(param::ArrayModelParameter, param_times::Vector, model_times::Vector)

    data = param.values.data
    ti = get_time_index_position(param)

    # first handle the back end
    model_last = last(model_times)
    param_last = last(param_times)

    if model_last < param_last # trim down the data

        trim_idx = findfirst(isequal(last(model_times)), param_times)
        idxs = repeat(Any[:], ndims(data))
        idxs[ti] = 1:trim_idx
        data = data[idxs...]

    elseif model_last > param_last # pad the data

        pad_length = length(model_times[findfirst(isequal(param_last), model_times)+1:end])
        dims = [size(data)...]
        dims[ti] = pad_length
        end_padding_rows = Array{Union{Missing, Number}}(missing, dims...)
        data = vcat(data, end_padding_rows)

    end

    # now handle the front end
    model_first = first(model_times)
    param_first = first(param_times)

    # note we do not allow for any trimming off the front end
    if model_first < param_first

        pad_length = length(model_times[1:findfirst(isequal(param_first), model_times)-1])
        dims = [size(data)...]
        dims[ti] = pad_length
        begin_padding_rows = Array{Union{Missing, Number}}(missing, dims...)
        data = vcat(begin_padding_rows, data)

    end

    return data
end

"""
    _get_param_times(param::ArrayModelParameter{TimestepArray{FixedTimestep{FIRST, STEP, LAST}, T, N, ti, S}})

Return the time labels that parameterize the `TimestepValue` which in turn parameterizes
the ArrayModelParameter `param`.
"""
function _get_param_times(param::ArrayModelParameter{TimestepArray{FixedTimestep{FIRST, STEP, LAST}, T, N, ti, S}}) where {FIRST, STEP, LAST, T, N, ti, S}
    return collect(FIRST:STEP:LAST)
end

"""
    _get_param_times(param::ArrayModelParameter{TimestepArray{VariableTimestep{TIMES}, T, N, ti, S}})

Return the time labels that parameterize the `TimestepValue` which in turn parameterizes
the ArrayModelParameter `param`.
"""
function _get_param_times(param::ArrayModelParameter{TimestepArray{VariableTimestep{TIMES}, T, N, ti, S}}) where {TIMES, T, N, ti, S}
    return [TIMES...]
end


text:
using Pkg, VisualStringDistances
using VisualStringDistances: KL, sinkhorn_divergence!
using StringDistances
using DataFrames
using Transducers

const DL = DamerauLevenshtein()

function get_all_package_names(registry_dir::AbstractString)
    packages = [x["name"] for x in values(Pkg.TOML.parsefile(joinpath(registry_dir, "Registry.toml"))["packages"])]
    sort!(packages)
    unique!(packages)
    return packages
end

names = get_all_package_names(expanduser("~/.julia/registries/General"))

filter!(x -> !endswith(x, "_jll"), names)
@info "Loaded list of non-JLL package names ($(length(names)) names)"

normalized_dl_cutoff = .2
dl_cutoff = 1
@info "Computing list of pairs of package names within $(dl_cutoff) in DL distance or $(normalized_dl_cutoff) in normalized DL distance"
@time df = DataFrame(tcollect( (name1=names[i],name2=names[j]) for i = 1:length(names) for j = 1:(i-1) if (normalize(DL)(names[i], names[j]) <= normalized_dl_cutoff) || DL(names[i], names[j]) <= dl_cutoff))

df.longest_length = max.(length.(df.name1), length.(df.name2))

function compute_sd(df, penalty, ϵ)
    tcollect( sinkhorn_divergence!(penalty, word_measure(df.name1[i]), word_measure(df.name2[i]), ϵ) for i = 1:size(df,1))
end

@info "Computing DL distances for pairs..."
col_dict = Dict{Any, String}()
col = "DL unnormalized"
col_dict[(name = :DL, normalization = :unnormalized, params=tuple())] = col
df[!, col] = DL.(df.name1, df.name2)

i = 0
for ϵ in (0.1, 0.5, 1.0), ρ in (1.0, 5.0, 10.0)
    global i += 1
    @info "($(i)/9) Computing sinkhorn divergences with ϵ=$ϵ, ρ=$ρ..."
    col = "SD ϵ=$(ϵ) ρ=$(ρ) unnormalized"
    col_dict[(name = :SD, normalization = :unnormalized, params=(ϵ=ϵ, ρ=ρ))] = col
    @time df[!, col] = compute_sd(df, KL(ρ), ϵ)
end

@info "Computing normalized distances..."
@time begin
    for ϵ in (0.1, 0.5, 1.0), ρ in (1.0, 5.0, 10.0)
        col = "SD ϵ=$(ϵ) ρ=$(ρ) normalized"
        col_dict[(name = :SD, normalization = :normalized, params=(ϵ=ϵ, ρ=ρ))] = col
        df[!,  col] = df[!, col_dict[(name = :SD, normalization = :unnormalized, params=(ϵ=ϵ, ρ=ρ))]] ./ df[!, :longest_length]

        col = "SD ϵ=$(ϵ) ρ=$(ρ) sqrt normalized"
        col_dict[(name = :SD, normalization = :sqrt_normalized, params=(ϵ=ϵ, ρ=ρ))] = col
        df[!,  col] = df[!, col_dict[(name = :SD, normalization = :unnormalized, params=(ϵ=ϵ, ρ=ρ))]] ./ sqrt.(df[!, :longest_length])
    end
end

col = "DL sqrt normalized"
col_dict[(name = :DL, normalization = :sqrt_normalized, params=tuple())] = col
df[!, col] =  df[!, col_dict[(name = :DL, normalization = :unnormalized, params=tuple())]] ./ sqrt.(df[!, :longest_length])

col = "DL normalized"
col_dict[(name = :DL, normalization = :normalized, params=tuple())] = col
df[!, col] =  df[!, col_dict[(name = :DL, normalization = :unnormalized, params=tuple())]] ./ df[!, :longest_length]

using Serialization
@info "Serializing..."
serialize(joinpath(@__DIR__, "names_df.jls"), Dict{Any, Any}(:df => df, :col_dict => col_dict))

@info "Done!"


text:
push!(LOAD_PATH,"../src/")
using Documenter, Stars

## References:
# https://github.com/Alexander-Barth/NCDatasets.jl/blob/master/docs/make.jl

CI = get(ENV, "CI", nothing) == "true"
makedocs(
    sitename="Stars.jl",
    format = Documenter.HTML(prettyurls = CI,),
    pages = [
        "Introduction"   => "Index.md",
        "Variables"      => "Variables.md",
        "IO" => "IO.md",
    ],
)

deploydocs(
    repo = "github.com/jl-spatial/Stars.jl.git",
)


text:
import Base: length, getindex, setindex!, size, copy, deepcopy, empty!,
             convert, transpose, ctranspose, names

"""
MIToS MSAs are subtypes of `AbstractMatrix{Residue}`,
because the most basic implementation of a MIToS MSA is a `Matrix` of `Residue`s.
"""
abstract AbstractMultipleSequenceAlignment <: AbstractMatrix{Residue}

"""
MIToS sequences are subtypes of `AbstractVector{Residue}`.
"""
abstract AbstractAlignedSequence <: AbstractVector{Residue}

# Multiple Sequence Alignment
# ===========================

"""
This MSA type include the `Matrix` of `Residue`s and the sequence names.
To allow fast indexing of MSAs using **sequence identifiers**,
they are saved as an `IndexedArray`.
"""
type MultipleSequenceAlignment <: AbstractMultipleSequenceAlignment
    id::IndexedArray{ASCIIString}
    msa::Matrix{Residue}
end

"""
...
"""
type AnnotatedMultipleSequenceAlignment <: AbstractMultipleSequenceAlignment
    id::IndexedArray{ASCIIString}
    msa::Matrix{Residue}
    annotations::Annotations
end

convert(::Type{MultipleSequenceAlignment}, msa::AnnotatedMultipleSequenceAlignment) = MultipleSequenceAlignment(msa.id, msa.msa)

# Aligned Sequence
# ================

type AlignedSequence <: AbstractAlignedSequence
    id::ASCIIString
    index::Int
    sequence::Vector{Residue}
end

type AnnotatedAlignedSequence <: AbstractAlignedSequence
    id::ASCIIString
    index::Int
    sequence::Vector{Residue}
    annotations::Annotations
end

convert(::Type{AlignedSequence}, seq::AnnotatedAlignedSequence) = AlignedSequence(seq.id, seq.index, seq.sequence)

# AbstractArray Interface
# -----------------------

for meth in (:size, :length)
    @eval $(meth)(msa::AbstractMultipleSequenceAlignment) = $(meth)(msa.msa)
    @eval $(meth)(seq::AbstractAlignedSequence) = $(meth)(seq.sequence)
end

for T in (:(AlignedSequence), :(AnnotatedAlignedSequence),
          :(MultipleSequenceAlignment), :(AnnotatedMultipleSequenceAlignment))
    @eval Base.linearindexing(::Type{$(T)}) = Base.LinearFast()
end

getindex(msa::AbstractMultipleSequenceAlignment, i::Int) = getindex(msa.msa, i)
getindex(seq::AbstractAlignedSequence, i::Int) = getindex(seq.sequence, i)
setindex!(msa::AbstractMultipleSequenceAlignment, value::Residue, i::Int) =  setindex!(msa.msa, value, i)
setindex!(seq::AbstractAlignedSequence, value::Residue, i::Int) = setindex!(seq.sequence, value, i)

# Transpose
# ---------

transpose(msa::AbstractMultipleSequenceAlignment)  = transpose(msa.msa)
ctranspose(msa::AbstractMultipleSequenceAlignment) = transpose(msa.msa) # transpose is ~ 0.00022 seconds faster than ctranspose for PF00085

# Selection without Mappings
# --------------------------

"""
Allows you to access the residues in a `Matrix{Residues}`/`Vector{Residues}` without annotations.
"""
getresidues(msa::AbstractMultipleSequenceAlignment) = msa.msa
getresidues(seq::AbstractAlignedSequence) = seq.sequence

"""
Gives you the number of sequences on the `MultipleSequenceAlignment`
"""
nsequences(msa::AbstractMultipleSequenceAlignment) = size(msa.msa, 1)
nsequences(msa::Matrix{Residue}) = size(msa, 1)

"""
`ncolumns` returns the MSA columns or sequence positions number.
"""
ncolumns(msa::AbstractMultipleSequenceAlignment) = size(msa.msa, 2)
ncolumns(msa::Matrix{Residue}) = size(msa, 2)
ncolumns(seq::AbstractAlignedSequence) = length(seq.sequence)
ncolumns(seq::Vector{Residue}) = length(seq)

"""
Gives you a `Vector{Vector{Residue}}` with all the sequences of the MSA without Annotations
"""
function getresiduesequences(msa::Matrix{Residue})
    nseq = nsequences(msa)
    tmsa = msa'
    sequences = Array(Vector{Residue}, nseq)
    for i in 1:nseq
        @inbounds sequences[i] = tmsa[:,i]
    end
    sequences
end

getresiduesequences(msa::AbstractMultipleSequenceAlignment) = getresiduesequences(msa.msa)

# Select sequence
# ---------------

# Gives you the annotations of the Sequence
function getsequence(data::Annotations, id::ASCIIString)
    GS = Dict{Tuple{ASCIIString,ASCIIString},ASCIIString}()
    GR = Dict{Tuple{ASCIIString,ASCIIString},ASCIIString}()
    if length(data.sequences) > 0 || length(data.residues) > 0
        for (key, value) in data.sequences
            if key[1] == id
                GS[key] = value
            end
        end
        for (key, value) in data.residues
            if key[1] == id
                GR[key] = value
            end
        end
        GS = sizehint!(GS, length(GS))
        GR = sizehint!(GR, length(GR))
    end
    Annotations(data.file, GS, data.columns, GR)
end

@doc """
`getsequence` takes an MSA and a sequence number or identifier to return a sequence.
If the MSA is n `AnnotatedMultipleSequenceAlignment`, it returns an `AnnotatedAlignedSequence` with sequence annotations.
From a `MultipleSequenceAlignment`, It returns an `AlignedSequence` object.
If an `Annotations` object  and a sequence identifier are used, this function returns the annotations related to the sequence.
""" getsequence

getsequence(msa::AnnotatedMultipleSequenceAlignment,
            i::Int) = AnnotatedAlignedSequence(msa.id[i], i,  vec(msa.msa[i,:]),
                                               getsequence(msa.annotations, msa.id[i]))

getsequence(msa::MultipleSequenceAlignment,
            i::Int) = AlignedSequence(msa.id[i], i, vec(msa.msa[i,:]))

getsequence(msa::Matrix{Residue}, i::Int) = vec(msa[i,:])

function getsequence(msa::AnnotatedMultipleSequenceAlignment, id::ASCIIString)
    i = findfirst(msa.id, id)
    AnnotatedAlignedSequence(id, i, vec(msa.msa[i,:]),
                             getsequence(msa.annotations, id))
end

function getsequence(msa::MultipleSequenceAlignment, id::ASCIIString)
    i = findfirst(msa.id, id)
    AlignedSequence(id, i, vec(msa.msa[i,:]))
end

getindex(msa::AbstractMultipleSequenceAlignment, id::ASCIIString) = getsequence(msa, id)

# Filters
# -------

"It's similar to `filtersequences!` but for a `Matrix` of `Residue`s"
filtersequences(msa::Matrix{Residue}, mask::AbstractVector{Bool}) = msa[mask, :]

"""
`filtersequences!(msa, mask[, annotate::Bool=true])`

Allows to filter `msa` sequences using a `AbstractVector{Bool}` `mask` (removes `false`s).
`AnnotatedMultipleSequenceAlignment` annotations are updated if `annotate` is `true` (default).
"""
function filtersequences!(msa::AnnotatedMultipleSequenceAlignment,
                          mask::AbstractVector{Bool}, annotate::Bool=true)
    msa.msa = filtersequences(msa.msa, mask)
    #msa.sequencemapping = msa.sequencemapping[ mask , : ]
    filtersequences!(msa.annotations, msa.id, mask)
    msa.id = IndexedArray(msa.id[ mask ])
    annotate && annotate_modification!(msa, string("filtersequences! : ", sum(~mask), " sequences have been deleted."))
    msa
end

function filtersequences!(msa::MultipleSequenceAlignment,
                          mask::AbstractVector{Bool}, annotate::Bool=false) # annotate is useful for calling this inside other functions
    msa.msa = filtersequences(msa.msa, mask)
    msa.id = IndexedArray(msa.id[ mask ])
    msa
end

"It's similar to `filtercolumns!` but for a `Matrix` or `Vector` of `Residue`s"
filtercolumns(msa::Matrix{Residue}, mask::AbstractVector{Bool}) = msa[ : , mask ]
filtercolumns(seq::Vector{Residue}, mask::AbstractVector{Bool}) = seq[ mask ]

"""
`filtercolumns!(msa, mask[, annotate::Bool=true])`

Allows to filter MSA columns/sequence positions using a `AbstractVector{Bool}` `mask`.
`AnnotatedMultipleSequenceAlignment`s or `AnnotatedAlignedSequence`s annotations are updated.
"""
function filtercolumns!(msa::AnnotatedMultipleSequenceAlignment,
                        mask::AbstractVector{Bool}, annotate::Bool=true)
    msa.msa = filtercolumns(msa.msa, mask)
    filtercolumns!(msa.annotations, mask)
    annotate && annotate_modification!(msa, string("filtercolumns! : ",
                                                   sum(~mask), " columns have been deleted."))
    msa
end

function filtercolumns!(msa::MultipleSequenceAlignment, mask::AbstractVector{Bool}, annotate::Bool=false) # annotate is useful for calling this inside other functions
    msa.msa = filtercolumns(msa.msa, mask)
    msa
end

function filtercolumns!(seq::AnnotatedAlignedSequence, mask::AbstractVector{Bool}, annotate::Bool=true)
    seq.sequence = filtercolumns(seq.sequence, mask)
    filtercolumns!(seq.annotations, mask)
    annotate && annotate_modification!(seq, string("filtercolumns! : ",
                                                   sum(~mask), " columns have been deleted."))
    seq
end

function filtercolumns!(seq::AlignedSequence, mask::AbstractVector{Bool}, annotate::Bool=false)
    seq.sequence = filtercolumns(seq.sequence, mask)
    seq
end

filtercolumns(seq::AbstractAlignedSequence, args...) = filtercolumns!(deepcopy(seq), args...)
filtercolumns(msa::AbstractMultipleSequenceAlignment, args...) = filtercolumns!(deepcopy(msa), args...)

# Copy, deepcopy, empty!
# ----------------------

for fun in (:copy, :deepcopy)
    @eval $(fun)(msa::AnnotatedMultipleSequenceAlignment) = AnnotatedMultipleSequenceAlignment($(fun)(msa.id), $(fun)(msa.msa), $(fun)(msa.annotations))
    @eval $(fun)(msa::MultipleSequenceAlignment) = MultipleSequenceAlignment($(fun)(msa.id), $(fun)(msa.msa))
    @eval $(fun)(seq::AnnotatedAlignedSequence) = AnnotatedAlignedSequence($(fun)(seq.id), $(fun)(seq.index), $(fun)(seq.sequence), $(fun)(seq.annotations))
    @eval $(fun)(seq::AlignedSequence) = AlignedSequence($(fun)(seq.id), $(fun)(seq.index), $(fun)(seq.sequence))
end

empty!(msa::AnnotatedMultipleSequenceAlignment) = (empty!(msa.id); empty!(msa.msa);
                                                   empty!(msa.annotations); msa)
empty!(msa::MultipleSequenceAlignment) = (empty!(msa.id); empty!(msa.msa); msa)
empty!(seq::AnnotatedAlignedSequence) = (empty!(seq.id); empty!(seq.index);
                                         v(seq.sequence); empty!(seq.annotations); seq)
empty!(seq::AlignedSequence) = (empty!(seq.id); empty!(seq.index);
                                empty!(seq.sequence); seq)

# Counting Gaps and Coverage
# --------------------------

"""
Calculates the fraction of gaps on the `Array` (alignment, sequence, column, etc.).
This function can take an extra `dim` argument for calculation of the gap fraction over the given dimension
"""
function gapfraction(x::AbstractArray{Residue})
    counter = 0
    len = 0
    for res in x
        counter += res == GAP ? 1 : 0
        len += 1
    end
    float(counter) / float(len)
end

gapfraction(x::AbstractArray{Residue},
              dim::Int) = vec( mapslices(gapfraction, x, dim) )

"""
Calculates the fraction of residues (no gaps) on the `Array` (alignment, sequence, column, etc.)
This function can take an extra `dim` argument for calculation of the residue fraction over the given dimension
"""
function residuefraction(x::AbstractArray{Residue})
    counter = 0
    len = 0
    for res in x
        counter += res == GAP ? 0 : 1
        len += 1
    end
    float(counter) / float(len)
end

residuefraction(x::AbstractArray{Residue},
                  dim::Int) = vec( mapslices(residuefraction, x, dim) )

"Coverage of the sequences with respect of the number of positions on the MSA"
coverage(msa::Matrix{Residue}) = residuefraction(msa, 2)
coverage(msa::AbstractMultipleSequenceAlignment) = coverage(msa.msa)

"Fraction of gaps per column/position on the MSA"
columngapfraction(msa::Matrix{Residue}) = gapfraction(msa, 1)
columngapfraction(msa::AbstractMultipleSequenceAlignment) = columngapfraction(msa.msa)

# Reference
# ---------

"""
Puts the sequence `i` as reference (as the first sequence) of the MSA.
This function swaps the sequences 1 and `i`, also an `id` can be used to select the sequence.
"""
function setreference!(msa::AnnotatedMultipleSequenceAlignment, i::Int, annotate::Bool=true)
    swap!(msa.id, 1, i)
    msa.msa[1, :], msa.msa[i, :] = msa.msa[i, :], msa.msa[1, :]
    annotate && annotate_modification!(msa, string("setreference! : Using ",
                                                   msa.id[1]," instead of ", msa.id[i]," as reference."))
    msa
end

function setreference!(msa::MultipleSequenceAlignment, i::Int, annotate::Bool=false) # annotate is useful for calling this inside other functions
    swap!(msa.id, 1, i)
    msa.msa[1, :], msa.msa[i, :] = msa.msa[i, :], msa.msa[1, :]
    msa
end

setreference!(msa::AbstractMultipleSequenceAlignment, id::ASCIIString,
              annotate::Bool=true) = setreference!(msa, findfirst(msa.id ,id), annotate)

function setreference!(msa::Matrix{Residue}, i::Int, annotate::Bool=false)
    msa[1, :], msa[i, :] = msa[i, :], msa[1, :]
    msa
end

"""
Creates a new Matrix{Residue}. This function deletes positions/columns of the MSA with gaps in the reference (first) sequence
"""
adjustreference(msa::Matrix{Residue},
                annotate::Bool=false) = msa[ : , vec(msa[1,:]) .!= GAP ] # annotate is useful for calling this inside other functions

"""
Removes positions/columns of the MSA with gaps in the reference (first) sequence
"""
adjustreference!(msa::AbstractMultipleSequenceAlignment,
                 annotate::Bool=true) = filtercolumns!(msa, vec(msa.msa[1,:]) .!= GAP, annotate)

"""
This functions deletes/filters sequences and columns/positions on the MSA on the following order:

 - Removes all the columns/position on the MSA with gaps on the reference sequence (first sequence)
 - Removes all the sequences with a coverage (with respect to the number of columns/positions on the MSA) **less** than a `coveragelimit` (default to `0.75`)
 - Removes all the columns/position on the MSA with **more** than a `gaplimit` (default to `0.5`: 50% of gaps)
"""
function gapstrip!(msa::AbstractMultipleSequenceAlignment, annotate::Bool=true;
                   coveragelimit::Float64=0.75, gaplimit::Float64=0.5)
    annotate && annotate_modification!(msa, string("gapstrip! : Deletes columns with gaps in the first sequence."))
    adjustreference!(msa, annotate)
    # Remove sequences with pour coverage of the reference sequence
    if ncolumns(msa) != 0
        annotate && annotate_modification!(msa, string("gapstrip! : Deletes sequences with a coverage less than ",
                                                       coveragelimit))
        filtersequences!(msa, coverage(msa) .>= coveragelimit, annotate)
    else
        throw("There are not columns in the MSA after the gap trimming")
    end
    # Remove columns with a porcentage of gap greater than gaplimit
    if nsequences(msa) != 0
        annotate && annotate_modification!(msa, string("gapstrip! : Deletes columns with more than ",
                                                       gaplimit, " gaps."))
        filtercolumns!(msa, columngapfraction(msa) .<= gaplimit, annotate)
    else
        throw("There are not sequences in the MSA after coverage filter")
    end
    msa
end

"""
Creates a new `Matrix{Residue}` with deleted sequences and columns/positions on the MSA:

 - Removes all the columns/position on the MSA with gaps on the reference sequence (first sequence)
 - Removes all the sequences with a coverage with respect to the number of columns/positions on the MSA **less** than a `coveragelimit` (default to `0.75`: sequences with 25% of gaps)
 - Removes all the columns/position on the MSA with **more** than a `gaplimit` (default to `0.5`: 50% of gaps)
"""
function gapstrip(msa::Matrix{Residue}; coveragelimit::Float64=0.75,
                  gaplimit::Float64=0.5)
    msa = adjustreference(msa)
    # Remove sequences with pour coverage of the reference sequence
    if ncolumns(msa) != 0
        msa = filtersequences(msa, coverage(msa) .>= coveragelimit )
    else
        throw("There are not columns in the MSA after the gap trimming")
    end
    if nsequences(msa) != 0
        msa = filtercolumns(msa, columngapfraction(msa) .<= gaplimit)
    else
        throw("There are not sequences in the MSA after coverage filter")
    end
    msa
end

# MSA getters
# -----------

"Returns the annotations of a MSA or a sequence."
annotations(msa::AnnotatedMultipleSequenceAlignment) = msa.annotations
annotations(seq::AnnotatedAlignedSequence) = seq.annotations

"""
`names(msa)`

Returns the `MSA` sequences names or identifiers as an `IndexedArray`.
"""
names(msa::AbstractMultipleSequenceAlignment) = msa.id
names(msa::Matrix{Residue}) = IndexedArrays.IndexedArray(ASCIIString[ string(i) for i in 1:nsequences(msa) ])

# Get annotations
# ---------------

for getter in [ :getannotcolumn, :getannotfile, :getannotresidue, :getannotsequence ]
    @eval $(getter)(msa::AnnotatedMultipleSequenceAlignment, args...) = $(getter)(msa.annotations, args...)
    @eval $(getter)(seq::AnnotatedAlignedSequence, args...) = $(getter)(seq.annotations, args...)
end

# Set annotations
# ---------------

for setter in [ :setannotcolumn!, :setannotfile!, :setannotresidue!, :setannotsequence!, :annotate_modification!, :delete_annotated_modifications!, :printmodifications ]
    @eval $(setter)(msa::AnnotatedMultipleSequenceAlignment, args...) = $(setter)(msa.annotations, args...)
    @eval $(setter)(seq::AnnotatedAlignedSequence, args...) = $(setter)(seq.annotations, args...)
end

# Used on AbstractMultipleSequenceAlignment methods
@inline annotate_modification!(msa::MultipleSequenceAlignment,
                               str::ASCIIString) = false # annotate_modification! is used on bool context: annotate && ...

# Show & Print
# ------------

"""
Gives an ASCIIString with the sequence number `seq` of the MSA
"""
asciisequence(msa::Matrix{Residue}, seq::Int) = ascii(convert(Vector{UInt8}, vec(msa[seq,:])))
asciisequence(msa::AbstractMultipleSequenceAlignment, seq::Int) = asciisequence(msa.msa, seq)
asciisequence(msa::AbstractMultipleSequenceAlignment, id::ASCIIString) = asciisequence(msa.msa, findfirst(msa.id, id))

# Mapping annotations
# ===================

"""
Converts a string of mappings into a vector of `Int`s

```
julia> _str2int_mapping(",,2,,4,5")
6-element Array{Int64,1}:
 0
 0
 2
 0
 4
 5

```
"""
function _str2int_mapping(mapping::ASCIIString)
    values = split(mapping, ',')
    len = length(values)
    intmap = Array(Int, len)
    @inbounds for i in 1:len
        value = values[i]
        intmap[i] = value == "" ? 0 : parse(Int, value)
    end
    intmap
end

getcolumnmapping(msa::AnnotatedMultipleSequenceAlignment) = _str2int_mapping(getannotfile(msa, "ColMap"))

"""
Returns the sequence coordinates as a `Vector{Int}` for an MSA sequence. That vector has one element for each MSA column.
If the number if `0` in the mapping, there is a gap in that column for that sequence.
"""
getsequencemapping(msa::AnnotatedMultipleSequenceAlignment,
                   seq_id::ASCIIString) = _str2int_mapping(getannotsequence(msa, seq_id, "SeqMap"))

getsequencemapping(msa::AnnotatedMultipleSequenceAlignment,
                   seq_num::Int) = getsequencemapping(msa, msa.id[seq_num])


text:
module BlackBoxOptim

using Distributions, StatsBase, Random, LinearAlgebra, Printf, Distributed, Compat
using SpatialIndexing
using Printf: @printf, @sprintf
using Compat: String, view

export  Optimizer, AskTellOptimizer, SteppingOptimizer, PopulationOptimizer,
        bboptimize, bbsetup, compare_optimizers,

        DiffEvoOpt, de_rand_1_bin, de_rand_1_bin_radiuslimited,
        adaptive_de_rand_1_bin, adaptive_de_rand_1_bin_radiuslimited,

        SeparableNESOpt, separable_nes,
        XNESOpt, xnes, dxnes,

        # Parameters
        DictChain, Parameters, ParamsDictChain, ParamsDict,
        chain, flatten,

        # Fitness
        FitnessScheme,
        ScalarFitnessScheme, ComplexFitnessScheme,
        MinimizingFitnessScheme, MaximizingFitnessScheme,
        IndexedTupleFitness, TupleFitnessScheme, ParetoFitnessScheme,
        EpsDominanceFitnessScheme, EpsBoxDominanceFitnessScheme,
        fitness_type, fitness_eltype, numobjectives,
        is_minimizing, nafitness, isnafitness,
        hat_compare, is_better, is_worse, same_fitness,
        aggregate,

        # Evaluator
        #ProblemEvaluator,

        # Problems
        Problems,
        OptimizationProblem, FunctionBasedProblem,
        minimization_problem,
        name, fitness_scheme_type, fitness_scheme, search_space, numdims, opt_value,
        fitness_is_within_ftol, objfunc, fitness,

        # Problem factory/family
        FunctionBasedProblemFamily, MinimizationProblemFamily,
        fixed_dim_problem, instantiate,

        save_fitness_history_to_csv_file,

        # Archive
        TopListArchive, EpsBoxArchive,
        best_fitness, best_candidate,
        last_top_fitness, delta_fitness, capacity,
        width_of_confidence_interval, fitness_improvement_potential,

        # OptimizationResults
        minimum, f_minimum, iteration_converged, parameters, population, pareto_frontier, params,
        archived_fitness,

        # OptController
        numruns, lastrun, problem,

        # Search spaces
        ParamBounds, Individual, SearchSpace, FixedDimensionSearchSpace, ContinuousSearchSpace,
        RangePerDimSearchSpace, symmetric_search_space,
        numdims, mins, maxs, deltas, ranges, range_for_dim, diameters,
        rand_individual, rand_individuals, rand_individuals_lhs,

        # Population
        FitPopulation,
        popsize,

        # Genetic operators
        GeneticOperator, MutationOperator, CrossoverOperator, EmbeddingOperator,
        NoMutation, MutationClock, GibbsMutationOperator, UniformMutation,
        PolynomialMutation,
        FixedGeneticOperatorsMixture, FAGeneticOperatorsMixture,
        RandomBound,
        SimpleSelector, RadiusLimitedSelector,
        SimulatedBinaryCrossover, SimplexCrossover, UnimodalNormalDistributionCrossover,
        ParentCentricCrossover,

        numparents, numchildren, apply!, adjust!,

        # Utilities
        FrequencyAdapter, update!, frequencies,
        name

module Utils
    using Random

    include("utilities/latin_hypercube_sampling.jl")
    include("utilities/assign_ranks.jl")
    include("utilities/halton_sequence.jl")
end

const SI = SpatialIndexing

include("search_space.jl")
include("parameters.jl")
include("fitness.jl")
include("ntuple_fitness.jl")

include("problem.jl")

include("frequency_adaptation.jl")

include("fit_individual.jl")
include("archive.jl")
include("archives/dominance_cone.jl")
include("archives/epsbox_archive.jl")

include("genetic_operators/genetic_operator.jl")

include("evaluator.jl")
include("parallel_evaluator.jl")

include("population.jl")
include("optimizer.jl")

# Different optimization algorithms/methods
include("random_search.jl")
include("differential_evolution.jl")
include("adaptive_differential_evolution.jl")
include("natural_evolution_strategies.jl")
include("dx_nes.jl")
include("resampling_memetic_search.jl")
include("simultaneous_perturbation_stochastic_approximation.jl")
include("generating_set_search.jl")
include("direct_search_with_probabilistic_descent.jl")

# multi-objective optimization algorithms
include("borg_moea.jl")

# End-user/top-level interface functions
include(joinpath("problems", "problem_family.jl"))
include("optimization_methods.jl")
include("default_parameters.jl")
include("optimization_result.jl")
include("opt_controller.jl")
include("bboptimize.jl")
include("compare_optimizers.jl")

# Problems for testing
include(joinpath("problems", "single_objective.jl"))
include(joinpath("problems", "multi_objective.jl"))

end # module BlackBoxOptim


text:
"""
    SingleBridgeOptimizer{BT<:AbstractBridge, OT<:MOI.ModelLike} <:
    AbstractBridgeOptimizer

The `SingleBridgeOptimizer` bridges any constraint supported by the bridge `BT`.
This is in contrast with the [`MathOptInterface.Bridges.LazyBridgeOptimizer`](@ref)
which only bridges the constraints that are unsupported by the internal model,
even if they are supported by one of its bridges.
"""
mutable struct SingleBridgeOptimizer{BT<:AbstractBridge,OT<:MOI.ModelLike} <:
               MOIB.AbstractBridgeOptimizer
    model::OT
    map::Map # index of bridged constraint -> constraint bridge
    con_to_name::Dict{MOI.ConstraintIndex,String}
    name_to_con::Union{Dict{String,MOI.ConstraintIndex},Nothing}
end

function SingleBridgeOptimizer{BT}(model::OT) where {BT,OT<:MOI.ModelLike}
    return SingleBridgeOptimizer{BT,OT}(
        model,
        Map(),
        Dict{MOI.ConstraintIndex,String}(),
        nothing,
    )
end

function bridges(bridge::MOI.Bridges.AbstractBridgeOptimizer)
    return EmptyMap()
end

function bridges(bridge::SingleBridgeOptimizer)
    return bridge.map
end

MOIB.supports_constraint_bridges(::SingleBridgeOptimizer) = true

# If `BT` bridges `MOI.Reals` (such as `Constraint.FunctionizeBridge` bridge,
# without this method, it creates a `StackOverflow` with
# `is_bridged`, `supports_bridging_constrained_variable`
# and `supports_add_constrained_variables`.
MOIB.is_bridged(::SingleBridgeOptimizer, ::Type{MOI.Reals}) = false

function MOIB.is_bridged(b::SingleBridgeOptimizer, S::Type{<:MOI.AbstractSet})
    return MOIB.supports_bridging_constrained_variable(b, S)
end

MOIB.is_bridged(::SingleBridgeOptimizer, ::MOI.VariableIndex) = false
function MOIB.is_bridged(
    b::SingleBridgeOptimizer,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,S},
) where {S}
    return MOIB.is_bridged(b, MOI.VariableIndex, S) &&
           haskey(Constraint.bridges(b), ci)
end
function MOIB.is_bridged(
    b::SingleBridgeOptimizer,
    ci::MOI.ConstraintIndex{MOI.VectorOfVariables,S},
) where {S}
    return MOIB.is_bridged(b, MOI.VectorOfVariables, S) &&
           haskey(Constraint.bridges(b), ci)
end

function MOIB.supports_bridging_constrained_variable(
    b::SingleBridgeOptimizer,
    S::Type{<:MOI.AbstractSet},
)
    return MOIB.supports_bridging_constraint(
        b,
        MOIU.variable_function_type(S),
        S,
    ) && MOI.supports_add_constrained_variables(b, MOI.Reals)
end

function MOIB.supports_bridging_constraint(
    ::SingleBridgeOptimizer{BT},
    F::Type{<:MOI.AbstractFunction},
    S::Type{<:MOI.AbstractSet},
) where {BT}
    return MOI.supports_constraint(BT, F, S)
end

function MOIB.is_bridged(
    b::SingleBridgeOptimizer,
    F::Type{<:MOI.AbstractFunction},
    S::Type{<:MOI.AbstractSet},
)
    return MOIB.supports_bridging_constraint(b, F, S)
end

function MOIB.is_bridged(
    ::SingleBridgeOptimizer,
    ::Type{<:MOI.AbstractScalarFunction},
)
    return false
end

function MOIB.bridge_type(
    ::SingleBridgeOptimizer{BT},
    ::Type{<:MOI.AbstractFunction},
    ::Type{<:MOI.AbstractSet},
) where {BT}
    return BT
end

MOIB.bridging_cost(::SingleBridgeOptimizer, args...) = 1.0

MOIB.recursive_model(b::SingleBridgeOptimizer) = b.model


text:
using CubicEoS
using CubicEoSDatabase: Data

c1c5mix = load(
    BrusilovskyEoSMixture,
    names = ("methane", "n-pentane"),
    component_dbs = (Data.martinez(), Data.brusilovsky_comp()),
    mix_eos_db = Data.brusilovsky_mix()
)

"""
Prints (Ση, T)-diagram (variables used in (Mikyska, 2012)) of single phase stability region.

Uses VT-stability solver.

The mixture of composition `χ_base` is loaded by `names` of its components.
Then calculation performed for each (T ∈ `T_range`, Ση_base ∈ `Ση_range`) point.

Warns through stderr if an exception happens.
"""
function vt_stability_η_temp(io::IO, χ_base, Ση_range, T_range)
    mixture = c1c5mix

    let
        prefix = "# "
        println(io, prefix, "VT-stability calculation of (Ση, T)-diagram")
        println(io, prefix, "Mixture: $(name(mixture))")
        println(io, prefix, "Composition (in molar parts) of base phase χ_base: ")
        for (c, χ) in zip(components(mixture), χ_base)
            println(io, prefix, "  $(name(c))\t\t$χ")
        end
        println(io, prefix, "Range of overall concentration Ση (mol m⁻³):")
        println(io, prefix, "  from: ", first(Ση_range))
        println(io, prefix, "    to: ", last(Ση_range))
        println(io, prefix, "Range of temperatures T (K):")
        println(io, prefix, "  from: ", first(T_range))
        println(io, prefix, "    to: ", last(T_range))

        columns = (
            ("T", "temperature (К)"),
            ("Ση", "overall concentration of mixture (mol m⁻³)"),
            # generator of tuples ("η_base_i", "description of η_base_i")
            (("η_base_$i", "concentration (mol m⁻³) of $(name(c)) in base phase") for (i, c) in enumerate(components(mixture)))...,
            ("is_fail", "if true, an error raised during calculation"),
            ("is_stable", "if true, base phase is stable (single phase state)"),
        )
        println(io, prefix, "Data representation")
        println(io, prefix, "  <column index>. <symbol>: <description>")
        for (i, col) in enumerate(columns)
            println(io, prefix, "  $i. ", col[1], ": ", col[2])
        end
    end

    sep = '\t'

    η_base = fill(NaN, length(mixture))

    for T in T_range
        for Ση in Ση_range
            η_base .= Ση .* χ_base
            RT = T * CubicEoS.GAS_CONSTANT_SI

            is_fail = false
            is_stable = false

            try
                is_stable, = vt_stability(mixture, η_base, 1.0, RT)
            catch e
                @warn "Exception occurs" exception=e η_base=repr(η_base) RT=RT
                is_fail = true
                is_stable = false
            finally
                println(io, join([T, Ση, η_base..., repr(is_fail), repr(is_stable)], sep))
            end
        end
    end
end

vt_stability_η_temp(χ_base, Ση_range, T_range) = vt_stability_η_temp(stdout, χ_base, Ση_range, T_range)

function vt_stability_η_temp(fname::AbstractString, χ_base, Ση_range, T_range)
    open(fname, "w") do io
        vt_stability_η_temp(io, χ_base, Ση_range, T_range)
    end
end

vt_stability_η_temp([0.547413, 0.452587], 100:10:15000, 300:1:450)

max_stars_repo_name: vvpisarev/CubicEoS

max_issues_repo_name: vvpisarev/CubicEoS

max_forks_repo_name: vvpisarev/CubicEoS


text:
export param, param0
using AutoGrad: Param
using Knet: atype

"""
    param(array; atype)
    param(dims...; init, atype)
    param0(dims...; atype)

The first form returns `Param(atype(array))`.

The second form Returns a randomly initialized `Param(atype(init(dims...)))`.

The third form `param0` is an alias for `param(dims...; init=zeros)`.

By default, `init` is `xavier_uniform` and `atype` is `Knet.atype()`.

"""
param,param0

# TODO: Knet.Param <: AutoGrad.Tracked as a separate type?
param(x::AbstractArray; atype=atype()) = Param(convert(atype,x))
param(d...; init=xavier_uniform, atype=atype())=Param(atype(init(d...)))
param0(d...; atype=atype())=param(d...; init=zeros, atype=atype)


text:
using Test, JuMP
const MOIT = MOI.Test
const MOIB = MOI.Bridges

function _model(optimizer::MOI.AbstractOptimizer)
    MOI.empty!(optimizer)
    return direct_model(optimizer)
end

function _model(optimizer_constructor)
    return Model(optimizer_constructor)
end

#"""
#    @test_suite setname subsets
#
#Defines a function `setname_test(model, config, exclude)` that runs the tests
#defined in the dictionary `setname_tests` with the model `model` and config
#`config` except the tests whose dictionary key is in `exclude`. If `subsets` is
#`true` then each test runs in fact multiple tests hence the `exclude` argument
#is passed as it can also contains test to be excluded from these subsets of
#tests.
#"""
macro test_suite(setname, subsets=false)
    testname = Symbol(string(setname) * "_test")
    testdict = Symbol(string(testname) * "s")
    if subsets
        runtest = :( f(model, config, exclude) )
    else
        runtest = :( f(model, config) )
    end
    esc(:(
      function $testname(model,
                         config::$MOI.Test.Config,
                         exclude::Vector{String} = String[])
            for (name,f) in $testdict
                if name in exclude
                    continue
                end
                @testset "$name" begin
                    $runtest
                end
            end
        end
    ))
end

# Utilities for building the mock `optimize!` from the solution of a solver
# Variables primal values for inner bridged model
function print_value(v, atol)
    i = round(v)
    if isapprox(v, i, atol=atol)
        print(float(i))
    else
        print(v)
    end
end
function inner_variable_value(model, atol=1e-4)
    inner = backend(model)
    println("optimize!(mock) = MOIU.mock_optimize!(mock,")
    println(JuMP.termination_status(model))
    if JuMP.primal_status(model) != MOI.NO_SOLUTION
        values = MOI.get(inner, MOI.VariablePrimal(),
                         MOI.get(inner, MOI.ListOfVariableIndices()))
        print("(MOI.FEASIBLE_POINT, [")
        for (i, v) in enumerate(values)
            if i > 1
                print(", ")
            end
            print_value(v, atol)
        end
        print("])")
    else
        print("MOI.NO_SOLUTION")
    end
    println(",")
    if JuMP.dual_status(model) != MOI.NO_SOLUTION
        for (F, S) in MOI.get(inner, MOI.ListOfConstraints())
            print("($F, $S) => [")
            for ci in MOI.get(inner, MOI.ListOfConstraintIndices{F, S}())
                print(MOI.get(inner, MOI.ConstraintDual(), ci))
                print(", ")
            end
            println("])")
        end
    end
    println(")")
end
# Constraint dual values for inner bridged model


text:
"""
Perform a union between two HyperRectangles.
"""
function union{T,N}(h1::HyperRectangle{N, T}, h2::HyperRectangle{N, T})
    m = min(minimum(h1), minimum(h2))
    mm = max(maximum(h1), maximum(h2))
    HyperRectangle{N, T}(m, mm-m)
end

"""
Perform a difference between two HyperRectangles.
"""
diff(h1::HyperRectangle, h2::HyperRectangle) = h1

"""
Perform a intersection between two HyperRectangles.
"""
function intersect{T,N}(h1::HyperRectangle{N, T}, h2::HyperRectangle{N, T})
    m = max(minimum(h1), minimum(h2))
    mm = min(maximum(h1), maximum(h2))
    HyperRectangle{N, T}(m, mm-m)
end

function intersect(a::SimpleRectangle, b::SimpleRectangle)
    min_n = max(minimum(a), minimum(b))
    max_n = min(maximum(a), maximum(b))
    w = max_n - min_n
    SimpleRectangle(min_n[1], min_n[2], w[1], w[2])
end


text:
struct DecisionBoundaryPQs <: ModelBasedPQs
    occ::SVDD.OCClassifier
end

decion_boundary_pqs(predictions) = -abs.(predictions)

function qs_score(qs::DecisionBoundaryPQs,
                  x::Array{T, 2},
                  labels::Dict{Symbol, Array{Int, 1}})::Array{Float64, 1} where {T <: Real}
    return decion_boundary_pqs(SVDD.predict(qs.occ, x))
end

function qs_score(qs::DecisionBoundaryPQs,
                  x::Array{T, 2},
                  labels::Dict{Symbol, Array{Int, 1}},
                  subspaces::Vector{Vector{Int}}) where {T <: Real}
    size(x,1) == size(qs.occ.data,1) || throw(DimensionMismatch("Number of dimensions in x must match the number of dimensions of the classifier."))
    return map(idx -> decion_boundary_pqs(SVDD.predict(qs.occ, x[subspaces[idx], :], idx)), eachindex(subspaces))
end