text:
"""
Creates lossy converter model between AC and DC grid
```
pconv_ac[i] + pconv_dc[i] == a + bI + cI^2
```
"""
function constraint_converter_losses(pm::AbstractACPModel, n::Int, cnd::Int, i::Int, a, b, c, plmax)
    pconv_ac = PowerModels.var(pm, n, cnd, :pconv_ac, i)
    pconv_dc = PowerModels.var(pm, n, cnd, :pconv_dc, i)
    iconv = PowerModels.var(pm, n, cnd, :iconv_ac, i)

    PowerModels.con(pm, n, cnd, :conv_loss)[i] = @NLconstraint(pm.model, pconv_ac + pconv_dc == a + b*iconv + c*iconv^2)
end
"""
Links converter power & current
```
pconv_ac[i]^2 + pconv_dc[i]^2 == vmc[i]^2 * iconv_ac[i]^2
```
"""
function constraint_converter_current(pm::AbstractACPModel, n::Int, cnd::Int, i::Int, Umax, Imax)
    vmc = PowerModels.var(pm, n, cnd, :vmc, i)
    pconv_ac = PowerModels.var(pm, n, cnd, :pconv_ac, i)
    qconv_ac = PowerModels.var(pm, n, cnd, :qconv_ac, i)
    iconv = PowerModels.var(pm, n, cnd, :iconv_ac, i)

    PowerModels.con(pm, n, cnd, :conv_i)[i] = @NLconstraint(pm.model, pconv_ac^2 + qconv_ac^2 == vmc^2 * iconv^2)
end
"""
Converter transformer constraints
```
p_tf_fr ==  g/(tm^2)*vm_fr^2 + -g/(tm)*vm_fr*vm_to * cos(va_fr-va_to) + -b/(tm)*vm_fr*vm_to*sin(va_fr-va_to)
q_tf_fr == -b/(tm^2)*vm_fr^2 +  b/(tm)*vm_fr*vm_to * cos(va_fr-va_to) + -g/(tm)*vm_fr*vm_to*sin(va_fr-va_to)
p_tf_to ==  g*vm_to^2 + -g/(tm)*vm_to*vm_fr  *    cos(va_to - va_fr)     + -b/(tm)*vm_to*vm_fr    *sin(va_to - va_fr)
q_tf_to == -b*vm_to^2 +  b/(tm)*vm_to*vm_fr  *    cos(va_to - va_fr)     + -g/(tm)*vm_to*vm_fr    *sin(va_to - va_fr)
```
"""
function constraint_conv_transformer(pm::AbstractACPModel, n::Int, cnd::Int, i::Int, rtf, xtf, acbus, tm, transformer)
    ptf_fr = PowerModels.var(pm, n, cnd, :pconv_tf_fr, i)
    qtf_fr = PowerModels.var(pm, n, cnd, :qconv_tf_fr, i)
    ptf_to = PowerModels.var(pm, n, cnd, :pconv_tf_to, i)
    qtf_to = PowerModels.var(pm, n, cnd, :qconv_tf_to, i)

    vm = PowerModels.var(pm, n, cnd, :vm, acbus)
    va = PowerModels.var(pm, n, cnd, :va, acbus)
    vmf = PowerModels.var(pm, n, cnd, :vmf, i)
    vaf = PowerModels.var(pm, n, cnd, :vaf, i)

    ztf = rtf + im*xtf
    if transformer
        ytf = 1/(rtf + im*xtf)
        gtf = real(ytf)
        btf = imag(ytf)
        gtf_sh = 0
        c1, c2, c3, c4 = ac_power_flow_constraints(pm.model, gtf, btf, gtf_sh, vm, vmf, va, vaf, ptf_fr, ptf_to, qtf_fr, qtf_to, tm)

        PowerModels.con(pm, n, cnd, :conv_tf_p_fr)[i] = c1
        PowerModels.con(pm, n, cnd, :conv_tf_q_fr)[i] = c2

        PowerModels.con(pm, n, cnd, :conv_tf_p_to)[i] = c3
        PowerModels.con(pm, n, cnd, :conv_tf_q_to)[i] = c4
    else
        PowerModels.con(pm, n, cnd, :conv_tf_p_fr)[i] = @constraint(pm.model, ptf_fr + ptf_to == 0)
        PowerModels.con(pm, n, cnd, :conv_tf_q_fr)[i] = @constraint(pm.model, qtf_fr + qtf_to == 0)
        @constraint(pm.model, va == vaf)
        @constraint(pm.model, vm == vmf)
    end
end
"constraints for a voltage magnitude transformer + series impedance"
function ac_power_flow_constraints(model, g, b, gsh_fr, vm_fr, vm_to, va_fr, va_to, p_fr, p_to, q_fr, q_to, tm)
    c1 = @NLconstraint(model, p_fr ==  g/(tm^2)*vm_fr^2 + -g/(tm)*vm_fr*vm_to * cos(va_fr-va_to) + -b/(tm)*vm_fr*vm_to*sin(va_fr-va_to))
    c2 = @NLconstraint(model, q_fr == -b/(tm^2)*vm_fr^2 +  b/(tm)*vm_fr*vm_to * cos(va_fr-va_to) + -g/(tm)*vm_fr*vm_to*sin(va_fr-va_to))
    c3 = @NLconstraint(model, p_to ==  g*vm_to^2 + -g/(tm)*vm_to*vm_fr  *    cos(va_to - va_fr)     + -b/(tm)*vm_to*vm_fr    *sin(va_to - va_fr))
    c4 = @NLconstraint(model, q_to == -b*vm_to^2 +  b/(tm)*vm_to*vm_fr  *    cos(va_to - va_fr)     + -g/(tm)*vm_to*vm_fr    *sin(va_to - va_fr))
    return c1, c2, c3, c4
end
"""
Converter reactor constraints
```
-pconv_ac == gc*vmc^2 + -gc*vmc*vmf*cos(vac-vaf) + -bc*vmc*vmf*sin(vac-vaf)
-qconv_ac ==-bc*vmc^2 +  bc*vmc*vmf*cos(vac-vaf) + -gc*vmc*vmf*sin(vac-vaf)
p_pr_fr ==  gc *vmf^2 + -gc *vmf*vmc*cos(vaf - vac) + -bc *vmf*vmc*sin(vaf - vac)
q_pr_fr == -bc *vmf^2 +  bc *vmf*vmc*cos(vaf - vac) + -gc *vmf*vmc*sin(vaf - vac)
```
"""
function constraint_conv_reactor(pm::AbstractACPModel, n::Int, cnd::Int, i::Int, rc, xc, reactor)
    pconv_ac = PowerModels.var(pm, n, cnd, :pconv_ac, i)
    qconv_ac = PowerModels.var(pm, n, cnd, :qconv_ac, i)
    ppr_fr = PowerModels.var(pm, n, cnd, :pconv_pr_fr, i)
    qpr_fr = PowerModels.var(pm, n, cnd, :qconv_pr_fr, i)

    vmf = PowerModels.var(pm, n, cnd, :vmf, i)
    vaf = PowerModels.var(pm, n, cnd, :vaf, i)
    vmc = PowerModels.var(pm, n, cnd, :vmc, i)
    vac = PowerModels.var(pm, n, cnd, :vac, i)

    zc = rc + im*xc
    if reactor
        yc = 1/(zc)
        gc = real(yc)
        bc = imag(yc)
        PowerModels.con(pm, n, cnd, :conv_pr_p)[i] = @NLconstraint(pm.model, -pconv_ac == gc*vmc^2 + -gc*vmc*vmf*cos(vac-vaf) + -bc*vmc*vmf*sin(vac-vaf))
        PowerModels.con(pm, n, cnd, :conv_pr_q)[i] = @NLconstraint(pm.model, -qconv_ac ==-bc*vmc^2 +  bc*vmc*vmf*cos(vac-vaf) + -gc*vmc*vmf*sin(vac-vaf))
        @NLconstraint(pm.model, ppr_fr ==  gc *vmf^2 + -gc *vmf*vmc*cos(vaf - vac) + -bc *vmf*vmc*sin(vaf - vac))
        @NLconstraint(pm.model, qpr_fr == -bc *vmf^2 +  bc *vmf*vmc*cos(vaf - vac) + -gc *vmf*vmc*sin(vaf - vac))
    else
        ppr_to = -pconv_ac
        qpr_to = -qconv_ac
        PowerModels.con(pm, n, cnd, :conv_pr_p)[i] = @constraint(pm.model, ppr_fr + ppr_to == 0)
        PowerModels.con(pm, n, cnd, :conv_pr_q)[i] = @constraint(pm.model, qpr_fr + qpr_to == 0)
        @constraint(pm.model, vac == vaf)
        @constraint(pm.model, vmc == vmf)

    end
end
"""
Converter filter constraints
```
ppr_fr + ptf_to == 0
qpr_fr + qtf_to +  (-bv) * filter *vmf^2 == 0
```
"""
function constraint_conv_filter(pm::AbstractACPModel, n::Int, cnd::Int, i::Int, bv, filter)
    ppr_fr = PowerModels.var(pm, n, cnd, :pconv_pr_fr, i)
    qpr_fr = PowerModels.var(pm, n, cnd, :qconv_pr_fr, i)
    ptf_to = PowerModels.var(pm, n, cnd, :pconv_tf_to, i)
    qtf_to = PowerModels.var(pm, n, cnd, :qconv_tf_to, i)

    vmf = PowerModels.var(pm, n, cnd, :vmf, i)

    PowerModels.con(pm, n, cnd, :conv_kcl_p)[i] = @constraint(pm.model,   ppr_fr + ptf_to == 0 )
    PowerModels.con(pm, n, cnd, :conv_kcl_q)[i] = @NLconstraint(pm.model, qpr_fr + qtf_to +  (-bv) * filter *vmf^2 == 0)
end
"""
LCC firing angle constraints
```
pconv_ac == cos(phi) * Srated
qconv_ac == sin(phi) * Srated
```
"""
function constraint_conv_firing_angle(pm::AbstractACPModel, n::Int, cnd::Int, i::Int, S, P1, Q1, P2, Q2)
    p = PowerModels.var(pm, n, cnd, :pconv_ac, i)
    q = PowerModels.var(pm, n, cnd, :qconv_ac, i)
    phi = PowerModels.var(pm, n, cnd, :phiconv, i)

    PowerModels.con(pm, n, cnd, :conv_cosphi)[i] = @NLconstraint(pm.model,   p == cos(phi) * S)
    PowerModels.con(pm, n, cnd, :conv_sinphi)[i] = @NLconstraint(pm.model,   q == sin(phi) * S)
end

function constraint_dc_droop_control(pm::AbstractACPModel, n::Int, cnd::Int, i::Int, busdc_i, vref_dc, pref_dc, k_droop)
    pconv_dc = PowerModels.var(pm, n, cnd, :pconv_dc, i)
    vdc = PowerModels.var(pm, n, cnd, :vdcm, busdc_i)

    PowerModels.con(pm, n, cnd, :conv_dc_droop)[i] = @constraint(pm.model, pconv_dc == pref_dc - sign(pref_dc) * 1 / k_droop * (vdc - vref_dc))
end


text:
if CUDA.functional()
    onGpu = gpu
else
    onGpu = identity
end
@testset "ConvFFT transform" begin
    @testset "ConvFFT 2D" begin
        originalSize = (10,10,1,2)
        tmp = zeros(originalSize);
        init = zeros(originalSize) |> onGpu
        init[5,5,1,2] = Float32(1)
        shears = ConvFFT(originalSize) |> onGpu
        res = shears(init);
        @test size(res) == (10,10,5,1,2)
        # TODO: for the other boundary conditions. This is just periodic
        function minimalTransform(shears, init)
            equivalent = zeros(10, 10, 5, 1, 2) |> onGpu
            for i=1:5
                equivalent[:,:,i,:,:] = irfft(rfft(init,(1,2)) .*
                                              shears.weight[:,:,i], 10, (1,2)) .+ shears.bias[:,i]
            end
            return equivalent
        end
        #@info "" minimum(res), minimum(minimalTransform(shears, init))
        @test minimalTransform(shears, init) ≈ res

        shears = ConvFFT(originalSize, 5, abs) |> onGpu
        res = shears(init);
        @test abs.(minimalTransform(shears, init)) ≈ res
    end
    @testset "ConvFFT 1D" begin
        originalSize = (10,1,2)
        tmp = zeros(originalSize); tmp
        init = zeros(originalSize) |> onGpu
        init[5,1,2] = Float32(1)
        shears = ConvFFT(originalSize, nConvDims=1, boundary=Pad(-1)) |> onGpu
        res = shears(init);
        @test size(res) == (10,5,1,2)
        # TODO: this is only padded
        function minimalTransform(shears, init)
            equivalent = zeros(16, 5, 1, 2) |> onGpu
            for i=1:5
                equivalent[:,i,:,:] = irfft(rfft(pad(init,
                                                     shears.bc.padBy),
                                                 (1,)) .*
                                            shears.weight[:,i], 16, (1,)) .+ shears.bias[:,i]
            end
            return equivalent[4:13, :, :, :]
        end
        @test minimalTransform(shears, init) ≈ res

        shears = ConvFFT(originalSize, 5, abs, nConvDims=1,
                         boundary=Pad(-1)) |> onGpu
        res = shears(init);
        @test abs.(minimalTransform(shears, init)) ≈ res
    end
end

ext: jl


text:
module CP

include("CP_Exact.jl")
include("CP_EM.jl")
include("CP_VB.jl")
include("CP_MC.jl")
include("BackwardKernel.jl")

using .CP_Exact, .CP_EM, .CP_VB, .CP_MC, .BackwardKernel

export generate, log_marginal
export standard_EM, dual_EM
export standard_VB, online_VB
export Particle, smc_weight, particle_mcmc
export EventQueue, sum, length, eltype, iterate

end

max_stars_repo_name: mehmetburakkurutmaz/bam

max_issues_repo_name: mehmetburakkurutmaz/bam

max_forks_repo_name: atcemgil/bam


text:
using TensorNetworkAD
using Test
using Zygote, OMEinsum
using TensorNetworkAD: model_tensor, tensorfromclassical, trg_svd

@testset "exampletensor" begin
    β = rand()
    @test model_tensor(Ising(),β) ≈ tensorfromclassical([β -β; -β β])
end


text:
using ..JuliaRL.FeatureCreators
import ..build_algorithm_dict
import Flux: Descent

"""
    Agent for learning online in the four rooms environment.
"""
mutable struct TCFourRoomsContAgent{O, P<:Resampling.AbstractPolicy} <: JuliaRL.AbstractAgent
    μ::P
    ER::ExperienceReplay
    WER::WeightedExperienceReplay
    opt::O
    tilecoder::TileCoder
    gvf::GVF
    training_gap::Int64
    buffer_size::Int64
    batch_size::Int64
    warm_up::Int64
    algo_dict::Dict{String, Resampling.LearningUpdate}
    sample_dict::Dict{String, String}
    value_dict::Dict{String, Array{Resampling.SparseLayer, 1}}
    counter::Int64
    state_t::Tuple{Array{Float64, 1}, Bool}
    action_t::Int64
    α_arr::Array{Float64, 1}
    max_is::Float64
    function TCFourRoomsContAgent(μ::P, gvf::GVF, opt::O, num_tilings, num_tiles, α_arr::Array{Float64, 1},
                                  training_gap, buffer_size, batch_size,
                                  warm_up, parsed, env_size; max_is_ratio=1.0) where {O, P<:Resampling.AbstractPolicy}

        # build tile coder....
        tilecoder = TileCoder(num_tilings, num_tiles, 2)
        num_features = num_tilings*(num_tiles+1)^2
        algo_dict, sample_dict, value_type_dict = build_algorithm_dict(parsed; max_is=max_is_ratio, size_features=num_features)
        # println(algo_dict)
        value_dict = Dict{String, Array{Resampling.SparseLayer, 1}}()
        for key in keys(algo_dict)
            if value_type_dict[key] == "State"
                value_dict[key] = [SparseLayer(num_features, 1) for α in α_arr]
            else
                value_dict[key] = [SparseLayer(num_features, 4) for α in α_arr]
            end
        end
        er_names = [:s_t, :s_tp1, :a_t, :a_tp1, :r, :γ_tp1, :terminal, :mu_t, :pi_t, :ρ]
        state_type = Array{Int64, 1}
        er_types = [state_type, state_type, Int64, Int64, Float64, Float64, Bool, Float64, Float64, Float64]
        ER = ExperienceReplay(buffer_size, er_types; column_names=er_names)
        WER = WeightedExperienceReplay(buffer_size, er_types; column_names=er_names)
        new{O, P}(μ, ER, WER, opt, tilecoder, gvf, training_gap, buffer_size, batch_size, warm_up,
                  algo_dict, sample_dict, value_dict, warm_up, ([0.0,0.0], false), 0, α_arr,
                  max_is_ratio)
    end
end

JuliaRL.get_action(agent::TCFourRoomsContAgent, s; rng=Random.GLOBAL_RNG) =
    StatsBase.sample(rng, agent.μ, s[1])

function JuliaRL.start!(agent::TCFourRoomsContAgent, env_s_tp1; rng=Random.GLOBAL_RNG, kwargs...)
    agent.counter = agent.warm_up
    agent.state_t = env_s_tp1
    agent.action_t = JuliaRL.get_action(agent, env_s_tp1; rng=rng)
    return agent.action_t
end

function JuliaRL.step!(agent::TCFourRoomsContAgent, env_s_tp1, r, terminal; rng=Random.GLOBAL_RNG, kwargs...)

    # After taking action at timestep t
    μ_prob = get(agent.μ, agent.state_t[1], agent.action_t, nothing, nothing, nothing)
    c, γ, π_prob = get(agent.gvf, agent.state_t, agent.action_t, env_s_tp1, nothing, nothing)
    action_tp1 = JuliaRL.get_action(agent, agent.state_t; rng=rng)

    experience = (create_features(agent.tilecoder, agent.state_t[1]./11.0),
                  create_features(agent.tilecoder, env_s_tp1[1]./11.0),
                  agent.action_t, action_tp1, c, γ, terminal, μ_prob, π_prob, (π_prob/μ_prob))

    add!(agent.ER, experience)
    add!(agent.WER, experience, π_prob/μ_prob)

    if agent.counter == 0
        agent.counter = agent.training_gap
        train_value_functions(agent; rng=rng)
    end
    agent.counter -= 1

    agent.state_t = env_s_tp1
    agent.action_t = action_tp1 #JuliaRL.get_action(agent, agent.state_t; rng=rng)
    return agent.action_t
end

function train_value_functions(agent::TCFourRoomsContAgent; rng=Random.GLOBAL_RNG)

    α_arr = agent.α_arr

    samp_er = sample(agent.ER, agent.batch_size; rng=rng)
    samp_wer = sample(agent.WER, agent.batch_size; rng=rng)

    arg_er = (samp_er[:ρ], samp_er[:s_t], samp_er[:s_tp1],
              samp_er[:r], samp_er[:γ_tp1], samp_er[:terminal],
              samp_er[:a_t], samp_er[:a_tp1], agent.gvf.policy)

    arg_wer = (ones(length(samp_wer[:ρ])), samp_wer[:s_t], samp_wer[:s_tp1],
               samp_wer[:r], samp_wer[:γ_tp1], samp_wer[:terminal],
               samp_wer[:a_t], samp_wer[:a_tp1], agent.gvf.policy)

    for (α_idx, α) in enumerate(α_arr)
        agent.opt.eta = α
        for key in keys(agent.algo_dict)
            if agent.sample_dict[key] == "IR"
                corr_term = 1.0
                if key == "BCIR"
                    corr_term = Resampling.total(agent.WER)/size(agent.WER)[1]
                end
                update!(agent.value_dict[key][α_idx], agent.opt, agent.algo_dict[key], arg_wer...; corr_term=corr_term)
            elseif agent.sample_dict[key] == "ER"
                if key == "WISBuffer"
                    corr_term = Resampling.total(agent.WER)/size(agent.WER)[1]
                    update!(agent.value_dict[key][α_idx], agent.opt, agent.algo_dict[key], arg_er...; corr_term=1.0/corr_term)
                elseif key == "NormIS"
                    update!(agent.value_dict[key][α_idx], agent.opt, agent.algo_dict[key], arg_er[1]./agent.max_is, arg_er[2:end]...;)
                else
                    update!(agent.value_dict[key][α_idx], agent.opt, agent.algo_dict[key], arg_er...;)
                end
            elseif agent.sample_dict[key] == "Optimal"
                samp_opt = agent.ER.buffer
                arg_opt = (samp_opt[:ρ], samp_opt[:s_t], samp_opt[:s_tp1],
                           samp_opt[:r], samp_opt[:γ_tp1], samp_opt[:terminal])
                update!(agent.value_dict[key][α_idx], agent.opt, agent.algo_dict[key], arg_opt...;)
            end
        end
    end
end

function predict(agent::TCFourRoomsContAgent, key, α_idx, eval_states::Array{Array{Float64, 1}, 1})
    states = [create_features(agent.tilecoder, eval_state./11.0) for eval_state in eval_states]
    return agent.value_dict[key][α_idx].(states)
end

function predict!(agent::TCFourRoomsContAgent, eval_states::Array{Array{Float64, 1}, 1}, predict_dict::Dict{String, Array{Array{Float64, 1}}})
    states = [create_features(agent.tilecoder, eval_state./11.0) for eval_state in eval_states]
    if length(keys(predict_dict)) == 0
        for key in keys(agent.value_dict)
            predict_dict[key] = [zeros(length(eval_states)) for i in 1:length(agent.α_arr)]
        end
    end
    for key in keys(agent.value_dict)
        for α_idx in 1:length(agent.α_arr)
            predict_dict[key][α_idx] .= agent.value_dict[key][α_idx].(states)
        end
    end
end


text:
module SingleFieldFESpacesTests

using Test
using Gridap.Arrays
using Gridap.TensorValues
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.FESpaces

domain =(0,1,0,1,0,1)
partition = (3,3,3)
model = CartesianDiscreteModel(domain,partition)

order = 2
grid_topology = get_grid_topology(model)
polytopes = get_polytopes(grid_topology)
reffes = [LagrangianRefFE(Float64,p,order) for p in polytopes]

dirichlet_tags = ["tag_24","tag_25"]
V0 = GradConformingFESpace(reffes,model,dirichlet_tags)
matvecdata = ([],[],[])
matdata = ([],[],[])
vecdata = ([],[])
test_single_field_fe_space(V0,matvecdata,matdata,vecdata)

cell_map = get_cell_map(V0)

f(x) = sin(4*pi*(x[1]-x[2]^2))+1

fh = interpolate_everywhere(f, V0)

fh = interpolate_dirichlet(f, V0)

dirichlet_values = compute_dirichlet_values_for_tags(V0,[1,2])

r=[2,2,1,1,2,2,1,1,2,2,2,2,1,2,1,1,1,1,2,2,2,2,1,2,1,1,1,1,2,2,1,1,2,2,1,2,1,1,2,2,1,2,1,1,2,2,1,2,1,1]

@test dirichlet_values == r

dirichlet_values = compute_dirichlet_values_for_tags(V0,2)
@test dirichlet_values == fill(2,num_dirichlet_dofs(V0))

free_values = zero_free_values(V0)
uh = FEFunction(V0,free_values,dirichlet_values)

#using Gridap.Visualization
#
#trian = get_triangulation(model)
#
#writevtk(model,"model")
#
#writevtk(trian,"trian",cellfields=["fh" => fh, "uh" => uh])

#
#
#strian = SkeletonTriangulation(model)
#
#fh_Γ = restrict(get_array(fh),strian)
#
#writevtk(strian,"strian",nsubcells=10,cellfields=["jump_fh" => jump(fh_Γ)])

end # module


text:
include("experimentSummaryInterface.jl")
include("GPInformationInterface.jl")
include("toolsInterface.jl")
include("evaluatorInterface.jl")
include("generatorInterface.jl")
include("selectorInterface.jl")
include("crossoverInterface.jl")
include("mutationInterface.jl")
include("replacementInterface.jl")
include("stopConditionsInterface.jl")

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
include("basic_binary_search_tree.jl")
mutable struct Splay{T}<:AbstractBinarySearchTree_arr where T
    n::Int
    root::Int
    par::Vector{Int}
    lch::Vector{Int}
    rch::Vector{Int}
    val::Vector{T}

    cnt::Vector{Int}
    size::Vector{Int}
end

max_stars_repo_name: KohRongSoon/Julia

max_issues_repo_name: KohRongSoon/Julia

max_forks_repo_name: KohRongSoon/Julia


text:
#
# Query input data.
#

typealias InputContext Dict{Symbol,Any}
typealias InputFrame Nullable{AbstractVector{Int}}
typealias InputSlotFlow Pair{Symbol, OutputFlow}
typealias InputSlotFlows Vector{InputSlotFlow}

immutable InputFlow <: AbstractVector{Any}
    ctx::InputContext
    ity::Input
    len::Int
    vals::AbstractVector
    frameoffs::InputFrame
    slotflows::InputSlotFlows

    function InputFlow(
            ctx,
            dom,
            vals::AbstractVector,
            frameoffs::InputFrame=InputFrame(),
            slotflows::InputSlotFlows=InputSlotFlow[])
        relative = !isnull(frameoffs)
        slots = InputSlot[n => output(pflow) for (n, pflow) in slotflows]
        ity = Input(dom) |> setrelative(relative) |> setslots(slots)
        len = length(vals)
        if !isnull(frameoffs)
            @assert length(get(frameoffs)) >= 1
            @assert get(frameoffs)[1] == 1
            @assert get(frameoffs)[end] == len+1
        end
        for (n, pflow) in slotflows
            @assert length(pflow) == len
        end
        return new(ctx, ity, len, vals, frameoffs, slotflows)
    end
end

InputFlow(ctx, dom, vals::AbstractVector, slotflows::InputSlotFlows) =
    InputFlow(ctx, dom, vals, InputFrame(), slotflows)

# Array interface.

size(flow::InputFlow) = (flow.len,)
length(flow::InputFlow) = flow.len

function getindex(flow::InputFlow, i::Int)
    if isnull(flow.frameoffs) && isempty(flow.slotflows)
        return flow.vals[i]
    end
    val =
        if isnull(flow.frameoffs)
            flow.vals[i]
        else
            frameoffs = get(flow.frameoffs)
            frame = 0
            while frameoffs[frame+1] <= i
                frame += 1
            end
            l = frameoffs[frame]
            r = frameoffs[frame+1]
            (i - l + 1, view(flow.vals, l:r-1))
        end
    if isempty(flow.slotflows)
        return val
    end
    return (val, ((n => pflow[i]) for (n, pflow) in flow.slotflows)...)
end

Base.array_eltype_show_how(::InputFlow) = (true, "")
summary(flow::InputFlow) = "InputFlow[$(flow.len) \ud7 $(flow.ity)]"

# Components and other properties.

context(flow::InputFlow) = flow.ctx

input(flow::InputFlow) = flow.ity
domain(flow::InputFlow) = domain(flow.ity)
isfree(flow::InputFlow) = isfree(flow.ity)
isrelative(flow::InputFlow) = isrelative(flow.ity)
slots(flow::InputFlow) = slots(flow.ity)

values(flow::InputFlow) = flow.vals
frameoffsets(flow::InputFlow) = get(flow.frameoffs)
slotflows(flow::InputFlow) = flow.slotflows

# Converting to a smaller input itynature.

function narrow(flow::InputFlow, ity::Input)
    if mode(flow.ity) == mode(ity)
        return flow
    end
    @assert fits(mode(flow.ity), mode(ity)) "fits($(flow.ity), $(ity))"
    frameoffs = !isrelative(ity) ? InputFrame() : flow.frameoffs
    slotflows =
        if length(flow.slotflows) == length(slots(ity))
            flow.slotflows
        elseif length(slots(ity)) == 0
            InputSlotFlow[]
        else
            keep = Set(n for (n, sl) in slots(ity))
            filter(p -> p.first in keep, flow.slotflows)
        end
    return InputFlow(flow.ctx, domain(ity), flow.vals, frameoffs, slotflows)
end

# Converting output flow to a new input flow.

function distribute(iflow::InputFlow, oflow::OutputFlow)
    @assert length(iflow) == length(oflow)
    return InputFlow(
        iflow.ctx,
        domain(oflow),
        values(oflow),
        !isnull(iflow.frameoffs) ?
            InputFrame(dist_frameoffs(get(iflow.frameoffs), offsets(oflow))) :
            iflow.frameoffs,
        !isempty(iflow.slotflows) ?
            dist_slotflows(iflow.slotflows, offsets(oflow)) :
            iflow.slotflows)
end

function dist_frameoffs(frameoffs::AbstractVector{Int}, offs::AbstractVector{Int})
    len = 1
    for i = 2:endof(frameoffs)
        l = offs[frameoffs[i-1]]
        r = offs[frameoffs[i]]
        if l < r
            len += 1
        end
    end
    distoffs = Vector{Int}(len)
    distoffs[1] = 1
    k = 1
    for i = 2:endof(frameoffs)
        l = offs[frameoffs[i-1]]
        r = offs[frameoffs[i]]
        if l < r
            k += 1
            distoffs[k] = r
        end
    end
    return distoffs
end

dist_frameoffs(frameoffs::AbstractVector{Int}, ::OneTo) = frameoffs

function dist_slotflows(slotflows::InputSlotFlows, offs::AbstractVector{Int})
    remap = Vector{Int}(offs[end]-1)
    for i = 1:endof(offs)-1
        l = offs[i]
        r = offs[i+1]
        for k = l:r-1
            remap[k] = i
        end
    end
    return InputSlotFlow[n => pflow[remap] for (n, pflow) in slotflows]
end

dist_slotflows(slotflows::InputSlotFlows, ::OneTo) = slotflows

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
@enum(BranchCost, BC_INFEASIBLE, BC_INTERVAL, BC_INTERVAL_REV, BC_INTERVAL_LP, BC_INTERVAL_LP_REV)

Base.@kwdef mutable struct BranchCostStorage{T<:Real}
    cost::BranchCost         = BC_INTERVAL
    𝛹n::Vector{T}           = T[]
    𝛹p::Vector{T}           = T[]
    δn::Vector{T}            = T[]
    δp::Vector{T}            = T[]
    ηn::Vector{T}            = T[]
    ηp::Vector{T}            = T[]
    μ1::T                    = 0.1
    μ2::T                    = 1.3
    μ3::T                    = 0.8
    β::T                     = 0.05
    μ_score::T               = 0.15
end
function initialize!(d::BranchCostStorage{T}, n::Int) where T <:AbstractFloat
    append!(d.𝛹n, ones(T,n));  append!(d.𝛹p, ones(T,n))
    append!(d.δn, zeros(T,n));  append!(d.δp, zeros(T,n))
    append!(d.ηn, zeros(T,n));  append!(d.ηp, zeros(T,n))
    return
end

#=
LP          -> COPY TO RELAXED SOLVER AND SOLVE
MILP        -> COPY TO RELAXED SOLVER AND SOLVE
SOCP        -> COPY TO RELAXED SOLVER AND SOLVE
MISOCP      -> COPY TO RELAXED SOLVER AND SOLVE
DIFF_CVX    -> COPY TO NLP SOLVER AND SOLVE (POTENTIAL MULTISTART)
NS_CVX      -> COPY TO NLP SOLVER AND SOLVE (POTENTIAL MULTISTART)
DIFF_NCVX   -> APPLY GLOBAL SOLVER (UNLESS USER REQUEST LOCAL SOLVE THEN NLP)
NS_NCVX     -> APPLY GLOBAL SOLVER (UNLESS USER REQUEST LOCAL SOLVE THEN NLP)
MINCVX      -> APPLY GLOBAL SOLVER (LOCAL SOLVE OPTION FUTURE FEATURE)
=#

abstract type AbstractProblemType end
struct LP <: AbstractProblemType end
struct MILP <: AbstractProblemType end
struct SOCP <: AbstractProblemType end
struct MISOCP <: AbstractProblemType end
struct DIFF_CVX <: AbstractProblemType end
struct MINCVX <: AbstractProblemType end

const ANY_PROBLEM_TYPE = Union{Nothing, LP, MILP, SOCP, MISOCP, DIFF_CVX, MINCVX}

@enum(GlobalEndState, GS_OPTIMAL, GS_INFEASIBLE, GS_NODE_LIMIT,
                      GS_ITERATION_LIMIT, GS_RELATIVE_TOL,
                      GS_ABSOLUTE_TOL, GS_TIME_LIMIT, GS_UNSET)

export EAGOParameters
"""
$(TYPEDEF)

Storage for parameters that do not change during a global solve.

$(TYPEDFIELDS)
"""
Base.@kwdef mutable struct EAGOParameters

    # Presolving options
    "Should EAGO attempt to remove type assert issues for user-defined functions (default = false)"
    presolve_scrubber_flag::Bool = false
    "Create and use DAG representations of user-defined function (default = false)."
    presolve_to_JuMP_flag::Bool = false
    "Rerranges the DAG using registered transformations (default = false)"
    presolve_flatten_flag::Bool = false

    # Conic reformulations
    "Attempt to bridge convex constraint to second order cone"
    conic_convert_quadratic::Bool = false

    # Iteration logging options
    "Turns logging on records global bounds, node count and run time. Additional
     options are available for recording information specific to subproblems (default = false)."
    log_on::Bool = false
    "Turns on logging of times and feasibility of subproblems (default = false)"
    log_subproblem_info::Bool = false
    "Log data every `log_interval` iterations (default = 1)."
    log_interval::Int = 1

    # Optimizer display options
    "The amount of information that should be printed to console while solving
    values range from 0 - 4: 0 is silent, 1 shows iteration summary statistics
    only, 2-4 show varying degrees of details about calculations within each
    iteration (default = 1)."
    verbosity::Int = 1
    "Display summary of iteration to console every `output_iterations` (default = 1000)"
    output_iterations::Int = 1000
    "Display header for summary to console every `output_iterations` (default = 10000)"
    header_iterations::Int = 100000

    # Node branching options
    "Convex coefficient used to select branch point. Branch point is given by
    `branch_cvx_factor*xmid + (1-branch_cvx_factor)*xsol` (default = 0.25)"
    branch_cvx_factor::Float64 = 0.25
    "Minimum distance from bound to have branch point normalized by width of
    dimension to branch on (default = 0.2)"
    branch_offset::Float64 = 0.15
    "Indicates that pseudocost branching should be used"
    branch_pseudocost_on::Bool = false
    "Variables to branch on (default is all nonlinear)."
    branch_variable::Vector{Bool} = Bool[]
    "[FUTURE FEATURE, NOT CURRENTLY IMPLEMENTED] Number of times repeat node
    processing priorto branching (default = 4)."
    branch_max_repetitions::Int = 4
    "[FUTURE FEATURE, NOT CURRENTLY IMPLEMENTED] Volume ratio tolerance required
    to repeat processing the current node (default = 0.9)"
    branch_repetition_tol::Float64 = 0.9

    # Termination limits
    "Maximum number of nodes (default = 1E-7)"
    node_limit::Int = 1*10^7
    "Maximum CPU time in seconds (default = 3600)"
    time_limit::Float64 = 3600.0
    "Maximum number of iterations (default 3E6)"
    iteration_limit::Int = 1E9
    "Absolute tolerance for termination (default = 1E-3)"
    absolute_tolerance::Float64 = 1E-3
    "Relative tolerance for termination (default = 1E-3)"
    relative_tolerance::Float64 = 1E-3
    "Absolute constraint feasibility tolerance"
    absolute_constraint_feas_tolerance::Float64 = 1E-8

    # Options for constraint propagation
    "Depth in B&B tree above which constraint propagation should be disabled (default = 1000)"
    cp_depth::Int = 0
    "Number of times to repeat forward-reverse pass routine (default = 3)"
    cp_repetitions::Int = 0
    "Disable constraint propagation if the ratio of new node volume to beginning node volume exceeds
    this number (default = 0.99)"
    cp_tolerance::Float64 = 0.99
    "Use only valid interval bounds during constraint propagation (default = false)"
    cp_interval_only::Bool = false

    # obbt options
    "Depth in B&B tree above which OBBT should be disabled (default = 6)"
    obbt_depth::Int = 6
    "Number of repetitions of OBBT to perform in preprocessing (default = 3)"
    obbt_repetitions::Int = 3
    "Turn aggresive OBBT on (default = false)"
    obbt_aggressive_on::Bool = true
    "Maximum iteration to perform aggresive OBBT (default = 2)"
    obbt_aggressive_max_iteration::Int = 2
    "Minimum dimension to perform aggresive OBBT (default = 2)"
    obbt_aggressive_min_dimension::Int = 2
    "Tolerance to consider bounds equal (default = 1E-10)"
    obbt_tolerance::Float64 = 1E-10

    # Options for linear bound tightening
    "Depth in B&B tree above which linear FBBT should be disabled (default = 1000)"
    fbbt_lp_depth::Int  = 1000
    "Number of repetitions of linear FBBT to perform in preprocessing (default = 3)"
    fbbt_lp_repetitions::Int  = 3

    # Duality-based bound tightening (DBBT) options
    "Depth in B&B tree above which duality-based bound tightening should be disabled (default = 1E10)"
    dbbt_depth::Int = 10^10
    "New bound is considered equal to the prior bound if within dbbt_tolerance (default = 1E-9)."
    dbbt_tolerance::Float64 = 1E-8

    # Subgradient tightening flag
    "Relax Tag used to specify type of McCormick operator"
    relax_tag::RelaxTag = NS()
    "Perform tightening of interval bounds using subgradients at each factor in
    each nonlinear tape during a forward pass (default = true)."
    subgrad_tighten::Bool = true
    "Perform tightening of interval bounds using subgradients at each factor in
    each nonlinear tape during a reverse pass (default = false)."
    reverse_subgrad_tighten::Bool = false
    "Outer round computed subgradient bounds by this amount"
    subgrad_tol::Float64 = 1E-10
    "Selects the type of relaxation to use for the bilinear term (multiplication): 0 corresponds to
     a standard McCormick arithmetic approach. Options 1-3 augment the standard McCormick relaxation
     with implied apriori relaxations, (1) results a subgradient based apriori relaxation approach, (2)
     results in a affine arithmetic based apriori approach, and (3) corresponds to a enumerative apriori
     relaxation based approach"
    mul_relax_style::Int = 0

    # Tolerance to add cuts and max number of cuts
    "Minimum number of cuts at each node to attempt (unsafe cuts not necessarily added)"
    cut_min_iterations::Int = 2
    "Maximum number of cuts at each node to attempt"
    cut_max_iterations::Int = 8
    "Absolute tolerance checked for continuing cut"
    cut_tolerance_abs::Float64 = 1E-6
    "Relative tolerance checked for continuing cut"
    cut_tolerance_rel::Float64 = 1E-3

    "Use tolerances to determine safe cuts in a Khajavirad 2018 manner"
    cut_safe_on::Bool = true
    "Lower tolerance for safe-lp cut, Khajavirad 2018"
    cut_safe_l::Float64 = 1E-7
    "Upper tolerance for safe-lp cut, Khajavirad 2018"
    cut_safe_u::Float64 = 1E7
    "Constant tolerance for safe-lp cut, Khajavirad 2018"
    cut_safe_b::Float64 = 1E9

    "Solve upper problem for every node with depth less than `upper_bounding_depth`
    and with a probabilityof (1/2)^(depth-upper_bounding_depth) otherwise (default = 8)"
    upper_bounding_depth::Int = 8

    # handling for domain violations
    "Amount about a domain violation to ignore when propagating bounds."
    domain_violation_guard_on::Bool = false
    "Amount about a domain violation to ignore when propagating bounds."
    domain_violation_ϵ::Float64 = 1E-9

    "If true, then EAGO forgos its default configuration process for subsolvers"
    user_solver_config::Bool = false

    "Absolute tolerance used to check for integrality of decision variables"
    integer_abs_tol::Float64 = 1E-9
    "Relative tolerance used to check for integrality of decision variables"
    integer_rel_tol::Float64 = 1E-9
end
const EAGO_PARAMETERS = fieldnames(EAGOParameters)

"""
$(TYPEDEF)

A structure used to hold objectives and constraints added to EAGO model.
The constraints generally aren't used for relaxations.
"""
Base.@kwdef mutable struct InputProblem

    # counts for number of variables and constraints
    _variable_count::Int = 0
    _constraint_count::Int = 0

    # constraint index to function and set dictionaries
    _vi_leq_constraints::Dict{CI{VI,LT}, Tuple{VI,LT}} = Dict{CI{VI,LT}, Tuple{VI,LT}}()
    _vi_geq_constraints::Dict{CI{VI,GT}, Tuple{VI,GT}} = Dict{CI{VI,GT}, Tuple{VI,GT}}()
    _vi_eq_constraints::Dict{CI{VI,ET}, Tuple{VI,ET}} = Dict{CI{VI,ET}, Tuple{VI,ET}}()
    _vi_it_constraints::Dict{CI{VI,IT}, Tuple{VI,IT}} = Dict{CI{VI,IT}, Tuple{VI,IT}}()
    _vi_zo_constraints::Dict{CI{VI,ZO}, Tuple{VI,ZO}} = Dict{CI{VI,ZO}, Tuple{VI,ZO}}()
    _vi_int_constraints::Dict{CI{VI,MOI.Integer}, Tuple{VI,MOI.Integer}} = Dict{CI{VI,MOI.Integer}, Tuple{VI,MOI.Integer}}()

    _linear_leq_constraints::Dict{CI{SAF,LT}, Tuple{SAF,LT}} = Dict{CI{SAF,LT}, Tuple{SAF,LT}}()
    _linear_geq_constraints::Dict{CI{SAF,GT}, Tuple{SAF,GT}} = Dict{CI{SAF,GT}, Tuple{SAF,GT}}()
    _linear_eq_constraints::Dict{CI{SAF,ET}, Tuple{SAF,ET}} = Dict{CI{SAF,ET}, Tuple{SAF,ET}}()

    _quadratic_leq_constraints::Dict{CI{SQF,LT}, Tuple{SQF,LT}} = Dict{CI{SQF,LT}, Tuple{SQF,LT}}()
    _quadratic_geq_constraints::Dict{CI{SQF,GT}, Tuple{SQF,GT}} = Dict{CI{SQF,GT}, Tuple{SQF,GT}}()
    _quadratic_eq_constraints::Dict{CI{SQF,ET}, Tuple{SQF,ET}} = Dict{CI{SQF,ET}, Tuple{SQF,ET}}()

    _conic_second_order::Dict{CI{VECOFVAR,SOC}, Tuple{VECOFVAR,SOC}} = Dict{CI{VECOFVAR,SOC}, Tuple{VECOFVAR,SOC}}()

    # constraint index to primal storage dictionary
    _linear_leq_primal::Dict{CI{SAF,LT},Float64} = Dict{CI{SAF,LT},Float64}()
    _linear_geq_primal::Dict{CI{SAF,GT},Float64} = Dict{CI{SAF,GT},Float64}()
    _linear_eq_primal::Dict{CI{SAF,ET},Float64} = Dict{CI{SAF,ET},Float64}()

    _quadratic_leq_primal::Dict{CI{SQF,LT},Float64} = Dict{CI{SQF,LT},Float64}()
    _quadratic_geq_primal::Dict{CI{SQF,GT},Float64} = Dict{CI{SQF,GT},Float64}()
    _quadratic_eq_primal::Dict{CI{SQF,ET},Float64} = Dict{CI{SQF,ET},Float64}()

    _linear_leq_prob_to_ip::Dict{CI{SAF,LT},CI{SAF,LT}} = Dict{CI{SAF,LT},CI{SAF,LT}}()
    _linear_geq_prob_to_ip::Dict{CI{SAF,GT},CI{SAF,GT}} = Dict{CI{SAF,GT},CI{SAF,GT}}()
    _linear_eq_prob_to_ip::Dict{CI{SAF,ET},CI{SAF,ET}} = Dict{CI{SAF,ET},CI{SAF,ET}}()

    _quadratic_leq_prob_to_ip::Dict{CI{SQF,LT},CI{SQF,LT}} = Dict{CI{SQF,LT},CI{SQF,LT}}()
    _quadratic_geq_prob_to_ip::Dict{CI{SQF,GT},CI{SQF,GT}} = Dict{CI{SQF,GT},CI{SQF,GT}}()
    _quadratic_eq_prob_to_ip::Dict{CI{SQF,ET},CI{SQF,ET}} = Dict{CI{SQF,ET},CI{SQF,ET}}()

    # nonlinear constraint storage
    _objective::Union{VI,SAF,SQF,Nothing} = nothing

    # nlp constraints (set by MOI.set(m, ::NLPBlockData...) in optimizer.jl)
    _nlp_data::Union{MOI.NLPBlockData,Nothing} = nothing

    # objective sense information (set by MOI.set(m, ::ObjectiveSense...) in optimizer.jl)
    _optimization_sense::MOI.OptimizationSense = MOI.MIN_SENSE
end

function MOI.empty!(ip::InputProblem)

    # empties all arrays and dictionaries in input problem
    for field in fieldnames(InputProblem)
        field_value = getfield(ip, field)
        if (field_value isa Array) || (field_value isa Dict)
            empty!(field_value)
        end
    end

    # resets all non-array and dictionaries in input problem to original value
    ip._variable_count = 0
    ip._constraint_count = 0
    ip._objective = nothing
    ip._nlp_data = nothing
    ip._optimization_sense = MOI.MIN_SENSE
    return
end

function Base.isempty(x::InputProblem)

    is_empty_flag = true
    new_input_problem = InputProblem()

    for field in fieldnames(InputProblem)
        field_value = getfield(x, field)
        if (field_value isa Array) || (field_value isa Dict)
            if !isempty(field_value)
                is_empty_flag = false
                break
            end
        elseif field_value isa Number
            if getfield(new_input_problem, field) != field_value
                is_empty_flag = false
                break
            end
        end
    end
    is_empty_flag &= isnothing(x._nlp_data)
    is_empty_flag &= isnothing(x._objective)
    return is_empty_flag
end

# Helper functions which simplify
_constraints(m::InputProblem, ::Type{VI}, ::Type{LT}) = m._vi_leq_constraints
_constraints(m::InputProblem, ::Type{VI}, ::Type{GT}) = m._vi_geq_constraints
_constraints(m::InputProblem, ::Type{VI}, ::Type{ET}) = m._vi_eq_constraints
_constraints(m::InputProblem, ::Type{VI}, ::Type{IT}) = m._vi_it_constraints
_constraints(m::InputProblem, ::Type{VI}, ::Type{ZO}) = m._vi_zo_constraints
_constraints(m::InputProblem, ::Type{VI}, ::Type{MOI.Integer}) = m._vi_int_constraints

_constraints(m::InputProblem, ::Type{SAF}, ::Type{LT}) = m._linear_leq_constraints
_constraints(m::InputProblem, ::Type{SAF}, ::Type{GT}) = m._linear_geq_constraints
_constraints(m::InputProblem, ::Type{SAF}, ::Type{ET}) = m._linear_eq_constraints

_constraints(m::InputProblem, ::Type{SQF}, ::Type{LT}) = m._quadratic_leq_constraints
_constraints(m::InputProblem, ::Type{SQF}, ::Type{GT}) = m._quadratic_geq_constraints
_constraints(m::InputProblem, ::Type{SQF}, ::Type{ET}) = m._quadratic_eq_constraints

_constraint_primal(m::InputProblem, ::Type{SAF}, ::Type{LT}) = m._linear_leq_primal
_constraint_primal(m::InputProblem, ::Type{SAF}, ::Type{GT}) = m._linear_geq_primal
_constraint_primal(m::InputProblem, ::Type{SAF}, ::Type{ET}) = m._linear_eq_primal

_constraint_primal(m::InputProblem, ::Type{SQF}, ::Type{LT}) = m._quadratic_leq_primal
_constraint_primal(m::InputProblem, ::Type{SQF}, ::Type{GT}) = m._quadratic_geq_primal
_constraint_primal(m::InputProblem, ::Type{SQF}, ::Type{ET}) = m._quadratic_eq_primal

"""
Extracts primal constraint value from local problem and saves result to appropriate
field of the global optimizer.
"""
function _extract_primal!(d, m::InputProblem, ::Type{F}, ::Type{S}) where {F,S}
    for (k, v) in _constraint_index_to_ip(m, F, S)
        _constraint_primal(m, F, S)[v] = MOI.get(d, MOI.ConstraintPrimal(), k)
    end
    return nothing
end

function _extract_primal_linear!(d, ip::InputProblem)
    _extract_primal!(d, ip, SAF, LT)
    _extract_primal!(d, ip, SAF, GT)
    _extract_primal!(d, ip, SAF, ET)
end
function _extract_primal_quadratic!(d, ip::InputProblem)
    _extract_primal!(d, ip, SQF, LT)
    _extract_primal!(d, ip, SQF, GT)
    _extract_primal!(d, ip, SQF, ET)
end

_constraint_index_to_ip(m::InputProblem, ::Type{SAF}, ::Type{LT}) = m._linear_leq_prob_to_ip
_constraint_index_to_ip(m::InputProblem, ::Type{SAF}, ::Type{GT}) = m._linear_geq_prob_to_ip
_constraint_index_to_ip(m::InputProblem, ::Type{SAF}, ::Type{ET}) = m._linear_eq_prob_to_ip

_constraint_index_to_ip(m::InputProblem, ::Type{SQF}, ::Type{LT}) = m._quadratic_leq_prob_to_ip
_constraint_index_to_ip(m::InputProblem, ::Type{SQF}, ::Type{GT}) = m._quadratic_geq_prob_to_ip
_constraint_index_to_ip(m::InputProblem, ::Type{SQF}, ::Type{ET}) = m._quadratic_eq_prob_to_ip

"""
Adds a constraint to the local problem storing the new constraint index and the associated
index in the input problem.
"""
function _add_constraint_store_ci!(d, m::InputProblem, ::Type{F}, ::Type{S}) where {F,S}
    for (ci_ip, fs) in _constraints(m, F, S)
        ci_wp = MOI.add_constraint(d, fs[1], fs[2])
        _constraint_index_to_ip(m, F, S)[ci_wp] = ci_ip
    end
    return nothing
end

function _add_constraint_store_ci_linear!(d, ip::InputProblem)
    _add_constraint_store_ci!(d, ip, SAF, LT)
    _add_constraint_store_ci!(d, ip, SAF, GT)
    _add_constraint_store_ci!(d, ip, SAF, ET)
end
function _add_constraint_store_ci_quadratic!(d, ip::InputProblem)
    _add_constraint_store_ci!(d, ip, SQF, LT)
    _add_constraint_store_ci!(d, ip, SQF, GT)
    _add_constraint_store_ci!(d, ip, SQF, ET)
end

"""
$(TYPEDEF)

A structure used to expressions and problem descriptions EAGO uses to formulate
relaxed problems.
"""
Base.@kwdef mutable struct ParsedProblem

    # Problem classification (set in parse_classify_problem!)
    _problem_type::ANY_PROBLEM_TYPE = nothing

    "_objective_saf stores the objective and is used for constructing linear affine cuts"
    _objective_saf::SAF = SAF(SAT[], 0.0)
    _objective::Union{VI,AffineFunctionIneq,BufferedQuadraticIneq,BufferedNonlinearFunction,Nothing} = nothing

    # objective sense information (set by convert_to_min in parse.jl)
    _optimization_sense::MOI.OptimizationSense = MOI.MIN_SENSE

    # non-single variable constraints (set in initial_parse)
    _saf_leq::Vector{AffineFunctionIneq} = AffineFunctionIneq[]
    _saf_eq::Vector{AffineFunctionEq} = AffineFunctionEq[]
    _sqf_leq::Vector{BufferedQuadraticIneq} = BufferedQuadraticIneq[]
    _sqf_eq::Vector{BufferedQuadraticEq} = BufferedQuadraticEq[]
    _conic_second_order::Vector{BufferedSOC} = BufferedSOC[]

    # nlp constraints
    _nlp_data::Union{MOI.NLPBlockData,Nothing} = nothing
    _nonlinear_constr::Vector{BufferedNonlinearFunction} = BufferedNonlinearFunction[]
    _relaxed_evaluator::Evaluator = Evaluator()

    # variables (set in initial_parse)
    _variable_info::Vector{VariableInfo{Float64}} = VariableInfo{Float64}[]
    _variable_count::Int = 0
end

function MOI.empty!(x::ParsedProblem)

    for field in fieldnames(ParsedProblem)
        field_value = getfield(x, field)
        if (field_value isa Array) || (field_value isa Dict)
            empty!(field_value)
        end
    end

    x._objective    = nothing
    x._problem_type = nothing
    x._nlp_data     = nothing

    x._optimization_sense = MOI.MIN_SENSE
    x._relaxed_evaluator = Evaluator()
    x._objective_saf = SAF(SAT[], 0.0)
    x._variable_count = 0
    return
end

function Base.isempty(x::ParsedProblem)

    is_empty_flag = true

    new_input_problem = ParsedProblem()
    for field in fieldnames(ParsedProblem)

        field_value = getfield(x, field)
        if field_value isa Array
            if !isempty(field_value)
                is_empty_flag = false
                break
            end

        elseif field_value isa Number
            if getfield(new_input_problem, field) != field_value
                is_empty_flag = false
                break
            end
        end
    end

    is_empty_flag &= isempty(x._objective_saf.terms)
    is_empty_flag &= iszero(x._objective_saf.constant)
    is_empty_flag &= isnothing(x._objective)
    is_empty_flag &= isnothing(x._nlp_data)

    return is_empty_flag
end

"""
    GlobalOptimizer

Optimizer internal to EAGO which performs branch-and-bound in order
to solve nonconvex MINLPs.
"""
Base.@kwdef mutable struct GlobalOptimizer{R,Q,S<:ExtensionType} <: MOI.AbstractOptimizer

    _subsolvers::SubSolvers{R,Q,S} = SubSolvers{R,Q,S}()
    _parameters::EAGOParameters = EAGOParameters()
    _input_problem::InputProblem = InputProblem()
    _working_problem::ParsedProblem = ParsedProblem()
    _auxillary_variable_info::Union{Nothing,_AuxVarData} = nothing

    obbt_variable_values::Vector{Bool} = Bool[]

    enable_optimize_hook::Bool = false
    ext

    # set as user-specified option
    _end_state::GlobalEndState = GS_UNSET
    _termination_status_code::MOI.TerminationStatusCode = MOI.OPTIMIZE_NOT_CALLED
    _result_status_code::MOI.ResultStatusCode = MOI.OTHER_RESULT_STATUS
    _obj_mult::Float64 = 1.0
    _obj_var_slack_added::Bool = false

    _stack::BinaryMinMaxHeap{NodeBB} = BinaryMinMaxHeap{NodeBB}()
    _current_node::NodeBB = NodeBB()

    _first_relax_point_set::Bool = false
    _current_xref::Vector{Float64} = Float64[]
    _candidate_xref::Vector{Float64} = Float64[]

    _use_prior_objective_xref::Bool = false
    _current_objective_xref::Vector{Float64} = Float64[]
    _prior_objective_xref::Vector{Float64} = Float64[]

    _user_branch_variables::Bool = false
    _fixed_variable::Vector{Bool} = Bool[]
    _branch_variable_count::Int = 0
    _branch_to_sol_map::Vector{Int} = Int[]
    _sol_to_branch_map::Vector{Int} = Int[]

    _continuous_solution::Vector{Float64} = Float64[]

    _preprocess_feasibility::Bool = true
    _preprocess_primal_status::MOI.ResultStatusCode = MOI.OTHER_RESULT_STATUS
    _preprocess_dual_status::MOI.ResultStatusCode = MOI.OTHER_RESULT_STATUS
    _preprocess_termination_status::MOI.TerminationStatusCode = MOI.OPTIMIZE_NOT_CALLED

    _lower_primal_status::MOI.ResultStatusCode = MOI.OTHER_RESULT_STATUS
    _lower_dual_status::MOI.ResultStatusCode = MOI.OTHER_RESULT_STATUS
    _lower_termination_status::MOI.TerminationStatusCode = MOI.OPTIMIZE_NOT_CALLED
    _lower_feasibility::Bool = true
    _lower_objective_value::Float64 = -Inf
    _lower_solution::Vector{Float64} = Float64[]
    _lower_lvd::Vector{Float64} = Float64[]
    _lower_uvd::Vector{Float64} = Float64[]

    _last_cut_objective::Float64 = -Inf

    _upper_result_status::MOI.ResultStatusCode = MOI.OTHER_RESULT_STATUS
    _upper_termination_status::MOI.TerminationStatusCode = MOI.OPTIMIZE_NOT_CALLED
    _upper_feasibility::Bool = true
    _upper_objective_value::Float64 = Inf
    _upper_variables::Vector{VI} =  VI[]
    _upper_solution::Vector{Float64} = Float64[]

    _postprocess_feasibility::Bool = true

    # set to time limit in initial_parse! in parse.jl, decremented throughout global_solve in optimize_nonconvex.jl
    _time_left::Float64 = 3600.0

    # set constructor reset on empty! and  to zero in initial parse! in parse.jl
    _start_time::Float64 = 0.0
    _run_time::Float64 = 0.0
    _parse_time::Float64 = 0.0
    _presolve_time::Float64 = 0.0
    _last_preprocess_time::Float64 = 0.0
    _last_lower_problem_time::Float64 = 0.0
    _last_upper_problem_time::Float64 = 0.0
    _last_postprocessing_time::Float64 = 0.0

    # reset in initial_parse! in parse.jl
    _min_converged_value::Float64 = Inf
    _global_lower_bound::Float64 = -Inf
    _global_upper_bound::Float64 = Inf
    _maximum_node_id::Int = 0
    _iteration_count::Int = 0
    _node_count::Int = 0

    # Storage for output, reset in initial_parse! in parse.jl
    _solution_value::Float64 = 0.0
    _feasible_solution_found::Bool = false
    _first_solution_node::Int = -1
    _best_upper_value::Float64 = Inf

    # Optimality-Based Bound Tightening (OBBT) Options
    _obbt_working_lower_index::Vector{Bool} = Bool[]
    _obbt_working_upper_index::Vector{Bool} = Bool[]
    _lower_indx_diff::Vector{Bool} = Bool[]
    _upper_indx_diff::Vector{Bool} = Bool[]
    _old_low_index::Vector{Bool} = Bool[]
    _old_upp_index::Vector{Bool} = Bool[]
    _new_low_index::Vector{Bool} = Bool[]
    _new_upp_index::Vector{Bool} = Bool[]
    _obbt_variables::Vector{VI} = VI[]
    _obbt_variable_count::Int = 0
    _obbt_performed_flag::Bool = false

    # Buffers for fbbt, set in presolve, used in preprocess
    _lower_fbbt_buffer::Vector{Float64} = Float64[]
    _upper_fbbt_buffer::Vector{Float64} = Float64[]
    _cp_improvement::Float64 = 0.0
    _cp_evaluation_reverse::Bool = false

    _cut_iterations::Int = 0
    _cut_add_flag::Bool = false
    _node_repetitions::Int = 0

    _log::Log = Log()

    _affine_relax_ci::Vector{CI{SAF,LT}} = CI{SAF,LT}[]
    _affine_objective_cut_ci::Union{CI{VI,LT},CI{SAF,LT},Nothing} = nothing

    _relaxed_variable_number::Int = 0
    _relaxed_variable_index::Vector{VI} = VI[]
    _relaxed_variable_et::Vector{Tuple{CI{VI, ET}, Int}} = Tuple{CI{VI, ET}, Int}[]
    _relaxed_variable_lt::Vector{Tuple{CI{VI, LT}, Int}} = Tuple{CI{VI, LT}, Int}[]
    _relaxed_variable_gt::Vector{Tuple{CI{VI, GT}, Int}} = Tuple{CI{VI, GT}, Int}[]
    _relaxed_variable_integer::Vector{CI{VI, MOI.Integer}} = CI{VI, MOI.Integer}[]

    _branch_variables::Vector{Bool} = Bool[]
    _nonbranching_int::Bool = false

    _new_eval_constraint::Bool = false
    _new_eval_objective::Bool = false

    _node_to_sv_leq_ci::Dict{Int,CI{VI,LT}} = Dict{Int,CI{VI,LT}}()
    _node_to_sv_geq_ci::Dict{Int,CI{VI,GT}} = Dict{Int,CI{VI,GT}}()
    _nonlinear_evaluator_created::Bool = false

    _branch_cost::BranchCostStorage{Float64} = BranchCostStorage{Float64}()
    _branch_variable_sparsity::SparseMatrixCSC{Bool,Int} = spzeros(Bool,1,1)
    _constraint_infeasiblity::Vector{Float64} = Float64[]
end

const EAGO_OPTIMIZER_ATTRIBUTES = Symbol[:relaxed_optimizer, :upper_optimizer,
                                         :enable_optimize_hook, :ext, :_parameters]
const EAGO_MODEL_STRUCT_ATTRIBUTES = Symbol[:_stack, :_log, :_current_node, :_working_problem, :_input_problem, :_branch_cost]
const EAGO_MODEL_NOT_STRUCT_ATTRIBUTES = setdiff(fieldnames(GlobalOptimizer), union(EAGO_OPTIMIZER_ATTRIBUTES, EAGO_MODEL_STRUCT_ATTRIBUTES))

function MOI.empty!(m::GlobalOptimizer{R,S,Q}) where {R,S,Q}

    # create a new empty optimizer and copy fields to m
    new_optimizer = GlobalOptimizer{R,S,Q}(_subsolvers = m._subsolvers,
                                           _parameters = m._parameters,
                                           _input_problem = m._input_problem,
                                           _working_problem = m._working_problem,
                                           ext = nothing)

    MOI.empty!(new_optimizer._subsolvers)
    MOI.empty!(new_optimizer._input_problem)
    MOI.empty!(new_optimizer._working_problem)
    for field in union(EAGO_MODEL_STRUCT_ATTRIBUTES, EAGO_MODEL_NOT_STRUCT_ATTRIBUTES)
        setfield!(m, field, getfield(new_optimizer, field))
    end

    return nothing
end

function MOI.is_empty(m::GlobalOptimizer{R,S,Q}) where {R,S,Q}

    is_empty_flag = uninitialized(_current_node(m))
    is_empty_flag &= isempty(m._stack)
    is_empty_flag &= isempty(m._log)
    is_empty_flag &= isempty(m._input_problem)
    is_empty_flag &= isempty(m._working_problem)

    new_optimizer = GlobalOptimizer{R,S,Q}(_subsolvers = m._subsolvers,
                                           _parameters = m._parameters,
                                           _input_problem = m._input_problem,
                                           _working_problem = m._working_problem,
                                           ext = nothing)
    for field in EAGO_MODEL_NOT_STRUCT_ATTRIBUTES
        if getfield(m, field) != getfield(new_optimizer, field)
            is_empty_flag = false
            break
        end
    end

    return is_empty_flag
end

MOI.get(m::GlobalOptimizer, ::MOI.ObjectiveBound) = _is_input_min(m) ? m._global_lower_bound : -m._global_upper_bound
MOI.get(m::GlobalOptimizer, ::MOI.ObjectiveValue) = _is_input_min(m) ? m._global_upper_bound : -m._global_lower_bound

_relaxed_optimizer(m::GlobalOptimizer{R,S,Q}) where {R,S,Q} = m._subsolvers.relaxed_optimizer
_upper_optimizer(m::GlobalOptimizer{R,S,Q})   where {R,S,Q} = m._subsolvers.upper_optimizer
_ext(m::GlobalOptimizer{R,S,Q})               where {R,S,Q} = m._subsolvers.ext

@inline function _is_input_min(m::GlobalOptimizer)
    return m._obj_mult == 1.0
end
@inline _branch_cost(m::GlobalOptimizer) = m._branch_cost.cost
@inline _cost_offset_β(m::GlobalOptimizer) = m._branch_cost.β
@inline _branch_cvx_α(m::GlobalOptimizer) =  m._parameters.branch_cvx_factor
@inline _branch_offset_β(m::GlobalOptimizer) = m._parameters.branch_offset
@inline _branch_pseudocost_on(m::GlobalOptimizer) = m._parameters.branch_pseudocost_on
@inline _cut_ϵ_abs(m::GlobalOptimizer) = m._parameters.cut_tolerance_abs
@inline _cut_ϵ_rel(m::GlobalOptimizer) = m._parameters.cut_tolerance_rel
@inline _cut_max_iterations(m::GlobalOptimizer) = m._parameters.cut_max_iterations

@inline _integer_abs_tol(m::GlobalOptimizer) = m._parameters.integer_abs_tol
@inline _integer_rel_tol(m::GlobalOptimizer) = m._parameters.integer_rel_tol

@inline _absolute_tol(m::GlobalOptimizer) = m._parameters.absolute_tolerance
@inline _relative_tol(m::GlobalOptimizer) = m._parameters.relative_tolerance
@inline _constraint_tol(m::GlobalOptimizer) = m._parameters.absolute_constraint_feas_tolerance

@inline _fbbt_lp_depth(m::GlobalOptimizer) = m._parameters.fbbt_lp_depth
@inline _fbbt_lp_repetitions(m::GlobalOptimizer) = m._parameters.fbbt_lp_repetitions

@inline _obbt_depth(m::GlobalOptimizer) = m._parameters.obbt_depth
@inline _obbt_repetitions(m::GlobalOptimizer) = m._parameters.obbt_repetitions
@inline _obbt_tolerance(m::GlobalOptimizer) = m._parameters.obbt_repetitions
@inline _obbt_aggressive_on(m::GlobalOptimizer) = m._parameters.obbt_aggressive_on
@inline _obbt_aggressive_max_iteration(m::GlobalOptimizer) = m._parameters.obbt_aggressive_max_iteration

@inline _user_solver_config(m::GlobalOptimizer) = m._parameters.user_solver_config
@inline _verbosity(m::GlobalOptimizer) = m._parameters.verbosity

@inline _cp_depth(m::GlobalOptimizer) = m._parameters.cp_depth
@inline _cp_repetitions(m::GlobalOptimizer) = m._parameters.cp_repetitions

@inline _iteration_count(m::GlobalOptimizer) = m._iteration_count
@inline _obj_var_slack_added(m::GlobalOptimizer) = m._obj_var_slack_added

@inline _bvi(m::GlobalOptimizer, i::Int) = m._branch_to_sol_map[i]
@inline _svi(m::GlobalOptimizer, i::Int) = m._sol_to_branch_map[i]

@inline _is_branch_var(m::GlobalOptimizer, i) = m._branch_variables[i]

@inline _current_node(m::GlobalOptimizer) = m._current_node
@inline _variable_info(m::GlobalOptimizer, i::Int) = m._input_problem._variable_info[i]
@inline _working_variable_info(m::GlobalOptimizer, i::Int) = m._working_problem._variable_info[i]

@inline _sparsity(::BranchVar, m::GlobalOptimizer, i::Int) = view(m._branch_variable_sparsity,1,:)

@inline _variable_num(::BranchVar, m::GlobalOptimizer) = m._branch_variable_count
@inline _variable_num(::FullVar, m::GlobalOptimizer) = m._working_problem._variable_count

@inline is_integer(::BranchVar, m::GlobalOptimizer, i::Int) = is_integer(_current_node(m), i)
@inline function is_integer(::FullVar, m::GlobalOptimizer, i::Int)
    if _is_branch_var(m,i)
        return is_integer(_current_node(m), _svi(m, i))
    end
    is_integer(_working_variable_info(m,i))
end

@inline _lower_bound(::BranchVar, m::GlobalOptimizer, i::Int) = lower_variable_bounds(_current_node(m), i)
@inline function _lower_bound(::FullVar, m::GlobalOptimizer, i::Int)
     if _is_branch_var(m,i)
         return lower_variable_bounds(_current_node(m), _svi(m, i))
     end
    lower_bound(_working_variable_info(m,i))
end

@inline _upper_bound(::BranchVar, m::GlobalOptimizer, i::Int) = upper_variable_bounds(_current_node(m), i)
@inline function _upper_bound(::FullVar, m::GlobalOptimizer, i::Int)
    if _is_branch_var(m,i)
        return upper_variable_bounds(_current_node(m), _svi(m, i))
    end
    upper_bound(_working_variable_info(m,i))
end

@inline _mid(::BranchVar, m::GlobalOptimizer, i::Int) = mid(_current_node(m), i)
@inline function _mid(::FullVar, m::GlobalOptimizer, i::Int)
    if _is_branch_var(m,i)
        return mid(_current_node(m), _svi(m, i))
    end
    mid(_working_variable_info(m,i))
end

@inline _diam(::BranchVar, m::GlobalOptimizer, i::Int) = diam(_current_node(m), i)
@inline function _diam(::FullVar, m::GlobalOptimizer, i::Int)
    if _is_branch_var(m,i)
        return diam(_current_node(m), _svi(m, i))
    end
    diam(_working_variable_info(m,i))
end

@inline _lower_solution(::BranchVar, m::GlobalOptimizer, i::Int) = m._lower_solution[_bvi(m, i)]
@inline _lower_solution(::FullVar, m::GlobalOptimizer, i::Int) = m._lower_solution[i]

@inline function _set_lower_solution!(::BranchVar, m::GlobalOptimizer, v::Float64, i::Int)
    m._lower_solution[_bvi(m, i)] = v
    return
end
@inline function _set_lower_solution!(::FullVar, m::GlobalOptimizer, v::Float64, i::Int)
    m._lower_solution[i] = v
    return
end

_constraint_primal(m::GlobalOptimizer, ::Type{F}, ::Type{S}) where {F,S} = _constraint_primal(m._input_problem, F, S)

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
# description and problem definition of the POMDP problem
module pomdp
dir_path = @__DIR__
push!(LOAD_PATH, dir_path * "/../sim") # simulation
using sim

@enum CollisionType NO_COLLISION=0 HITTING=1 BEING_HIT=2
@enum DriverType WEAK=1 MEDIUM=2 STRONG=3
@enum DriverAction BRAKE=1 ACC=2 NOTHING=3

const ACTIONS = [BRAKE, ACC, NOTHING]
const DRIVERS = [WEAK, MEDIUM, STRONG]

struct Collision
  d::Float64
  t::Float64
  ctype::Enum

  Collision(d, t, ctype) = new(d, t, ctype)
  Collision() = new(0.0, 0.0, NO_COLLISION)
end

const max_t = 2.0
const max_d = 10.0
function reward(x, u, nx, agent, world)
  c = agent.custom

  base_r = abs(x[3]) > 0.35 * world.road.width ? -1e9 : nx[2]^3
  coll_r = 0.0
  if c.ctype != NO_COLLISION
    t = c.t < 1.0 ? 1.0 : c.t
    coll_r = -1e9 * (max_d - c.d) / 10.0 * (max_t - t)
    coll_r = c.ctype == BEING_HIT ? coll_r * 1e-5 : coll_r
  end

  return base_r + coll_r
end

function P_adv(o::Enum, c::Collision, s::Enum)
  if c.ctype != NO_COLLISION
    P = 0.0
    if s == WEAK
      P = o == BRAKE ? 0.9 : P
      P = o == ACC ? 0.05 : P
      P = o == NOTHING ? 0.05 : P
    elseif s == MEDIUM
      if c.ctype == HITTING
        P = o == BRAKE ? 0.9 : P
        P = o == ACC ? 0.05 : P
        P = o == NOTHING ? 0.05 : P
      elseif c.ctype == BEING_HIT
        P = o == BRAKE ? 0.05 : P
        P = o == ACC ? 0.90 : P
        P = o == NOTHING ? 0.05 : P
      end
    elseif s == STRONG
      P = o == BRAKE ? 0.05 : P
      P = o == ACC ? 0.05 : P
      P = o == NOTHING ? 0.9 : P
    end

    return P
  else
    if o == BRAKE || o == ACC
      return 0.0
    elseif o == NOTHING
      return 1.0
    end
  end
end

function update_belief(b::AbstractArray{Float64, 1}, o::Enum, c::Collision)
  bp = similar(b)
  bp_sum = 0.0
  for i in 1:length(b)
    #println("i -> $(i), b[i] -> $(b[i]), P -> $(P_adv(o, c, DRIVERS[i]))")
    bp[i] = b[i] * P_adv(o, c, DRIVERS[i])
    bp_sum += bp[i]
  end
  bp ./= sum(bp)

  return bp_sum != 0.0 ? bp : b
end

function sample_adv_a(dt::Enum, c::Collision)
  if c.ctype != NO_COLLISION
    # set cP
    cP = Float64[]
    if dt == WEAK
      cP = [0.9, 0.95, 1.0]
    elseif dt == MEDIUM
      if c.ctype == HITTING
        cP = [0.9, 0.95, 1.0]
      elseif c.ctype == BEING_HIT
        cP = [0.05, 0.95, 1.0]
        #cP = [0.05, 0.1, 1.0]
      end
    elseif dt == STRONG
      cP = [0.05, 0.10, 1.0]
    end

    # sample
    r = rand()
    for i in 1:length(ACTIONS)
      if cP[i] >= r
        return ACTIONS[i]
      end
    end
  else
    return NOTHING
  end
end

function adv_controller!(u::AbstractArray{Float64, 1},
                         x::AbstractArray{Float64, 1},
                         dx::AbstractArray{Float64, 1},
                         agent_world::Pair{Agent, World},
                         t::Float64)
  agent = agent_world.first

  v_track = agent.custom[1]
  p_track = agent.custom[2]
  o = agent.custom[3]
  ctrl_d = agent.custom[4]

  u[1] = 0.1 * (v_track - x[2])
  u[2] = 1.0 * (p_track - x[3])
  if o == BRAKE
    u[1] = sign(x[2]) * ctrl_d.mi[1]
  elseif o == ACC
    u[1] = sign(x[2]) * ctrl_d.mx[1]
  elseif o == NOTHING
  end
end

end


text:
using Docile
using Compat
import Base: minimum, maximum, mean, var, std, count, sum, hist
import Base: haskey, keys, values, show

include("varstats.jl")
include("groupby.jl")


text:
# tests for segment selection and StatsModels @formula+df interface
using Random
@testset "segment/model selection" begin

datapath = joinpath(dirname(@__FILE__), "data")

# NOTE: 1. we use intercept=false because of a StatsModels issue where predict with drop_intercept model didn't work
# NOTE: 2. we skip the poisson because it does not converge without an intercept
@testset "$family" for (family, dist, link) in (("gaussian", Normal(), IdentityLink()), ("binomial", Binomial(), LogitLink()))
    data = CSV.read(joinpath(datapath,"gamlr.$family.data.csv"); header=[:y, :x1, :x2, :x3])
    offset = fill(0.001,length(data.y))

    @testset "$L" for L in [LassoModel, GammaLassoModel]
        R = Lasso.pathtype(L)

        @testset "$f" for (f,intercept) in ((@formula(y~x1+x2+x3), true), (@formula(y~0+x1+x2+x3), false))
            path = fit(R, f, data, dist, link; intercept=intercept, offset=offset)

            @testset "$(typeof(select))" for select in [MinAIC(), MinAICc(), MinBIC(), MinCVmse(path), MinCV1se(path)]
                Random.seed!(421)
                m = fit(L, f, data, dist, link; select=select, intercept=intercept, offset=offset)

                Random.seed!(421)
                pathcoefs = coef(path, select)

                Random.seed!(421)
                pathpredict = Lasso.predict(path, data; select=select, offset=offset)

                @test pathcoefs == coef(m)
                if isa(dist, Normal)
                    @test pathpredict == GLM.predict(m, data) + offset
                else
                    @test pathpredict == GLM.predict(m, data; offset=offset)
                end
            end
        end
    end
end

end


text:
using Turing, ContinuousBenchmarks

include(joinpath(ContinuousBenchmarks.STAN_DATA_DIR, "MoC-stan.data.jl"))
include(joinpath(ContinuousBenchmarks.STAN_MODELS_DIR, "MoC.model.jl"))
#include(joinpath(ContinuousBenchmarks.BENCH_DIR, "MoC-stan.run.jl"))

setadbackend(:reverse_diff)

data = map(ind -> nbstandata[1][ind], ["K", "V", "M", "N", "z", "w", "doc", "alpha", "beta"])
@tbenchmark(HMC(20, 0.01, 5), nbmodel, (data...))
bench_res = @tbenchmark(HMC(5000, 0.01, 5), nbmodel, (data...))
# bench_res = @tbenchmark(HMCDA(1000, 0.65, 0.3), nbmodel, (data...))
# bench_res = @tbenchmark(NUTS(2000, 0.65), nbmodel, (data...))
bench_res[4].names = ["phi[1]", "phi[2]", "phi[3]", "phi[4]"]
logd = build_logd("Mixture-of-Categorical", bench_res...)

#logd["stan"] = stan_d
#logd["time_stan"] = nb_time

print_log(logd)
send_log(logd)


text:
""" Convert wiring diagrams to and from syntactic expressions.

Wiring diagrams are a graphical syntax for morphisms in a monoidal category. As
mathematical objects, they are intermediate between the morphisms themselves and
expressions in the textual syntax: a single morphism may correspond to many
wiring diagrams, and a single wiring diagram may correspond to many syntactic
expressions.

It is trivial to convert a morphism expression to a wiring diagram, but not to
go the other way around.
"""
module WiringDiagramExpressions
export to_hom_expr, to_wiring_diagram

using Compat
using LightGraphs

using ...Syntax
using ...Doctrines: id, compose, otimes, munit, mcopy, mmerge, create, delete
using ...Doctrines.Permutations
using ..WiringDiagramCore, ..WiringLayers
import ..WiringLayers: to_wiring_diagram
using ..WiringDiagramAlgorithms: Junction, add_junctions!,
  crossing_minimization_by_sort

# Expression -> Diagram
#######################

""" Convert a morphism expression into a wiring diagram.

The morphism expression should belong to the doctrine of symmetric monoidal
categories, possibly with diagonals and codiagonals. Thus, the doctrines of
cartesian, cocartesian, and biproduct categories are supported.
"""
function to_wiring_diagram(expr::GATExpr)
  functor((Ports, WiringDiagram), expr;
    terms = Dict(
      :Ob => expr -> Ports([first(expr)]),
      :Hom => expr -> singleton_diagram(Box(expr)),
    )
  )
end

# Diagram -> Expression
#######################

""" Convert a wiring diagram into a morphism expression.
"""
function to_hom_expr(Ob::Type, Hom::Type, d::WiringDiagram)
  box_to_expr(v::Int) = to_hom_expr(Ob, Hom, box(d,v))

  # Step 0: Reduction: Add junction nodes to ensure any wiring layer between
  # two boxes is a permutation.
  d = copy(d)
  add_junctions!(d)

  # Special case: no boxes.
  if nboxes(d) == 0
    return hom_expr_between(Ob, d, input_id(d), output_id(d))
  end

  # Step 1. Transitive reduction.
  transitive_reduction!(Ob, d)

  # Step 2. Alternating series- and parallel-reduction.
  n = nboxes(d)
  while true
    parallel_reduction!(Ob, Hom, d)
    series_reduction!(Ob, Hom, d)
    # Repeat until the box count stops decreasing.
    nboxes(d) >= n && break
    n = nboxes(d)
  end

  # Termination: process all remaining boxes (always at least one, and exactly
  # one if there is no creation or deletion).
  if nboxes(d) > 1
    product = otimes(map(box_to_expr, box_ids(d)))
    encapsulate!(d, box_ids(d), discard_boxes=true, value=product)
  end
  v = first(box_ids(d))
  foldl(compose_simplify_id, [
    hom_expr_between(Ob, d, input_id(d), v),
    box_to_expr(v),
    hom_expr_between(Ob, d, v, output_id(d)),
  ])
end
function to_hom_expr(Syntax::Module, d::WiringDiagram)
  to_hom_expr(Syntax.Ob, Syntax.Hom, d)
end

function to_hom_expr(Ob::Type, Hom::Type, box::Box)
  if box.value isa Hom
    box.value
  else
    dom = otimes(ports_to_obs(Ob, input_ports(box)))
    codom = otimes(ports_to_obs(Ob, output_ports(box)))
    Hom(box.value, dom, codom)
  end
end
function to_hom_expr(Ob::Type, Hom::Type, junction::Junction)
  junction_to_expr(Ob, junction)
end

""" All possible parallel reductions of a wiring diagram.
"""
function parallel_reduction!(Ob::Type, Hom::Type, d::WiringDiagram)
  parallel = Vector{Int}[]
  parallel_unsorted = find_parallel(graph(d), source=input_id(d), sink=output_id(d))
  products = map(collect(parallel_unsorted)) do ((source, target), vs)
    exprs = Hom[ to_hom_expr(Ob, Hom, box(d,v)) for v in vs ]
    σ = crossing_minimization_by_sort(d, vs,
      sources = isnothing(source) ? Int[] : [source],
      targets = isnothing(target) ? Int[] : [target])
    push!(parallel, vs[σ])
    otimes(exprs[σ])
  end

  if !isempty(parallel)
    encapsulate!(d, parallel, discard_boxes=true, values=products)
  end
  return d
end

""" All possible series reductions of a wiring diagram.
"""
function series_reduction!(Ob::Type, Hom::Type, d::WiringDiagram)
  box_to_expr(v::Int) = to_hom_expr(Ob, Hom, box(d,v))

  series = find_series(graph(d), source=input_id(d), sink=output_id(d))
  composites = map(series) do vs
    exprs = Hom[ box_to_expr(vs[1]) ]
    for i in 2:length(vs)
      layer_expr = hom_expr_between(Ob, d, vs[i-1], vs[i])
      append!(exprs, [layer_expr, box_to_expr(vs[i])])
    end
    foldl(compose_simplify_id, exprs)
  end

  if !isempty(series)
    encapsulate!(d, series, discard_boxes=true, values=composites)
  end
  return d
end

""" All possible transitive reductions of a wiring diagram.
"""
function transitive_reduction!(Ob::Type, d::WiringDiagram)
  reduced = transitive_reduction!(copy(graph(d)))
  for edge in collect(edges(graph(d)))
    if !has_edge(reduced, edge)
      for wire in wires(d, edge)
        value = port_value(d, wire.source) # =?= port_value(d, wire.target)
        v = add_box!(d, Junction(value, 1, 1))
        add_wire!(d, Wire(wire.source => Port(v,InputPort,1)))
        add_wire!(d, Wire(Port(v,OutputPort,1) => wire.target))
        rem_wire!(d, wire)
      end
    end
  end
  return d
end

""" Morphism expression for wires between two boxes.

Assumes that the wires form a permutation morphism.
"""
function hom_expr_between(Ob::Type, diagram::WiringDiagram, v1::Int, v2::Int)
  layer = wiring_layer_between(diagram, v1, v2)
  inputs = ports_to_obs(Ob, output_ports(diagram, v1))
  outputs = ports_to_obs(Ob, input_ports(diagram, v2))

  σ = to_permutation(layer)
  @assert !isnothing(σ) "Conversion of non-permutation not implemented"
  @assert inputs == outputs[σ]
  permutation_to_expr(σ, inputs)
end

coerce_ob(Ob::Type, x) = x isa Ob ? x : Ob(Ob,x)
ports_to_obs(Ob::Type, ports) = Ob[ coerce_ob(Ob, port) for port in ports ]

function compose_simplify_id(f::GATExpr, g::GATExpr)
  if head(f) == :id; g
  elseif head(g) == :id; f
  else compose(f,g) end
end

""" Convert a wiring layer into a permutation, if it is one.

Otherwise, return nothing.
"""
function to_permutation(layer::WiringLayer)::Union{Nothing,Vector{Int}}
  nin, nout = layer.ninputs, layer.noutputs
  if nin != nout; return nothing end
  σ = zeros(Int, nin)
  for i in 1:nin
    wires = out_wires(layer, i)
    if length(wires) == 1
      σ[i] = last(first(wires))
    else
      return nothing
    end
  end
  return σ
end

""" Convert a junction node to a morphism expression.
"""
function junction_to_expr(Ob::Type, junction::Junction)
  ob = coerce_ob(Ob, junction.value)
  compose_simplify_id(
    mmerge_foldl(ob, junction.ninputs),
    mcopy_foldl(ob, junction.noutputs)
  )
end

# FIXME: These functions belong elsewhere, probably in the standard library.
function mcopy_foldl(A, n::Int)
  @assert n >= 0
  if n > 2
    compose(mcopy(A), otimes(mcopy_foldl(A, n-1), id(A)))
  elseif n == 2
    mcopy(A)
  elseif n == 1
    id(A)
  else # n == 0
    delete(A)
  end
end
function mmerge_foldl(A, n::Int)
  @assert n >= 0
  if n > 2
    compose(otimes(mmerge_foldl(A, n-1), id(A)), mmerge(A))
  elseif n == 2
    mmerge(A)
  elseif n == 1
    id(A)
  else # n == 0
    create(A)
  end
end

# Graph operations
##################

""" Find parallel compositions in a directed graph.
"""
function find_parallel(g::DiGraph; source=nothing, sink=nothing)
  parallel = Dict{Pair{Union{Int,Nothing},Union{Int,Nothing}},Vector{Int}}()
  for v in 1:nv(g)
    if v in (source, sink); continue end
    # Note: The definition in the literature on series-parallel digraphs
    # requires that `length(v_in) == 1 && length(v_out) == 1`.
    # Replacing the conjunction with a disjunction allows more general
    # parallel reductions while still ensuring that no vertex is counted as
    # parallel in more than one place.
    v_in, v_out = inneighbors(g,v), outneighbors(g,v)
    if length(v_in) == 1 || length(v_out) == 1
      for u in (isempty(v_in) ? [nothing] : v_in)
        for w in (isempty(v_out) ? [nothing] : v_out)
          push!(get!(parallel, u => w, Int[]), v)
        end
      end
    end
  end
  filter(pair -> length(last(pair)) > 1, parallel)
end

""" Find series compositions in a directed graph.
"""
function find_series(g::DiGraph; source=nothing, sink=nothing)::Vector{Vector{Int}}
  series_graph = DiGraph(nv(g))
  for edge in edges(g)
    if (length(outneighbors(g,src(edge))) == 1 &&
        length(inneighbors(g,dst(edge))) == 1 &&
        src(edge) != source && dst(edge) != sink)
      add_edge!(series_graph, edge)
    end
  end
  series = Vector{Int}[]
  for component in weakly_connected_components(series_graph)
    if length(component) > 1
      sub, vmap = induced_subgraph(series_graph, component)
      push!(series, Int[ vmap[v] for v in topological_sort_by_dfs(sub) ])
    end
  end
  series
end

""" Transitive reduction of a directed acyclic graph.

Algorithm described at CS.SE: https://cs.stackexchange.com/a/29133
"""
function transitive_reduction!(g::DiGraph)::DiGraph
  for u in 1:nv(g)
    for v in outneighbors(g, u)
      for w in neighborhood(g, v, nv(g))
        if w != v
          rem_edge!(g, u, w)
        end
      end
    end
  end
  return g
end

end


text:
# Moved to "old" versions of run_random_neutral_walk() to the end of this file.  2/9/21.
# Generates walks (in parallel) for circuits in a circult list c_list
# Does not compute a combined evolvability for the multiple walks.
# Instead computes average evolability of the walks.
# This makes sense because the circuits in c_list do not map to the same phenotype
function run_random_neutral_walk( c_list::Vector{Chromosome}, steps::Int64, maxsteps::Int64, nreps::Int64;
    csvfile::String="" )
  df = DataFrame()
  df.complexity = Float64[]
  df.avg_complexity = Float64[]
  df.all_count = Float64[]
  df.robust_count = Float64[]
  df.evolvable_count = Float64[]
  df.evolvability = Float64[]
  df.circuit_dist = Float64[]
  println("size(df): ",size(df))
  p = c_list[1].params
  for c in c_list
    @assert c.params == p
  end
  row_avg_list = pmap( c->avg_random_neutral_walks( c, steps, maxsteps, nreps ), c_list )
  #row_avg_list = map( c->avg_random_neutral_walks( c, steps, maxsteps, nreps ), c_list )
  for row_avg in row_avg_list
    push!(df,row_avg)
  end
  hostname = chomp(open("/etc/hostname") do f read(f,String) end)
  println("# host: ",hostname," with ",nprocs()-1,"  processes: " )
  println("# date: ",Dates.now())
  println("# funcs: ", funcs)
  println("# steps: ",steps)
  println("# maxsteps: ",maxsteps)
  println("# nreps: ",nreps)
  open( csvfile, "w" ) do f
    println(f,"# date : ",Dates.now())
    println(f,"# host: ",hostname," with ",nprocs()-1,"  processes: " )
    println(f,"# funcs: ", funcs)
    print_parameters(f,p,comment=true),
    println(f,"# steps: ",steps)
    println(f,"# maxsteps: ",maxsteps)
    println(f,"# nreps: ",nreps)
    CSV.write( f, df, append=true, writeheader=true )
  end
  df
end

# Helper function for the above run_random_neutral_walk()
function avg_random_neutral_walks( c::Chromosome, steps::Int64, maxsteps::Int64, nreps::Int64 )
    df_row = Tuple[]
    for i = 1:nreps
      push!(df_row, random_neutral_walk( c, steps, maxsteps, evo_count_only=false ) )
      #println("df_row: ",df_row[end])
    end
    df_row_avg = [ mean(df_row[k][j] for k = 1:length(df_row)) for j = 1:length(df_row[end]) ]
    #println("df_row_avg: ",df_row_avg)
  df_row_avg
end

# Returns an evolvable_count or a dataframe row depending on whether evo_count_only is true or false
# Does NOT use mut_evolve().  Instead implements basic version of neutral walk algorithm.
function random_neutral_walk( c::Chromosome, steps::Int64, maxsteps::Int64; evo_count_only::Bool=false )
  orig_complexity =  complexity5(c)
  c_orig = deepcopy(c)
  all_outputs_list = Goal[]
  all_unique_list = Goal[]
  avg_complexity = zeros(Float64,steps)   # mean complexity of c over steps
  all_count = zeros(Int64,steps)
  robust_count = zeros(Int64,steps)
  evolvable_count = zeros(Int64,steps)
  evolvability = zeros(Float64,steps)
  circuit_dist = zeros(Float64,steps)
  evolvable_count = zeros(Int64,steps)
  evolvability = zeros(Float64,steps)
  funcs = default_funcs(c.params.numinputs)
  goal = output_values(c)
  total_steps = 0  # Number of mutation steps
  #println("goal: ",goal,"  complexity: ",complexity5(c))
  for i = 1:steps
    j = 1
    while j <= maxsteps   # terminated by a break statement
      sav_c = deepcopy(c)
      (c,active) = mutate_chromosome!( c, funcs )
      outputs = output_values(c)
      if outputs == goal
        #println("successful step for i= ",i,"  outputs: ",outputs)
        break
      end
      c = sav_c
      j += 1
    end
    if step == maxsteps
      error("Failed to find neutral mutation")
    end
    total_steps += j
    avg_complexity[i] = complexity5(c)
    all_outputs = mutate_all( c, funcs )
    robustness = length(filter(x->x==goal,all_outputs))
    robust_count[i] =  i == 1 ? robustness : robust_count[i-1] + robustness
    all_count[i] =  i == 1 ? length(all_outputs) : all_count[i-1] + length(all_outputs)
    all_unique_list = unique( vcat( all_unique_list, unique( all_outputs)))
    evolvable_count[i] = length(all_unique_list)
    evolvability[i] = evolvable_count[i]/all_count[i]
    circuit_dist[i] = circuit_distance( c, c_orig )
    #println("k: ",k,"  all_count[i]: ",all_count[i],"  robust_count[i]: ",robust_count[i],"  evolvable_count[i]: ",evolvable_count[i] )
  end
  println("total steps: ",total_steps)
  if evo_count_only
    evolvable_count[i]
  else
    ( orig_complexity, avg_complexity[steps], all_count[steps], robust_count[steps], evolvable_count[steps],
      evolvability[steps], circuit_dist[steps] )
  end
end

function binned_circuit_complexities( p::Parameters, n_circuits::Int64, max_bin_size::Int64,
    complexity_lower_bounds::Vector{Float64} )
  funcs = default_funcs(p.numinputs)
  complexity_bins = [ Chromosome[] for k = 1:length(complexity_lower_bounds) ]
  for i = 1:n_circuits
    c = random_chromosome( p, funcs )
    cmplx = complexity5(c)
    cmplx = cmplx >= 0.0 ? cmplx : 0.0
    k = 1
    while k < length(complexity_lower_bounds) && cmplx >= complexity_lower_bounds[k+1]
      k += 1
    end
    #println("cmplx: ",cmplx,"  k: ",k)
    if length(complexity_bins[k]) < max_bin_size
      push!(complexity_bins[k], c )
    end
  end
  complexity_bins
end

# Seems to be superceded by later versions of the same function.
function run_random_neutral_walk( p::Parameters, ngoals::Int64, steps::Int64, maxsteps::Int64 )
  df = DataFrame()
  df.goal = Goal[]
  df.evolvable_count = Int64[]
  df.complexity = Float64[]
  funcs = default_funcs( p.numinputs )
  for gc = 1:ngoals
    g = randgoal( p.numinputs, p.numoutputs )
    (c,step,worse,same,better,output,matched_goals,matched_goals_list) = mut_evolve_repeat( 10, p, [g], funcs, maxsteps)
    res = mut_evolve_repeat( 10, p, [g], funcs, maxsteps)
    if res == nothing
      error("run_random_neutral_walk: mut_evolve_repeat() failed.  Try increasing maxsteps.")
    end
    (c,step,worse,same,better,output,matched_goals,matched_goals_list) = res
    @assert output_values(c) == g
    complexity = complexity5(c)
    ev_count = random_neutral_walk( c, steps, maxsteps, evo_count_only=true )
    push!(df, (g, ev_count, complexity ))
  end
  df
end

# Starts with a goal.  Then evolves a circuit to map to that goal
# See diary10_13.txt for example runs
function run_random_neutral_walk( p::Parameters, g::Goal, steps::Int64, maxsteps::Int64, nreps::Int64=1; skip_steps::Int64=1 )
  funcs = default_funcs(p.numinputs)
  ev_count = 0
  for _ = 1:nreps
    # evolve a chromosome that outputs goal g.
    res = mut_evolve_repeat( 10, p, [g], funcs, maxsteps)
    if res == nothing
      error("run_random_neutral_walk: mut_evolve_repeat() failed.  Try increasing maxsteps.")
    end
    (c,step,worse,same,better,output,matched_goals,matched_goals_list) = res
    #(c,step,worse,same,better,output,matched_goals,matched_goals_list) = mut_evolve_repeat( 10, p, [g], funcs, maxsteps)
    @assert output_values(c) == g
    ev_count += random_neutral_walk( c, steps, maxsteps, skip_steps=skip_steps, evo_count_only=true )
  end
  ev_count/nreps
end


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
using DelimitedFiles
using SparseArrays
Random.seed!(233)

mpi_init()
sp = sprand(10,10,0.3)
SP = mpi_SparseTensor(sp)
SPt = SP'
sess = Session(); init(sess)
Array(run(sess, SP))-Array(run(sess, SPt))'


text:
@defcomp impactaggregation begin
    regions = Index()

    eloss = Variable(index=[time,regions])
    sloss = Variable(index=[time,regions])
    loss = Variable(index=[time,regions])

    income = Parameter(index=[time,regions])

    water = Parameter(index=[time,regions])
    forests = Parameter(index=[time,regions])
    heating = Parameter(index=[time,regions])
    cooling = Parameter(index=[time,regions])
    agcost = Parameter(index=[time,regions])
    drycost = Parameter(index=[time,regions])
    protcost = Parameter(index=[time,regions])
    entercost = Parameter(index=[time,regions])
    hurrdam = Parameter(index=[time,regions])
    extratropicalstormsdam = Parameter(index=[time,regions])
    species = Parameter(index=[time,regions])
    deadcost = Parameter(index=[time,regions])
    morbcost = Parameter(index=[time,regions])
    wetcost = Parameter(index=[time,regions])
    leavecost = Parameter(index=[time,regions])

    # Other economic losses
    eloss_other = Parameter(index=[time,regions])
    # Other non-economic losses
    sloss_other = Parameter(index=[time,regions])

    function run_timestep(p, v, d, t)

        if is_first(t)
            for r in d.regions
                v.eloss[t, r] = 0.0
                v.sloss[t, r] = 0.0
            end
        else
            for r in d.regions
                v.eloss[t, r] = min(
                    0.0 -
                    p.water[t, r] -
                    p.forests[t, r] -
                    p.heating[t, r] -
                    p.cooling[t, r] -
                    p.agcost[t, r] +
                    p.drycost[t, r] +
                    p.protcost[t, r] +
                    p.entercost[t, r] +
                    p.hurrdam[t, r] +
                    p.extratropicalstormsdam[t, r] +
                    p.eloss_other[t,r],
                    p.income[t, r])

                v.sloss[t, r] = 0.0 +
                    p.species[t, r] +
                    p.deadcost[t, r] +
                    p.morbcost[t, r] +
                    p.wetcost[t, r] +
                    p.leavecost[t, r] +
                    p.sloss_other[t,r]

                v.loss[t, r] = (v.eloss[t, r] + v.sloss[t, r]) * 1000000000.0
            end
        end
    end
end


text:
import LinearAlgebra: I, Diagonal
import FFTW: fftfreq, fft, ifft

"""
    compute_RCWA_method(η::Function, μ::Function, L::Float64, N::Int, k::Float64, i::Int)

Compute the fields using Luke's RCWA method.
The output fields are on a grid compatible with the BIE method.

Arguments:
η:: A periodic function returning the metasurface magnetic conductance
μ:: A periodic function returning the metasurface electric resistance
L:: the length of the unit cell (some multiple of the period of η, μ)
N:: the number of discretization points along the unit cell
k:: the wavenumber/frequency of the incident planewave (assumed TM polarization)
i:: the index of the incident planewave in the Fourier basis

Returns:
NamedTuple with fields:
x:: the positions of the grid points
kz₁:: the values of the wavevector in the z component for y<0
kz₂:: the values of the wavevector in the z component for y>0
Ĩ₁:: an array of the incident fields in the Fourier basis for y<0
Õ₁:: an array of the outgoing fields in the Fourier basis for y<0
Ĩ₂:: an array of the incident fields in the Fourier basis for y>0
Õ₂:: an array of the outgoing fields in the Fourier basis for y>0
"""
function compute_RCWA_method(η::Function, μ::Function, L::Float64, N::Int, k::Float64, i::Int)
    # discretization parameters
    dx = L/N
    x⃗ = range(0, step=dx, length=N)
    k⃗ = 2*pi*fftfreq(N, 1/dx)
    β⃗ = @. sqrt(Complex(k^2 - k⃗^2))
    # construct Fourier-basis operators
    β̂ = Diagonal(β⃗)
    η̃ = ifft(fft(Matrix(Diagonal(η.(x⃗))), 2), 1)
    μ̃ = ifft(fft(Matrix(Diagonal(μ.(x⃗))), 2), 1)
    # build GSTC matrices and scattering matrix
    A = [
        -β̂-k*η̃    -β̂-k*η̃;
        -β̂-k*μ̃     β̂+k*μ̃
    ]
    B = [
        -β̂+k*η̃    -β̂+k*η̃;
        -β̂+k*μ̃     β̂-k*μ̃
    ]
    S = A\B
    # solve
    Ĩ₁ = [n == i ? 1 : 0 for n in 1:N]
    Ĩ₂ = zeros(N)
    inc = [Ĩ₁; Ĩ₂]
    scat = S*inc
    Õ₁, Õ₂ = scat[1:N], scat[N+1:2N]
    (x=x⃗, kz₁=β⃗, kz₂=β⃗, Ĩ₁=Ĩ₁, Õ₁=Õ₁, Ĩ₂=Ĩ₂, Õ₂=Õ₂)
end


text:
using DynamicalSystemsBase
using DelayEmbeddings
using Test
using DelimitedFiles

println("\nTesting mdop_embedding.jl...")
@testset "Nichkawde method MDOP" begin

# For comparison reasons using Travis CI we carry out the integration on a UNIX
# OS and save the resulting time series
# # solve Mackey-Glass-Delay Diff.Eq. as in the Paper
# function mackey_glass(du,u,h,p,t)
#   beta,n,gamma,tau = p
#   hist = h(p, t-tau)[1]
#   du[1] = (beta*hist)/(1+hist^n) - gamma * u[1]
# end
# # set parameters
# h(p,t) = 0
# tau_d = 44
# n = 10
# β = 0.2
# γ = 0.1
# δt = 0.5
# p = (β,n,γ,tau_d)
#
# # time span
# tspan = (0.0, 12000.0)
# u0 = [1.0]
#
# prob = DDEProblem(mackey_glass,u0,h,tspan,p; constant_lags=tau_d)
# alg = MethodOfSteps(Tsit5())
# sol = solve(prob,alg; adaptive=false, dt=δt)
#
# s = [u[1] for u in sol.u]
# s = s[4001:end]

s = readdlm(joinpath(tsfolder, "1.csv"))
s = vec(s)
Y = Dataset(s)

theiler = 57

@testset "beta statistic" begin
    ## Test beta_statistic (core algorithm of mdop_embedding)

    taus = 0:100
    β = @inferred DelayEmbeddings.beta_statistic(Y, s, taus, theiler)
    maxi, max_idx = findmax(β)

    @test maxi>4.1
    @test taus[max_idx]>=50

    # # display results as in Fig. 3 of the paper
    # using Plots
    # plot(taus, β, linewidth = 3, label = "1st embedding cycle")
    # plot!(title = "β-statistic for Mackey Glass System as in Fig. 3 in the Paper")
    # xlabel!("τ")
    # ylabel!("β-Statistic")

    # # display results as in Fig. 3 of the paper
    # using PyPlot
    # pygui(true)
    # figure()
    # plot(taus, β, linewidth = 3, label = "1st embedding cycle")
    # scatter(max_idx-1,maxi, c="red")
    # title("β-statistic for Mackey Glass System as in Fig. 3 in the Paper")
    # xlabel("τ")
    # ylabel("β-Statistic")
    # xticks(0:10:100)
    # grid()

    # test different tau range
    taus2 = 1:4:100
    β2 = @inferred beta_statistic(Y, s, taus2, theiler)
    maxi2, max_idx2 = findmax(β2)

    @test maxi2>4.1
    @test taus2[max_idx2]>=40

    # # display results as in Fig. 3 of the paper
    # using Plots
    # plot(taus2, β2, linewidth = 3, label = "1st embedding cycle")
    # plot!(title = "coarse β-statistic for Mackey Glass System as in Fig. 3 in the Paper")
    # xlabel!("τ")
    # ylabel!("β-Statistic")

end

@testset "mdop_embedding univariate" begin

    taus = 0:100
    β = DelayEmbeddings.beta_statistic(Y, s, taus, theiler)
    Y, τ_vals, ts_vals, FNNs, betas = mdop_embedding(s; τs = taus, w = theiler)
    # for different τs
    taus2 = 1:4:100
    Y2, τ_vals2, ts_vals2, FNNs2, betas2 = mdop_embedding(s; τs = taus2, w = theiler)

    @test round.(β, digits=6) == round.(betas[:,1], digits=6)
    @test size(Y,2) == 5
    @test size(Y,2) == size(Y2,2)
    @test sum(findall(x -> x != 1, ts_vals))==0
    @test sum(abs.(diff(τ_vals)) .< 10) == 0

    # # display results as in Fig. 3 of the paper
    # using Plots
    # # Figure as in Fig.3 in the paper
    # plot(taus, betas[:,1], linewidth = 3, label = "embedding cycle 1")
    # plot!([taus[τ_vals[2]+1]],[betas[τ_vals[2]+1,1]], seriestype = :scatter, color="red", label = "")
    # for i = 2:size(betas,2)
    #   plot!(taus, betas[:,i], linewidth = 3, label = "embedding cycle $i")
    #   plot!([taus[τ_vals[i+1]+1]],[betas[τ_vals[i+1]+1,i]], seriestype = :scatter, color="red", label = "")
    # end
    # plot!(title = "β-statistic's for each embedding cycle of Mackey Glass System as in Fig. 3")
    # xlabel!("delay τ")
    # ylabel!("log10 β(τ)")
    #
    # # Figure of coarse grained analysis
    # taus22 = zeros(length(taus2))
    # [taus22[i]=taus2[i] for i = 1: length(taus2)]
    # plot(taus22, betas2[:,1], linewidth = 3, label = "embedding cycle 1")
    # trueind = findall(x -> x == τ_vals2[2], taus22)
    # plot!(taus22[trueind],[betas2[trueind,1]], seriestype = :scatter, color="red", label = "")
    # for i = 2:size(betas2,2)
    #   plot!(taus22, betas2[:,i], linewidth = 3, label = "embedding cycle $i")
    #   trueind = findall(x -> x == τ_vals2[i+1], taus22)
    #   plot!(taus22[trueind],[betas2[trueind,i]], seriestype = :scatter, color="red", label = "")
    # end
    # plot!(title = "β-statistic's for each embedding cycle of Mackey Glass System as in Fig. 3")
    # xlabel!("delay τ")
    # ylabel!("log10 β(τ)")

end

@testset "estimate τ max (Roessler)" begin
    # For comparison reasons using Travis CI we carry out the integration on a UNIX
    # OS and save the resulting time series
    # roe = Systems.roessler([1.0, 0, 0]; a=0.2, b=0.2, c=5.7)
    # sroe = trajectory(roe, 500; dt = 0.05, Ttr = 100.0)
    # writedlm("2.csv", sroe)

    sroe = readdlm(joinpath(tsfolder, "2.csv"))
    tws = 25:32

    τ_m, L = @inferred mdop_maximum_delay(sroe[:, 2], tws)
    @test τ_m == 26
    τ_m, Ls = @inferred mdop_maximum_delay(Dataset(sroe[:, 1:2]), tws)
    @test τ_m == 26

    # # reproduce Fig.2 of the paper
    # tws = 1:2:101
    # τ_m, L = DelayEmbeddings.mdop_maximum_delay(s[:,2]; tw = tws, samplesize=1.0)
    #
    # using Plots
    # twss = zeros(length(tws))
    # [twss[cnt] = i for (cnt,i) in enumerate(tws)]
    # plot(twss,L, label="")
    # xlabel!("time window")
    # ylabel!("L")

    # tw=1:4:200
    # tau_max, LL = DelayEmbeddings.mdop_maximum_delay(sroe; tw=tw)
    #
    # using Plots
    # gui()
    # twss = zeros(length(tw))
    # [twss[cnt] = i for (cnt,i) in enumerate(tw)]
    # plot(twss,LL, label="")
    # xlabel!("time window")
    # ylabel!("L")
end

@testset "mdop_embedding multivariate" begin
    # For comparison reasons using Travis CI we carry out the integration on a UNIX
    # OS and save the resulting time series
    # roe = Systems.roessler([1.0, 0, 0]; a=0.2, b=0.2, c=5.7)
    # sroe = trajectory(roe, 500; dt = 0.05, Ttr = 100.0)
    # writedlm("2.csv", sroe)

    sroe = readdlm(joinpath(tsfolder, "2.csv"))
    tra = Dataset(sroe)
    w1 = estimate_delay(sroe[:,1], "mi_min")
    w2 = estimate_delay(sroe[:,2], "mi_min")

    theiler = w2
    taus = 0:26
    mc = 10

    Y, τ_vals, ts_vals, FNNs, betas =  mdop_embedding(sroe[:,1]; τs = taus, w = theiler, max_num_of_cycles = mc)

    max_idx, ts_number = @inferred DelayEmbeddings.choose_optimal_tau2(betas)
    @test ts_number == 1
    @test taus[max_idx] == τ_vals[2]

    Y2, τ_vals2, ts_vals2, FNNs2, betas2 = mdop_embedding(tra; τs = taus, w = theiler, max_num_of_cycles = mc)
    ttra = standardize(tra)
    b1 = DelayEmbeddings.beta_statistic(Dataset(ttra[:,ts_vals2[1]]), ttra[:,1], taus, theiler)
    b2 = DelayEmbeddings.beta_statistic(Dataset(ttra[:,ts_vals2[1]]), ttra[:,2], taus, theiler)
    b3 = DelayEmbeddings.beta_statistic(Dataset(ttra[:,ts_vals2[1]]), ttra[:,3], taus, theiler)

    @test betas2[1][:,1] == b1
    @test betas2[1][:,2] == b2
    @test betas2[1][:,3] == b3

    @test size(Y2,2) == 3
    @test τ_vals2[1] == τ_vals2[2] == 0
    @test τ_vals2[3] == maximum(taus)

    @test ts_vals2[2] == ts_vals2[3] == 2
    @test ts_vals2[1] == 3

end

end


text:
# This file includes DDESystems

import DifferentialEquations: MethodOfSteps, Tsit5
import UUIDs: uuid4

"""
    @def_dde_system ex

where `ex` is the expression to define to define a new AbstractDDESystem component type. The usage is as follows:
```julia
@def_dde_system mutable struct MyDDESystem{T1,T2,T3,...,TN,OP,RH,RO,ST,IP,OP} <: AbstractDDESystem
    param1::T1 = param1_default                 # optional field
    param2::T2 = param2_default                 # optional field
    param3::T3 = param3_default                 # optional field
        ⋮
    paramN::TN = paramN_default                 # optional field
    constlags::CL = constlags_default           # mandatory field
    depslags::DL = depslags_default             # mandatory field
    righthandside::RH = righthandside_function  # mandatory field
    history::HST = history_function             # mandatory field
    readout::RO = readout_function              # mandatory field
    state::ST = state_default                   # mandatory field
    input::IP = input_defauult                  # mandatory field
    output::OP = output_default                 # mandatory field
end
```
Here, `MyDDESystem` has `N` parameters. `MyDDESystem` is represented by the `righthandside` and `readout` function. `state`, `input` and `output` is the state, input port and output port of `MyDDESystem`.

!!! warning
    `righthandside` must have the signature
    ```julia
    function righthandside(dx, x, u, t, args...; kwargs...)
        dx .= .... # update dx
    end
    ```
    and `readout` must have the signature
    ```julia
    function readout(x, u, t)
        y = ...
        return y
    end
    ```

!!! warning
    New DDE system must be a subtype of `AbstractDDESystem` to function properly.

# Example
```julia
julia> _delay_feedback_system_cache = zeros(1)
1-element Array{Float64,1}:
 0.0

julia> _delay_feedback_system_tau = 1.
1.0

julia> _delay_feedback_system_constlags = [1.]
1-element Array{Float64,1}:
 1.0

julia> _delay_feedback_system_history(cache, u, t) = (cache .= 1.)
_delay_feedback_system_history (generic function with 1 method)

julia> function _delay_feedback_system_rhs(dx, x, h, u, t,
           cache=_delay_feedback_system_cache, τ=_delay_feedback_system_tau)
           h(cache, u, t - τ)  # Update cache
           dx[1] = cache[1] + x[1]
       end
_delay_feedback_system_rhs (generic function with 3 methods)

julia> @def_dde_system mutable struct MyDDESystem{RH, HST, RO, IP, OP} <: AbstractDDESystem
           constlags::Vector{Float64} = _delay_feedback_system_constlags
           depslags::Nothing = nothing
           righthandside::RH = _delay_feedback_system_rhs
           history::HST = _delay_feedback_system_history
           readout::RO = (x, u, t) -> x
           state::Vector{Float64} = rand(1)
           input::IP = nothing
           output::OP = Outport(1)
       end

julia> ds = MyDDESystem();
```
"""
macro def_dde_system(ex)
    checksyntax(ex, :AbstractDDESystem)
    appendcommonex!(ex)
    foreach(nex -> appendex!(ex, nex), [
        :( alg::$ALG_TYPE_SYMBOL = Causal.MethodOfSteps(Causal.Tsit5()) ),
        :( integrator::$INTEGRATOR_TYPE_SYMBOL = Causal.construct_integrator(
            Causal.DDEProblem, input, (righthandside, history), state, t, modelargs, solverargs;
            alg=alg, modelkwargs=(;
            zip(
                (keys(modelkwargs)..., :constant_lags, :dependent_lags),
                (values(modelkwargs)..., constlags, depslags))...
                ),
            solverkwargs=solverkwargs, numtaps=3) )
        ])
    quote
        Base.@kwdef $ex
    end |> esc
end

##### Define DDE system library.

"""
    DDESystem(; constantlags, depslags, righthandside, history, readout, state, input, output)

Construct a generic DDE system
"""
@def_dde_system mutable struct DDESystem{CL, DL, RH, HST, RO, ST, IP, OP} <: AbstractDDESystem
    constlags::CL
    depslags::DL
    righthandside::RH
    history::HST
    readout::RO
    state::ST
    input::IP
    output::OP
end

"""
    DDESystem(; constantlags, depslags, righthandside, history, readout, state, input, output)

Constructs DelayFeedbackSystem
"""
@def_dde_system mutable struct DelayFeedbackSystem{RH, HST, RO, IP, OP} <: AbstractDDESystem
    constlags::Vector{Float64} = Causal._delay_feedback_system_constlags
    depslags::Nothing = nothing
    righthandside::RH = Causal._delay_feedback_system_rhs
    history::HST = Causal._delay_feedback_system_history
    readout::RO = (x, u, t) -> x
    state::Vector{Float64} = rand(1)
    input::IP = nothing
    output::OP = Outport(1)
end

_delay_feedback_system_cache = zeros(1)
_delay_feedback_system_tau = 1.
_delay_feedback_system_constlags = [1.]
_delay_feedback_system_history(cache, u, t) = (cache .= 1.)
function _delay_feedback_system_rhs(dx, x, h, u, t,
    cache=Causal._delay_feedback_system_cache, τ=Causal._delay_feedback_system_tau)
    h(cache, u, t - τ)  # Update cache
    dx[1] = cache[1] + x[1]
end

##### Pretty-printing

show(io::IO, ds::DDESystem) = print(io,
    "DDESystem(righthandside:$(ds.righthandside), readout:$(ds.readout), state:$(ds.state), t:$(ds.t), ",
    "input:$(ds.input), output:$(ds.output))")
show(io::IO, ds::DelayFeedbackSystem) = print(io,
    "DelayFeedbackSystem(state:$(ds.state), t:$(ds.t), input:$(ds.input), output:$(ds.output))")


text:
module QuadOsc

export quadosc

include("SeriesAccelerations.jl")
using .SeriesAccelerations
using QuadGK

@doc raw"""
    quadosc(fn, a, Inf, fnzeros; ...)

Integrate the function `fn(x)` from `a` to `Inf`. The function `fnzeros(n)`
takes an integer `n` and is such that `fn(fnzeros(n)) == 0`. The algorithm
works by integrating between successive zeros, and accelerating the alternating
series.

The argument `pren` is the number of intervals to integrate before applying the
series acceleration.

`atol` and `rtol` specify the absolute and relative tolerances for determining
convergence.

`order` is passed on to `quadgk()` of the
[QuadGK](https://github.com/JuliaMath/QuadGK.jl) package.

`nconvergences` is the number of iterations before convergence is declared.

See `?QuadOsc.accel_cohen_villegas_zagier` for details on the series
acceleration.
"""
function quadosc(f::Function, a::Number, b::Number, zerosf::Function; pren=2,
		 atol=zero(Float64), rtol=sqrt(eps(Float64)), order=7,
                 #nconvergences=ceil(Int,-1.31*log10(rtol)),
                 #accelerator=accel_cohen_villegas_zagier)
                 nconvergences=5,
                 accelerator=accel_wynn_eps)
    @assert b == Inf
    T = Float64

    i1 = findfirst(n -> zerosf(n) - a >= 0, 1:typemax(Int))
    z1 = zerosf(i1)
    Ipre, Epre = quadgk(f, a, z1; atol=atol, rtol=rtol, order=order)

    z0 = z1
    for i=1:pren
        i1 += 1
        z1 = zerosf(i1)
        I, E = quadgk(f, z0, z1; atol=atol, rtol=rtol, order=order)
        Ipre += I
        Epre += E
        z0 = z1
    end

    I = T(0)
    oldI = T(0)
    ak = T[]
    ek = T[]
    while nconvergences > 0
        i1 += 1
        z1 = zerosf(i1)
        I, E = quadgk(f, z0, z1; atol=atol, rtol=rtol)
        push!(ak, I)
        push!(ek, E)
        z0 = z1
        I = accelerator(ak)

        adiff = abs(I - oldI)
        rdiff = adiff * 2 / abs(I + oldI)
        if adiff <= atol || rdiff <= rtol
            nconvergences -= 1
        end
        oldI = I
    end
    I = Ipre + I
    E = Epre + abs(accelerator(ek))
    return I, E
end

end # module

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


text:
function calc_grad( pw::PWGrid, Potentials, Focc, W::Array{ComplexF64,2} )
    Ngwx    = size(W)[1]
    Nstates = size(W)[2]
    Ω = pw.Ω
    Ns = pw.Ns
    #
    grad = zeros( ComplexF64, Ngwx, Nstates )

    H_psi = op_H( pw, Potentials, W )
    for i = 1:Nstates
        grad[:,i] = H_psi[:,i]
        for j = 1:Nstates
            grad[:,i] = grad[:,i] - dot( W[:,j], H_psi[:,i] ) * W[:,j]
        end
        grad[:,i] = Focc[i]*grad[:,i]
    end

    #println("Pass here ...")

    F = Matrix(Diagonal(Focc))
    HW = op_H( pw, Potentials, W )
    ℍ = W' * HW
    HFH = ℍ*F - F*ℍ
    #denom = ones(Nstates,Nstates)*2.0
    #ℚ = HFH ./ denom
    ℚ = 0.5*HFH
    grad[:,:] = grad[:,:] + W*ℚ

    return grad
end

#function calc_grad( pw::PWGrid, Potentials, Focc, W::Array{ComplexF64,2} )
#  Ngwx    = size(W)[1]
#  Nstates = size(W)[2]
#  Ω = pw.Ω
#  Ns = pw.Ns
#  #
#  grad = zeros( ComplexF64, Ngwx, Nstates )
#
#  F = diagm(Focc)
#  HW = op_H( pw, Potentials, W )
#  ℍ = W' * HW
#  HFH = ℍ*F - F*ℍ
#  denom = ones(Nstates,Nstates)*2.0
#  ℚ = HFH ./ denom
#  grad = (HW - W * ℍ)*F + W*ℚ
#
#  return grad
#end
#