text:
begin
    using CSV
    using Revise
    using BNPRegressionGGA2021
    using DataFrames
    using DataStructures
    using Distributions
    using Random
    using RCall
    using Statistics
    using StatsBase
    using LinearAlgebra
    const BNP = BNPRegressionGGA2021
end;

function preprocess(dy, x0, x1)
    # Simulate the responses
    y0 = @. rand(dy(x0))
    y1 = LinRange(-3, 3, 50) |> collect

    # Expand the grid
    grid = rcopy(R"""expand.grid(y1 = $y1, x1 = $x1)""")
    y1 = grid[!, :y1]
    x1 = grid[!, :x1]

    # Expand x0-x1 using splines
    x0min, x0max = extrema(x0)
    R"""
    X0 <- splines::bs($x0, df = 6, Boundary.knots = c($x0min, $x0max))
    """
    @rget X0

    # Standardise responses and predictors
    mX0, sX0 = mean_and_std(X0, 1)
    for col in 1:size(X0, 2)
        if sX0[col] > 0
            X0[:, col] = (X0[:, col] .- mX0[col]) ./ sX0[col]
        end
    end
    my0, sy0 = mean_and_std(y0, 1)
    y0 = (y0 .- my0) ./ sy0

    # Add a constant term to the design matrices
    X0d = rand([0, 1], size(X0, 1), 4)
    X0d = ((X0d .- mean(X0d, dims = 1)) ./ std(X0d, dims = 1))
    X0 = [ones(size(X0, 1)) X0 X0d]
    X1 = zeros(0, size(X0, 2))

    # Generate the true (standardized) density over the grid
    f1 = pdf.(dy.(x1), my0 .+ y1 .* sy0) .* sy0

    # Set the mapping
    mapping = [[1], collect(2:7), [8], [9], [10], [11]]

    # Return the preprocessed results
    return y0, y1, x1, X0, X1, f1, mapping
end;

# Example 1: Continuous predictor (linear regression, normal distribution)
begin
    Random.seed!(1)
    N0, N1, Nrep, Niter = 200, 0, 20, 100
    dy(x) = Normal(x, 1)
    x1raw = LinRange(-1, 1, N1)
    x0raw = repeat(LinRange(-2, 2, N0 ÷ Nrep), Nrep)
    df = zeros(Bool, Niter, 6)
    # for iter in 1:Niter
    #     y0, y1, x1, X0, X1, f1, mapping = preprocess(dy, x0raw, x1raw)
    #     smpl = BNP.DGSBPNormalDependent(; y0, X0, y1, X1, mapping)
    #     _, _, chaing = BNP.sample!(smpl; mcmcsize = 10000)
    #     df[iter, :] .= begin
    #         chaing |>
    #         DataStructures.counter |>
    #         collect |>
    #         x -> sort(x, by = x -> x[2], rev = true) |>
    #         x -> getindex(x, 1) |>
    #         x -> getindex(x, 1)
    #     end
    # end
    # filename = "data/simulated-best-normal-bnp-ex-01-2022-03-31.csv"
    # CSV.write(filename, DataFrame(df, :auto))

    # Save the 100 simulated samples
    df = map(1:Niter) do iter
        y0, y1, x1, X0, X1, f1, mapping = preprocess(dy, x0raw, x1raw)
        df = DataFrame(X0, :auto)
        df[!, :iter] .= iter
        df[!, :y0] = y0
        df
    end;
    df = reduce(vcat, df)
    CSV.write("data/encargo-2022-03-31-data-ex1.csv", df);
end;

# Example 2: Continuous predictor (linear regression, skewnormal distribution)
begin
    Random.seed!(1)
    N0, N1, Nrep, Niter = 200, 0, 20, 100
    dy(x) = SkewNormal(2x, 4, 4)
    x1raw = LinRange(-1, 1, N1)
    x0raw = repeat(LinRange(-2, 2, N0 ÷ Nrep), Nrep)
    # df = zeros(Bool, Niter, 6)
    # for iter in 1:Niter
    #     y0, y1, x1, X0, X1, f1, mapping = preprocess(dy, x0raw, x1raw)
    #     smpl = BNP.DGSBPNormalDependent(; y0, X0, y1, X1, mapping)
    #     _, _, chaing = BNP.sample!(smpl; mcmcsize = 10000)
    #     df[iter, :] .= begin
    #         chaing |>
    #         DataStructures.counter |>
    #         collect |>
    #         x -> sort(x, by = x -> x[2], rev = true) |>
    #         x -> getindex(x, 1) |>
    #         x -> getindex(x, 1)
    #     end
    # end
    # filename = "data/simulated-best-normal-bnp-ex-02-2022-03-31.csv"
    # CSV.write(filename, DataFrame(df, :auto))

    # Save the 100 simulated samples
    df = map(1:Niter) do iter
        y0, y1, x1, X0, X1, f1, mapping = preprocess(dy, x0raw, x1raw)
        df = DataFrame(X0, :auto)
        df[!, :iter] .= iter
        df[!, :y0] = y0
        df
    end;
    df = reduce(vcat, df)
    CSV.write("data/encargo-2022-03-31-data-ex2.csv", df);
end;

# Example 3: Continuous predictor (cubic regression, normal distribution)
begin
    Random.seed!(1)
    N0, N1, Nrep, Niter = 200, 0, 20, 100
    dy(x) = Normal(0.2x^3, 0.5)
    x1raw = LinRange(-1, 1, N1)
    x0raw = repeat(LinRange(-2, 2, N0 ÷ Nrep), Nrep)
    df = zeros(Bool, Niter, 6)
    # for iter in 1:Niter
    #     y0, y1, x1, X0, X1, f1, mapping = preprocess(dy, x0raw, x1raw)
    #     smpl = BNP.DGSBPNormalDependent(; y0, X0, y1, X1, mapping)
    #     _, _, chaing = BNP.sample!(smpl; mcmcsize = 10000)
    #     df[iter, :] .= begin
    #         chaing |>
    #         DataStructures.counter |>
    #         collect |>
    #         x -> sort(x, by = x -> x[2], rev = true) |>
    #         x -> getindex(x, 1) |>
    #         x -> getindex(x, 1)
    #     end
    # end
    # filename = "data/simulated-best-normal-bnp-ex-03-2022-03-31.csv"
    # CSV.write(filename, DataFrame(df, :auto))

    # Save the 100 simulated samples
    df = map(1:Niter) do iter
        y0, y1, x1, X0, X1, f1, mapping = preprocess(dy, x0raw, x1raw)
        df = DataFrame(X0, :auto)
        df[!, :iter] .= iter
        df[!, :y0] = y0
        df
    end;
    df = reduce(vcat, df)
    CSV.write("data/encargo-2022-03-31-data-ex3.csv", df);
end;

# Example 4: Continuous predictor (cubic regression, mixture distribution)
begin
    Random.seed!(1)
    N0, N1, Nrep, Niter = 200, 0, 20, 100
    dy(x) = SkewNormal(0.4x^3, 4, 4)
    x1raw = LinRange(-1, 1, N1)
    x0raw = repeat(LinRange(-2, 2, N0 ÷ Nrep), Nrep)
    df = zeros(Bool, Niter, 6)
    # for iter in 1:Niter
    #     y0, y1, x1, X0, X1, f1, mapping = preprocess(dy, x0raw, x1raw)
    #     smpl = BNP.DGSBPNormalDependent(; y0, X0, y1, X1, mapping)
    #     _, _, chaing = BNP.sample!(smpl; mcmcsize = 10000)
    #     df[iter, :] .= begin
    #         chaing |>
    #         DataStructures.counter |>
    #         collect |>
    #         x -> sort(x, by = x -> x[2], rev = true) |>
    #         x -> getindex(x, 1) |>
    #         x -> getindex(x, 1)
    #     end
    # end
    # filename = "data/simulated-best-normal-bnp-ex-04-2022-03-31.csv"
    # CSV.write(filename, DataFrame(df, :auto))

    # Save the 100 simulated samples
    df = map(1:Niter) do iter
        y0, y1, x1, X0, X1, f1, mapping = preprocess(dy, x0raw, x1raw)
        df = DataFrame(X0, :auto)
        df[!, :iter] .= iter
        df[!, :y0] = y0
        df
    end;
    df = reduce(vcat, df)
    CSV.write("data/encargo-2022-03-31-data-ex4.csv", df);
end;


text:
struct SearchProblem <: AbstractProblem
    patternset::PatternSet
    trace::Vector{Float64}
    simulator::AbstractSimulator
    solver::AbstractSolver
end

ext: jl


text:
module TestCopy

using Test

import MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities

const DoubleDicts = MathOptInterface.Utilities.DoubleDicts

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

###
### Test helpers
###

abstract type AbstractDummyModel <: MOI.ModelLike end

function MOI.empty!(::AbstractDummyModel) end

function MOI.copy_to(dest::AbstractDummyModel, src::MOI.ModelLike)
    return MOIU.default_copy_to(dest, src)
end

MOI.supports(::AbstractDummyModel, ::MOI.ObjectiveSense) = true

function MOI.supports(
    ::AbstractDummyModel,
    ::MOI.ConstraintPrimalStart,
    ::Type{<:MOI.ConstraintIndex},
)
    return true
end

function MOI.supports_constraint(
    ::AbstractDummyModel,
    ::Type{MOI.VariableIndex},
    ::Type{MOI.EqualTo{Float64}},
)
    return true
end

function MOI.supports_constraint(
    ::AbstractDummyModel,
    ::Type{MOI.VectorOfVariables},
    ::Type{MOI.Zeros},
)
    return true
end

struct DummyModel <: AbstractDummyModel end

# Implements add_variable and add_constraint
struct DummyModelWithAdd <: AbstractDummyModel end

MOI.add_variable(::DummyModelWithAdd) = MOI.VariableIndex(0)

MOI.add_variables(::DummyModelWithAdd, n) = fill(MOI.VariableIndex(0), n)

function MOI.add_constraint(
    ::DummyModelWithAdd,
    ::MOI.VariableIndex,
    ::MOI.EqualTo{Float64},
)
    return MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{Float64}}(0)
end

struct DummyEvaluator <: MOI.AbstractNLPEvaluator end

###
### The tests
###

function _test_identity_index_map(::Type{T}) where {T}
    model = MOI.Utilities.Model{T}()
    x, y = MOI.add_variables(model, 2)
    cx = MOI.add_constraint(model, x, MOI.EqualTo(one(T)))
    cy = MOI.add_constraint(model, y, MOI.EqualTo(zero(T)))
    c = MOI.add_constraint(model, one(T) * x + y, MOI.LessThan(zero(T)))
    index_map = MOI.Utilities.identity_index_map(model)
    @test x == index_map[x]
    @test y == index_map[y]
    @test cx == index_map[cx]
    @test cy == index_map[cy]
    @test c == index_map[c]
    return
end

function test_identity_index_map()
    _test_identity_index_map(Int)
    _test_identity_index_map(Float64)
    return
end

function test_IndexMap()
    map = MOIU.IndexMap()
    @test length(map) == 0
    @test occursin("Utilities.IndexMap()", sprint(show, map))
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    cx = MOI.ConstraintIndex{MOI.VariableIndex,MOI.Integer}(1)
    cy = MOI.ConstraintIndex{MOI.VariableIndex,MOI.Integer}(2)
    map = MOIU.IndexMap(Dict(x => y), DoubleDicts.IndexDoubleDict())
    map[cx] = cy
    @test length(map) == 2
    return
end

function test_AUTOMATIC()
    src = DummyModel()
    dest = DummyModel()
    @test_throws ErrorException MOIU.default_copy_to(dest, src)
    try
        @test_throws ErrorException MOIU.default_copy_to(dest, src)
    catch err
        @test sprint(showerror, err) ==
              "Model DummyModel does not support copy with names."
    end
end

"""
    test_issue_849()

Create variables in same ordering when NLPBlock is used (#849)
"""
function test_issue_849()
    model = MOIU.UniversalFallback(MOIU.Model{Float64}())
    a = MOI.add_variable(model)
    b, c = MOI.add_variables(model, 2)
    x, cx = MOI.add_constrained_variable(model, MOI.GreaterThan(0.0))
    y, cy = MOI.add_constrained_variables(model, MOI.Nonnegatives(1))
    nlp_data = MOI.NLPBlockData(
        [MOI.NLPBoundsPair(0.0, 1.0) for i in 1:5],
        DummyEvaluator(),
        false,
    )
    MOI.set(model, MOI.NLPBlock(), nlp_data)
    copy = MOIU.UniversalFallback(MOIU.Model{Float64}())
    index_map = MOIU.default_copy_to(copy, model)
    for vi in [a, b, c, x, y[1]]
        @test index_map[vi] == vi
    end
    for ci in [cx, cy]
        @test index_map[ci] == ci
    end
    return
end

###
### ConstrainedVariablesModel
###

"""
    ConstrainedVariablesModel

Model that distinguish variables created constrained
"""
struct ConstrainedVariablesModel <: MOI.ModelLike
    added_constrained::Vector{Bool}
end

MOI.empty!(model::ConstrainedVariablesModel) = empty!(model.added_constrained)

MOI.supports_incremental_interface(::ConstrainedVariablesModel) = true

function MOI.copy_to(
    dest::ConstrainedVariablesModel,
    src::MOI.ModelLike;
    kwargs...,
)
    return MOIU.default_copy_to(dest, src; kwargs...)
end

function MOI.add_variables(model::ConstrainedVariablesModel, n)
    m = length(model.added_constrained)
    for _ in 1:n
        push!(model.added_constrained, false)
    end
    return MOI.VariableIndex.(m .+ (1:n))
end

function MOI.add_constrained_variables(
    model::ConstrainedVariablesModel,
    set::MOI.AbstractVectorSet,
)
    m = length(model.added_constrained)
    for _ in 1:MOI.dimension(set)
        push!(model.added_constrained, true)
    end
    ci = MOI.ConstraintIndex{MOI.VectorOfVariables,typeof(set)}(m + 1)
    return MOI.VariableIndex.(m .+ (1:MOI.dimension(set))), ci
end

function MOI.add_constraint(
    ::ConstrainedVariablesModel,
    func::MOI.VectorOfVariables,
    set::MOI.AbstractVectorSet,
)
    return MOI.ConstraintIndex{typeof(func),typeof(set)}(
        func.variables[1].value,
    )
end

function test_ConstrainedVariablesModel()
    src = MOIU.Model{Int}()
    x = MOI.add_variables(src, 3)
    cx = MOI.add_constraint(src, [x[1], x[3], x[1], x[2]], MOI.Nonnegatives(4))
    y, cy = MOI.add_constrained_variables(src, MOI.Nonpositives(3))
    dest = ConstrainedVariablesModel(Bool[])
    idxmap = MOI.copy_to(dest, src)
    for vi in x
        @test !dest.added_constrained[idxmap[vi].value]
    end
    for vi in y
        @test dest.added_constrained[idxmap[vi].value]
    end
    return
end

###
### AbstractConstrainedVariablesModel
###

abstract type AbstractConstrainedVariablesModel <: MOI.ModelLike end

mutable struct OrderConstrainedVariablesModel <:
               AbstractConstrainedVariablesModel
    constraintIndices::Array{MOI.ConstraintIndex}
    inner::MOIU.Model{Float64}
    function OrderConstrainedVariablesModel()
        return new(MOI.ConstraintIndex[], MOIU.Model{Float64}())
    end
end
mutable struct ReverseOrderConstrainedVariablesModel <:
               AbstractConstrainedVariablesModel
    constraintIndices::Array{MOI.ConstraintIndex}
    inner::MOIU.Model{Float64}
    function ReverseOrderConstrainedVariablesModel()
        return new(MOI.ConstraintIndex[], MOIU.Model{Float64}())
    end
end

function MOI.add_variables(model::AbstractConstrainedVariablesModel, n)
    return MOI.add_variables(model.inner, n)
end

function MOI.add_variable(model::AbstractConstrainedVariablesModel)
    return MOI.add_variable(model.inner)
end

function MOI.add_constraint(
    model::AbstractConstrainedVariablesModel,
    f::MOI.AbstractFunction,
    s::MOI.AbstractSet,
)
    ci = MOI.add_constraint(model.inner, f, s)
    push!(model.constraintIndices, ci)
    return ci
end

function MOI.copy_to(
    dest::AbstractConstrainedVariablesModel,
    src::MOI.ModelLike;
    kwargs...,
)
    return MOIU.default_copy_to(dest, src; kwargs...)
end

MOI.supports_incremental_interface(::AbstractConstrainedVariablesModel) = true

function MOI.empty!(model::AbstractConstrainedVariablesModel)
    model.constraintIndices = MOI.ConstraintIndex[]
    MOI.empty!(model.inner)
    return
end

function MOI.supports_constraint(
    ::OrderConstrainedVariablesModel,
    ::Type{MOI.VectorOfVariables},
    ::Type{MOI.Nonnegatives},
)
    return false
end

function MOI.supports_add_constrained_variables(
    ::OrderConstrainedVariablesModel,
    ::Type{MOI.Nonnegatives},
)
    return true
end

function MOI.supports_constraint(
    ::OrderConstrainedVariablesModel,
    ::Type{MOI.VectorOfVariables},
    ::Type{MOI.Nonnegatives},
)
    return true
end

function MOI.supports_add_constrained_variables(
    ::OrderConstrainedVariablesModel,
    ::Type{MOI.Nonpositives},
)
    return false
end

function MOI.supports_constraint(
    ::ReverseOrderConstrainedVariablesModel,
    ::Type{MOI.VectorOfVariables},
    ::Type{MOI.Nonnegatives},
)
    return true
end

function MOI.supports_add_constrained_variables(
    ::ReverseOrderConstrainedVariablesModel,
    ::Type{MOI.Nonnegatives},
)
    return false
end

function MOI.supports_constraint(
    ::ReverseOrderConstrainedVariablesModel,
    ::Type{MOI.VectorOfVariables},
    ::Type{MOI.Nonnegatives},
)
    return false
end

function MOI.supports_add_constrained_variables(
    ::ReverseOrderConstrainedVariablesModel,
    ::Type{MOI.Nonpositives},
)
    return true
end

function MOI.supports_constraint(
    ::OrderConstrainedVariablesModel,
    ::Type{MOI.VectorAffineFunction{Float64}},
    ::Type{MOI.Nonnegatives},
)
    return true
end

function MOI.supports_constraint(
    ::ReverseOrderConstrainedVariablesModel,
    ::Type{MOI.VectorAffineFunction{Float64}},
    ::Type{MOI.Nonpositives},
)
    return true
end

function MOI.supports_constraint(
    ::OrderConstrainedVariablesModel,
    ::Type{MOI.VariableIndex},
    ::Type{<:MOI.GreaterThan},
)
    return true
end

function MOI.supports_add_constrained_variable(
    ::OrderConstrainedVariablesModel,
    ::Type{<:MOI.GreaterThan},
)
    return false
end

function MOI.supports_constraint(
    ::OrderConstrainedVariablesModel,
    ::Type{MOI.VariableIndex},
    ::Type{<:MOI.LessThan},
)
    return false
end

function MOI.supports_add_constrained_variable(
    ::OrderConstrainedVariablesModel,
    ::Type{<:MOI.LessThan},
)
    return true
end

function MOI.supports_constraint(
    ::ReverseOrderConstrainedVariablesModel,
    ::Type{MOI.VariableIndex},
    ::Type{<:MOI.GreaterThan},
)
    return false
end

function MOI.supports_add_constrained_variable(
    ::ReverseOrderConstrainedVariablesModel,
    ::Type{<:MOI.GreaterThan},
)
    return true
end

function MOI.supports_constraint(
    ::ReverseOrderConstrainedVariablesModel,
    ::Type{MOI.VariableIndex},
    ::Type{<:MOI.LessThan},
)
    return true
end

function MOI.supports_add_constrained_variable(
    ::ReverseOrderConstrainedVariablesModel,
    ::Type{<:MOI.LessThan},
)
    return false
end

"""
    test_create_variables_using_supports_add_constrained_variable()

From Issue #987
"""
function test_create_variables_using_supports_add_constrained_variable()
    # With vectors
    src = MOIU.Model{Float64}()
    a, c1 = MOI.add_constrained_variables(src, MOI.Nonpositives(3))
    c2 = MOI.add_constraint(src, a, MOI.Nonnegatives(3))

    dest = OrderConstrainedVariablesModel()
    index_map = MOI.copy_to(dest, src)
    @test typeof(c1) == typeof(dest.constraintIndices[2])
    @test typeof(c2) == typeof(dest.constraintIndices[1])

    dest = ReverseOrderConstrainedVariablesModel()
    index_map = MOI.copy_to(dest, src)
    @test typeof(c1) == typeof(dest.constraintIndices[1])
    @test typeof(c2) == typeof(dest.constraintIndices[2])

    b, cb = MOI.add_constrained_variables(src, MOI.Nonnegatives(2))
    c3 = MOI.add_constraint(src, b, MOI.Zeros(2))

    d, cd = MOI.add_constrained_variables(src, MOI.Zeros(2))
    c4 = MOI.add_constraint(src, d, MOI.Nonpositives(2))

    dest = OrderConstrainedVariablesModel()
    bridged_dest = MOI.Bridges.full_bridge_optimizer(dest, Float64)
    @test MOIU.sorted_variable_sets_by_cost(bridged_dest, src) ==
          Type[MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives]
    @test MOI.supports_add_constrained_variables(bridged_dest, MOI.Nonnegatives)
    @test MOI.get(bridged_dest, MOI.VariableBridgingCost{MOI.Nonnegatives}()) ==
          0.0
    @test MOI.supports_constraint(
        bridged_dest,
        MOI.VectorOfVariables,
        MOI.Nonnegatives,
    )
    @test MOI.get(
        bridged_dest,
        MOI.ConstraintBridgingCost{MOI.VectorOfVariables,MOI.Nonnegatives}(),
    ) == 0.0
    @test MOI.supports_add_constrained_variables(bridged_dest, MOI.Nonpositives)
    @test MOI.get(bridged_dest, MOI.VariableBridgingCost{MOI.Nonpositives}()) ==
          1.0
    @test MOI.supports_constraint(
        bridged_dest,
        MOI.VectorOfVariables,
        MOI.Nonpositives,
    )
    @test MOI.get(
        bridged_dest,
        MOI.ConstraintBridgingCost{MOI.VectorOfVariables,MOI.Nonpositives}(),
    ) == 1.0
    @test MOI.supports_add_constrained_variables(bridged_dest, MOI.Zeros)
    @test MOI.get(bridged_dest, MOI.VariableBridgingCost{MOI.Zeros}()) == 1.0
    @test MOI.supports_constraint(
        bridged_dest,
        MOI.VectorOfVariables,
        MOI.Zeros,
    )
    @test MOI.get(
        bridged_dest,
        MOI.ConstraintBridgingCost{MOI.VectorOfVariables,MOI.Zeros}(),
    ) == 2.0
    index_map = MOI.copy_to(bridged_dest, src)
    @test length(dest.constraintIndices) == 4

    dest = ReverseOrderConstrainedVariablesModel()
    bridged_dest = MOI.Bridges.full_bridge_optimizer(dest, Float64)
    @test MOIU.sorted_variable_sets_by_cost(bridged_dest, src) ==
          Type[MOI.Zeros, MOI.Nonpositives, MOI.Nonnegatives]
    @test MOI.supports_add_constrained_variables(bridged_dest, MOI.Nonnegatives)
    @test MOI.get(bridged_dest, MOI.VariableBridgingCost{MOI.Nonnegatives}()) ==
          2.0
    @test MOI.supports_constraint(
        bridged_dest,
        MOI.VectorOfVariables,
        MOI.Nonnegatives,
    )
    @test MOI.get(
        bridged_dest,
        MOI.ConstraintBridgingCost{MOI.VectorOfVariables,MOI.Nonnegatives}(),
    ) == 1.0
    @test MOI.supports_add_constrained_variables(bridged_dest, MOI.Nonpositives)
    @test MOI.get(bridged_dest, MOI.VariableBridgingCost{MOI.Nonpositives}()) ==
          0.0
    @test MOI.supports_constraint(
        bridged_dest,
        MOI.VectorOfVariables,
        MOI.Nonpositives,
    )
    @test MOI.get(
        bridged_dest,
        MOI.ConstraintBridgingCost{MOI.VectorOfVariables,MOI.Nonpositives}(),
    ) == 1.0
    @test MOI.supports_add_constrained_variables(bridged_dest, MOI.Zeros)
    @test MOI.get(bridged_dest, MOI.VariableBridgingCost{MOI.Zeros}()) == 1.0
    @test MOI.supports_constraint(
        bridged_dest,
        MOI.VectorOfVariables,
        MOI.Zeros,
    )
    @test MOI.get(
        bridged_dest,
        MOI.ConstraintBridgingCost{MOI.VectorOfVariables,MOI.Zeros}(),
    ) == 3.0
    index_map = MOI.copy_to(bridged_dest, src)
    @test length(dest.constraintIndices) == 4

    # With single variables
    src = MOIU.Model{Float64}()
    a, c1 = MOI.add_constrained_variable(src, MOI.GreaterThan{Float64}(5.0))
    c2 = MOI.add_constraint(src, a, MOI.LessThan{Float64}(1.0))

    dest = OrderConstrainedVariablesModel()
    index_map = MOI.copy_to(dest, src)
    @test typeof(c1) == typeof(dest.constraintIndices[2])
    @test typeof(c2) == typeof(dest.constraintIndices[1])

    dest = ReverseOrderConstrainedVariablesModel()
    index_map = MOI.copy_to(dest, src)
    @test typeof(c1) == typeof(dest.constraintIndices[1])
    @test typeof(c2) == typeof(dest.constraintIndices[2])

    dest = OrderConstrainedVariablesModel()
    bridged_dest = MOI.Bridges.full_bridge_optimizer(dest, Float64)
    @test MOI.get(
        bridged_dest,
        MOI.VariableBridgingCost{MOI.LessThan{Float64}}(),
    ) == 0.0
    @test MOI.get(
        bridged_dest,
        MOI.ConstraintBridgingCost{MOI.VariableIndex,MOI.LessThan{Float64}}(),
    ) == 2.0
    @test MOI.get(
        bridged_dest,
        MOI.VariableBridgingCost{MOI.GreaterThan{Float64}}(),
    ) == 1.0
    @test MOI.get(
        bridged_dest,
        MOI.ConstraintBridgingCost{MOI.VariableIndex,MOI.GreaterThan{Float64}}(),
    ) == 0.0
    index_map = MOI.copy_to(bridged_dest, src)
    @test typeof(c1) == typeof(dest.constraintIndices[2])
    @test typeof(c2) == typeof(dest.constraintIndices[1])

    dest = ReverseOrderConstrainedVariablesModel()
    bridged_dest = MOI.Bridges.full_bridge_optimizer(dest, Float64)
    @test MOI.get(
        bridged_dest,
        MOI.VariableBridgingCost{MOI.LessThan{Float64}}(),
    ) == 1.0
    @test MOI.get(
        bridged_dest,
        MOI.ConstraintBridgingCost{MOI.VariableIndex,MOI.LessThan{Float64}}(),
    ) == 0.0
    @test MOI.get(
        bridged_dest,
        MOI.VariableBridgingCost{MOI.GreaterThan{Float64}}(),
    ) == 0.0
    @test MOI.get(
        bridged_dest,
        MOI.ConstraintBridgingCost{MOI.VariableIndex,MOI.GreaterThan{Float64}}(),
    ) == 2.0
    index_map = MOI.copy_to(bridged_dest, src)
    @test typeof(c1) == typeof(dest.constraintIndices[1])
    @test typeof(c2) == typeof(dest.constraintIndices[2])
end

function test_filtering_copy()
    # Create a basic model.
    src = MOIU.Model{Float64}()
    x = MOI.add_variable(src)
    c1 = MOI.add_constraint(src, x, MOI.LessThan{Float64}(1.0))
    c2 = MOI.add_constraint(src, x, MOI.GreaterThan{Float64}(0.0))

    # Filtering function: the default case where this function always returns
    # true is already well-tested by the above cases.
    # Only keep the constraint c1.
    f(::Any) = true
    f(cidx::MOI.ConstraintIndex) = cidx == c1

    # Perform the copy.
    dst = OrderConstrainedVariablesModel()
    index_map = MOI.copy_to(dst, MOI.Utilities.ModelFilter(f, src))

    @test typeof(c1) == typeof(dst.constraintIndices[1])
    @test length(dst.constraintIndices) == 1
    return
end

function test_filtering_copy_empty()
    src = MOI.Utilities.Model{Float64}()
    MOI.set(src, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    model = MOI.Utilities.ModelFilter(src) do item
        return item !== MOI.ObjectiveSense()
    end
    @test MOI.is_empty(model)
    MOI.add_variable(src)
    @test !MOI.is_empty(model)
    MOI.empty!(model)
    @test MOI.is_empty(model)
    MOI.set(src, MOI.Name(), "Model")
    @test !MOI.is_empty(model)
    return
end

function test_filtering_copy_AbstractModelAttribute()
    src = MOI.Utilities.Model{Float64}()
    MOI.set(src, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    dest = MOI.Utilities.Model{Float64}()
    filtered_src = MOI.Utilities.ModelFilter(src) do item
        if item isa MOI.AbstractModelAttribute
            return item != MOI.ObjectiveSense()
        end
        return true
    end
    @test MOI.is_empty(filtered_src)
    MOI.copy_to(dest, filtered_src)
    @test MOI.get(dest, MOI.ObjectiveSense()) == MOI.FEASIBILITY_SENSE
    return
end

function test_filtering_copy_AbstractVariableAttribute()
    src = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}())
    x = MOI.add_variable(src)
    MOI.set(src, MOI.VariableName(), x, "x")
    MOI.set(src, MOI.VariablePrimalStart(), x, 1.0)
    dest = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}())
    filtered_src = MOI.Utilities.ModelFilter(src) do item
        if item isa MOI.AbstractVariableAttribute
            return item != MOI.VariableName()
        end
        return true
    end
    @test !MOI.is_empty(filtered_src)
    index_map = MOI.copy_to(dest, filtered_src)
    @test MOI.get(dest, MOI.VariableName(), index_map[x]) == ""
    @test MOI.get(dest, MOI.VariablePrimalStart(), index_map[x]) == 1.0
    return
end

function test_filtering_copy_AbstractConstraintAttribute()
    src = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}())
    x = MOI.add_variable(src)
    c = MOI.add_constraint(src, MOI.VectorOfVariables([x]), MOI.Nonnegatives(1))
    MOI.set(src, MOI.ConstraintName(), c, "c")
    MOI.set(src, MOI.ConstraintDualStart(), c, [1.0])
    dest = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}())
    filtered_src = MOI.Utilities.ModelFilter(src) do item
        if item isa MOI.AbstractConstraintAttribute
            return item != MOI.ConstraintName()
        end
        return true
    end
    @test !MOI.is_empty(filtered_src)
    index_map = MOI.copy_to(dest, filtered_src)
    @test MOI.get(dest, MOI.ConstraintName(), index_map[c]) == ""
    @test MOI.get(dest, MOI.ConstraintDualStart(), index_map[c]) == [1.0]
    return
end

function test_filtering_copy_ListOfConstraintIndices()
    src = MOI.Utilities.Model{Float64}()
    x = MOI.add_variables(src, 4)
    MOI.add_constraint.(src, x, MOI.GreaterThan(0.0))
    dest = MOI.Utilities.Model{Float64}()
    filtered_src = MOI.Utilities.ModelFilter(src) do item
        if item isa MOI.ConstraintIndex
            return isodd(item.value)
        end
        return true
    end
    @test !MOI.is_empty(filtered_src)
    index_map = MOI.copy_to(dest, filtered_src)
    @test MOI.get(dest, MOI.NumberOfVariables()) == 4
    @test MOI.get(
        dest,
        MOI.NumberOfConstraints{MOI.VariableIndex,MOI.GreaterThan{Float64}}(),
    ) == 2
    @test MOI.get(
        filtered_src,
        MOI.NumberOfConstraints{MOI.VariableIndex,MOI.GreaterThan{Float64}}(),
    ) == 2
    for xi in x
        @test haskey(index_map, xi)
        ci = MOI.ConstraintIndex{MOI.VariableIndex,MOI.GreaterThan{Float64}}(
            xi.value,
        )
        @test haskey(index_map, ci) == isodd(xi.value)
    end
    return
end

function test_filtering_copy_ListOfConstraintTypesPresent()
    src = MOI.Utilities.Model{Float64}()
    x = MOI.add_variables(src, 4)
    MOI.add_constraint.(src, x, MOI.GreaterThan(0.0))
    MOI.add_constraint.(src, x, MOI.ZeroOne())
    dest = MOI.Utilities.Model{Float64}()
    index_map = MOI.copy_to(
        dest,
        MOI.Utilities.ModelFilter(src) do item
            if item isa Tuple{Type,Type}
                return item == (MOI.VariableIndex, MOI.ZeroOne)
            end
            return true
        end,
    )
    @test MOI.get(dest, MOI.NumberOfVariables()) == 4
    @test MOI.get(
        dest,
        MOI.NumberOfConstraints{MOI.VariableIndex,MOI.GreaterThan{Float64}}(),
    ) == 0
    @test MOI.get(
        dest,
        MOI.NumberOfConstraints{MOI.VariableIndex,MOI.ZeroOne}(),
    ) == 4
    return
end

###
### BoundModel
###
mutable struct BoundModel <: MOI.ModelLike
    # Type of model that only supports ≤ bound constraints. In particular, it
    # does not support integrality constraints.
    inner::MOIU.Model{Float64}
    BoundModel() = new(MOIU.Model{Float64}())
end

MOI.add_variable(model::BoundModel) = MOI.add_variable(model.inner)

function MOI.add_constraint(
    model::BoundModel,
    f::MOI.AbstractFunction,
    s::MOI.LessThan{Float64},
)
    return MOI.add_constraint(model.inner, f, s)
end

function MOI.supports_constraint(
    ::BoundModel,
    ::Type{MOI.VariableIndex},
    ::MOI.LessThan{Float64},
)
    return true
end

MOI.supports_incremental_interface(::BoundModel) = true

function MOI.copy_to(dest::BoundModel, src::MOI.ModelLike; kwargs...)
    return MOIU.default_copy_to(dest, src; kwargs...)
end

MOI.empty!(model::BoundModel) = MOI.empty!(model.inner)

MOI.supports(::BoundModel, ::Type{MOI.NumberOfConstraints}) = true

function MOI.get(model::BoundModel, attr::MOI.NumberOfConstraints)
    return MOI.get(model.inner, attr)
end

function test_BoundModel_filtering_copy()
    # Create a basic model.
    src = MOIU.Model{Float64}()
    x = MOI.add_variable(src)
    c1 = MOI.add_constraint(src, x, MOI.LessThan{Float64}(10.0))
    c2 = MOI.add_constraint(src, x, MOI.Integer())

    # Filtering function: get rid of integrality constraint.
    f(::Any) = true
    f(::MOI.ConstraintIndex{F,MOI.Integer}) where {F} = false

    # Perform the unfiltered copy. This should throw an error (i.e. the
    # implementation of BoundModel
    # should be correct).
    dst = BoundModel()
    @test_throws(
        MOI.UnsupportedConstraint{MOI.VariableIndex,MOI.Integer},
        MOI.copy_to(dst, src),
    )

    # Perform the filtered copy. This should not throw an error.
    dst = BoundModel()
    MOI.copy_to(dst, MOI.Utilities.ModelFilter(f, src))
    @test MOI.get(
        dst,
        MOI.NumberOfConstraints{MOI.VariableIndex,MOI.LessThan{Float64}}(),
    ) == 1
    @test MOI.get(
        dst,
        MOI.NumberOfConstraints{MOI.VariableIndex,MOI.Integer}(),
    ) == 0
end

# We create a `OnlyCopyConstraints` that don't implement `add_constraint` but
# implements `pass_nonvariable_constraints` to check that this is passed accross
# all layers without falling back to `pass_nonvariable_constraints_fallback`
# which calls `add_constraint`.

struct OnlyCopyConstraints{F,S} <: MOI.ModelLike
    constraints::MOIU.VectorOfConstraints{F,S}
    function OnlyCopyConstraints{F,S}() where {F,S}
        return new{F,S}(MOIU.VectorOfConstraints{F,S}())
    end
end

function MOI.Utilities._add_variable(::OnlyCopyConstraints) end

MOI.empty!(model::OnlyCopyConstraints) = MOI.empty!(model.constraints)

function MOI.supports_constraint(
    model::OnlyCopyConstraints,
    F::Type{<:MOI.AbstractFunction},
    S::Type{<:MOI.AbstractSet},
)
    return MOI.supports_constraint(model.constraints, F, S)
end

function MOIU.pass_nonvariable_constraints(
    dest::OnlyCopyConstraints,
    src::MOI.ModelLike,
    idxmap::MOIU.IndexMap,
    constraint_types,
)
    return MOIU.pass_nonvariable_constraints(
        dest.constraints,
        src,
        idxmap,
        constraint_types,
    )
end

function _test_copy_of_constraints_passed_as_copy_accross_layers(
    ::Type{T},
) where {T}
    F = MOI.ScalarAffineFunction{T}
    S = MOI.EqualTo{T}
    S2 = MOI.GreaterThan{T}
    src = MOIU.Model{T}()
    x = MOI.add_variable(src)
    MOI.add_constraint(src, T(1) * x, MOI.EqualTo(T(1)))
    MOI.add_constraint(src, T(2) * x, MOI.EqualTo(T(2)))
    MOI.add_constraint(src, T(3) * x, MOI.GreaterThan(T(3)))
    MOI.add_constraint(src, T(4) * x, MOI.GreaterThan(T(4)))
    dest = MOIU.CachingOptimizer(
        MOI.Bridges.full_bridge_optimizer(
            MOIU.UniversalFallback(
                MOIU.GenericOptimizer{
                    T,
                    MOI.Utilities.ObjectiveContainer{T},
                    MOI.Utilities.VariablesContainer{T},
                    OnlyCopyConstraints{F,S},
                }(),
            ),
            T,
        ),
        MOIU.AUTOMATIC,
    )
    MOI.copy_to(dest, src)
    voc = dest.model_cache.model.model.constraints.constraints
    @test MOI.get(voc, MOI.NumberOfConstraints{F,S}()) == 2
    @test !haskey(dest.model_cache.model.constraints, (F, S))
    @test MOI.get(dest, MOI.NumberOfConstraints{F,S2}()) == 2
    @test haskey(dest.model_cache.model.constraints, (F, S2))
    return
end

function test_copy_of_constraints_passed_as_copy_accross_layers()
    _test_copy_of_constraints_passed_as_copy_accross_layers(Int)
    _test_copy_of_constraints_passed_as_copy_accross_layers(Float64)
    return
end

end  # module

TestCopy.runtests()


text:
export initRemoteChannel

"""
function ref = initRemoteChannel(f,pid,args...)

Runs the function f with arguments args on worker

input:
       func::Union{Function,Type} -- function or type constructor
        pid::Int64 -- worker that will run f and store result
       args::varargs -- Arguments of f. These must be available on worker pid
       kwargs::varargs -- Key word arguments of f. These must be available on worker pid

output:
        ref::RemoteChannel -- Reusable reference to output of f.

"""
function initRemoteChannel(func::Union{Function,Type}, pid::Int64, args...; kwargs...)
  return RemoteChannel(()->initChannel(func,args,kwargs), pid)
end

function initChannel(func::Union{Function,Type},args::Tuple,kwargs)
  obj = func(args...; kwargs...)
  chan = Channel{typeof(obj)}(1)
  put!(chan,obj)
  return chan
end


text:
module Chordal

using SparseArrays, SuiteSparse, LinearAlgebra
using JuMP
using DataStructures

include("types.jl")
include("utils.jl")
include("chordal_decomp.jl")
include("chordal_extension.jl")
include("chordal_graph.jl")
include("clique_graph.jl")
include("clique_tree.jl")
include("completion.jl")
include("edm.jl")
include("jump.jl")

export sparsity_pattern, get_chordal_extension
export generate_clique_graph, merge_cliques!, get_cliques
export get_selectors, make_selectors_from_cliques, make_selectors_from_clique_graph
export maxdet_completion, minrank_completion
export build_constraints_lmi!
export get_sparsity_pattern_from_cliques, reconstruct_from_sparse_varref
export edm_completion

end


text:
@naiveVariationalRule(:node_type     => Softmax,
                      :outbound_type => Message{Categorical},
                      :inbound_types => (Nothing, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution),
                      :name          => VBSoftmaxOut)

@naiveVariationalRule(:node_type     => Softmax,
                      :outbound_type => Message{GaussianWeightedMeanPrecision},
                      :inbound_types => (ProbabilityDistribution, Nothing, ProbabilityDistribution, ProbabilityDistribution),
                      :name          => VBSoftmaxIn1)

@naiveVariationalRule(:node_type     => Softmax,
                      :outbound_type => Message{Function},
                      :inbound_types => (ProbabilityDistribution, ProbabilityDistribution, Nothing, ProbabilityDistribution),
                      :name          => VBSoftmaxXi)

@naiveVariationalRule(:node_type     => Softmax,
                      :outbound_type => Message{Function},
                      :inbound_types => (ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, Nothing),
                      :name          => VBSoftmaxA)


text:
#function liftinputalg(s::DiscreteLinearControlSystem)
#    n = HybridSystems.dim(s)
#    z = find(i -> iszero(sum(abs.(s.B[i,:]))), 1:n)
#    # TODO the affine space may not be parallel to classical axis
#    DiscreteLinearAlgebraicSystem(s.A[z, :], (eye(n))[z, :])
#end
#
#import MultivariatePolynomials
#const MP = MultivariatePolynomials
#
## Unbounded input
#function liftinputalg(s::ConstrainedSystem{<:AbstractSystem, LightAutomaton{P, MP.FullSpace}}, c) where P
#    ConstrainedSystem(liftinputalg(s.s), s.a), c
#end
#function liftinputalg(s::ConstrainedSystem, c)
#    n = HybridSystems.dim(s.s)
#    nu = size(s.s.B, 2)
#    A = [s.s.A s.s.B;
#         zeros(nu, n+nu)]
#    B = [zeros(n, nu); eye(nu)]
#    newc = [[c[v]; zeros(nu)] for v in vertices(s.a.G)]
#    newinv = [s.a.inv[v] * s.a.inputbounds[v] for v in vertices(s.a.G)]
#    newinb = [MP.FullSpace() for v in vertices(s.a.G)]
#    newa = Automaton(s.a.G, newinv, newinb)
#    liftinputalg(ConstrainedSystem(LinearControlSystem(A, B), newa), newc)
#end
#
#function liftmat(Q, o, c)
#    [o*o'-Q c-o; c'-o' 1]
#end
#function liftmat(Q, o, c, A)
#    @show Q
#    @show c
#    @show A
#    liftmat(A*Q*A', A*o, A*c)
#end
#
#using Convex
#using SCS
#using Polyhedra
#export getcis
#function getcis(sin::ConstrainedSystem, c, uc=false)
#    s, c = liftinputalg(sin, c)
#    as = s.s
#    @show as.A
#    @show as.E
#    g = s.a.G
#    n = LightGraphs.nv(g)
#    d = HybridSystems.dim(as)
#    Q = [Semidefinite(d) for v in vertices(g)]
#    if uc
#        C = [Variable(d) for v in vertices(g)]
#    else
#        C = [zeros(d) for v in vertices(g)]
#    end
#
##   Sol = [diagm([1, 1, 2])
##   constr1 = Q >= eye(dim(s))
#    constrs = Constraint[]
#    for u in vertices(g)
#        m = outdegree(g, u)
#        constr = liftmat(Q[u], c[u], C[u], as.A)
##       test = liftmat(eye(d), c[u], as.A)
##       display(liftmat(eye(d), c[u], as.A))
#        for v in out_neighbors(g, u)
#            constr -= 1/m * liftmat(Q[v], c[v], C[v], as.E)
##           display(liftmat(eye(d), c[v], as.E))
##           test -= 1/m * liftmat(eye(d), c[v], as.E)
#        end
##       @show constr
##       println("test")
##       display(test)
##       println()
##       @show eigmin(test)
#        push!(constrs, constr >= 0)
#        P = liftmat(Q[u], c[u], C[u])
#        for hs in ineqs(s.a.inv[u])
#            a = [hs.a; hs.β]
#            push!(constrs, a' * P * a >= 0)
#        end
#    end
#
#    problem = maximize(sum(logdet(Q[i]) for i in 1:n), constrs)
#
#    solve!(problem, SCSSolver())
#
#    @show problem.status
#    @show problem.optval
#    @show sum(i -> log(det(Q[i].value)), 1:n)
#
#    for i in vertices(g)
#        display(Q[i].value)
#        if uc
#            display(liftmat(Q[i].value, c[i], C[i].value))
#        end
#    end
#    if uc
#        [Ellipsoid(inv(Q[u].value), c[u] + 2*reshape(C[u].value, d)) for u in vertices(g)]
#    else
#        [Ellipsoid(inv(Q[u].value), c[u]) for u in vertices(g)]
#    end
#end


text:
module NeuralSimulation

using DiffEqBase
export CableSection, HHCableFunction, BEHHCable, HHCableProblem, HinesMatrix, init, solve, step!

abstract type AbstractNeuronalAlgorithm <: DiffEqBase.DEAlgorithm end

# Eventually we might just use `ImplicitEuler` as the algorithm type, and just bake a different integrator based
# on the unique method signature
struct BackwardsEuler <: AbstractNeuronalAlgorithm end
export BackwardsEuler

# Placeholder geometry type, this will be it's own major portion of the package/API eventually
abstract type AbstractGeometry end

struct CableSection{T,N} <: AbstractGeometry
    radius::T
    len::T
    rl::T
    cm::T
    num::N
end

mutable struct HinesMatrix{S, Alg}
    L::S
    D::S
    U::S
end

# TODO: other constructors for more elaborate tree structure types
# this version is specific to backwards Euler
function HinesMatrix(cable::CableSection, alg::BackwardsEuler, dt)
        rl = cable.rl
        a = cable.radius
        cm = cable.cm
        n = cable.num
        dx = cable.len/n

        L = Vector{typeof(cable.len)}(undef, n)
        D = similar(L)
        U = similar(L)

        # L and U are static
        L[1] = 0.0
        L[2:end-1] .= -(a*dt)/(2*rl*cm*dx^2)
        L[end] = -(a*dt)/(rl*cm*dx^2)

        U[1] = -(a*dt)/(rl*cm*dx^2)
        U[2:end-1] .= -(a*dt)/(2*rl*cm*dx^2)
        U[end] = 0.0

        return HinesMatrix{typeof(L), typeof(alg)}(L, D, U)
end

abstract type AbstractNeuronalFunction{iip} <: DiffEqBase.AbstractDiffEqFunction{iip} end

# A NeuronalSystem containing equations should be rewritten into more efficient forms for certain methods.
# e.g. rewrites for solving backwards Euler in fewer steps (as we have done below)
struct HHCableFunction{iip, F1, F2} <: AbstractNeuronalFunction{iip}
    f1::F1  # in-place update of coefficients for Hines matrix
    f2::F2  # in-place Channel gating dynamics
    function HHCableFunction{iip}(f1, f2) where {iip}
        new{iip, typeof(f1), typeof(f2)}(f1,f2)
    end
end

abstract type AbstractNeuronalIntegrator{Alg, IIP, U, T} <: DiffEqBase.DEIntegrator{Alg, IIP, U, T} end

mutable struct BackwardsEulerHHCableIntegrator{IIP, S, T, P, F, HinesMatrix} <: AbstractNeuronalIntegrator{BackwardsEuler, IIP, S, T}
    f::F     # equations
    uprev::S # Previous state
    u::S     # Current state
    tprev::T # Previous time step
    t::T     # Current time step
    t0::T    # Initial time step, only for re-initialization
    dt::T    # Step size
    p::P     # Parameters container
    H::HinesMatrix # Hines matrix
end

const BEHHCable = BackwardsEulerHHCableIntegrator
DiffEqBase.isinplace(::HHCableFunction{IIP}) where {IIP} = IIP
DiffEqBase.isinplace(::BEHHCable{IIP}) where {IIP} = IIP

abstract type AbstractNeuronalProblem{uType, tType, isinplace} <: DiffEqBase.DEProblem end

struct HHCableProblem{uType, tType, isinplace, P, F, G} <: AbstractNeuronalProblem{uType, tType, isinplace}
    f::F
    u0::uType
    tspan::tType
    p::P # constants to be supplied as second arg of `f`
    geometry::G
    function HHCableProblem(f::AbstractNeuronalFunction, cable::AbstractGeometry, u0, tspan, p)
        new{typeof(u0), typeof(tspan), DiffEqBase.isinplace(f), typeof(p), typeof(f), typeof(cable)}(f, u0, tspan, p, cable)
    end
end

DiffEqBase.isinplace(::HHCableProblem{U, T, IIP, P, F}) where {U, T, IIP, P, F} = IIP

# Initialization

function DiffEqBase.solve(prob::AbstractNeuronalProblem, args...; kwargs...)
    DiffEqBase.__solve(prob,args...; kwargs...)
end

function DiffEqBase.init(prob::AbstractNeuronalProblem, args...; kwargs...)
    DiffEqBase.__init(prob,args...; kwargs...)
end

function DiffEqBase.__init(prob::HHCableProblem, alg::BackwardsEuler; dt = error("dt is required for this algorithm"))
    cable = prob.geometry
    Hines = HinesMatrix(cable, alg, dt)
    hhcable_init(prob.f,
                 DiffEqBase.isinplace(prob),
                 prob.u0,
                 prob.tspan[1],
                 dt,
                 prob.p,
                 Hines)
end

function DiffEqBase.__solve(prob::HHCableProblem, alg::BackwardsEuler; dt = error("dt is required for this algorithm"))
    u0    = prob.u0
    tspan = prob.tspan
    ts    = Array(tspan[1]:dt:tspan[2])
    n     = length(ts)
    us    = Vector{typeof(u0)}(undef, n)
    cable = prob.geometry
    # With fixed-step methods we can pre-allocate the history of results
    us[1] = copy(u0)
    for j in 2:n
        us[j] = similar(u0)
    end
    integ = DiffEqBase.__init(prob, alg, dt = dt)
    for i in 2:n
        step!(integ)
        us[i] .= integ.u
    end
    # for now just return the vector series of Vm values
    # sol = DiffEqBase.build_solution(prob, alg, ts, us, calculate_error = false)
    # return sol
    return us
end

@inline function hhcable_init(f::F, IIP::Bool, u0::S, t0::T, dt::T, p::P, H::HinesMatrix) where
    {F, P, T, S<:AbstractArray{T}}
    integ = BEHHCable{IIP, S, T, P, F, typeof(H)}(f, copy(u0), copy(u0), t0, t0, t0, dt, p, H)
    return integ
end

# Solves an unbranched cable (tridiagonal matrix); needs modification for generalization to branched cable case
# In branched structures, references to the parent compartment (idx - 1) are substituted with a lookup table
# NOTE: This solver is independent of the DE method (its just slightly-modified LAPACK--see `dgtsv()`)
function thomas!(Hines, V)
    N = length(V)
    L = Hines.L
    D = Hines.D
    U = Hines.U
    #Backward pass
    for i in N:-1:2
        f = U[i] / D[i]
        D[i-1] -= f * L[i]
        V[i-1] -= f * V[i]
    end
    # Root solve
    V[1] /= D[1]
    # Forward pass
    for j in 2:N
        V[j] -= L[j]*V[j-1]
        V[j] /= D[j]
    end
    return nothing
end

# Stepping
# In-place stepping
@inline function DiffEqBase.step!(integ::BEHHCable{true, S, T}) where {T, S}
    integ.uprev        .= integ.u
    f                   = integ.f
    p                   = integ.p
    t                   = integ.t
    dt                  = integ.dt
    update_hines!       = f.f1
    solve_gates!        = f.f2
    Hines               = integ.H
    integ.tprev = t
    integ.t += dt
    # Update coefficients of Hines matrix
    update_hines!(integ.H, integ.u, p, dt)
    # Solve Hines for membrane voltage
    thomas!(integ.H, view(integ.u, :, 1))
    # Update dimensionless channel gating vars
    solve_gates!(integ.u, p, dt)
    return nothing
end

# TODO: Interpolation
function (integ::BEHHCable)(t::T) where T
end

end # module

max_stars_repo_name: wsphillips/FutureNeuralSimulation

max_issues_repo_name: wsphillips/FutureNeuralSimulation

max_forks_repo_name: wsphillips/FutureNeuralSimulation


text:
using WORLD
using Compat
using Base.Test

@show WORLD.version

# Check consistency of the results between WORLD and WORLD.jl.
# Due to the results of the WORLD were dumped (see ./data) on linux,
# we test consistency only on linux.
@static is_linux() ? include("consistency.jl") : nothing

function test_dio(x, fs::Int=44100, period::Float64=5.0)
    println("test_dio: fs=$(fs), period=$(period)")
    opt = DioOption(71.0, 800.0, 2, period, 1)
    f0, timeaxis = dio(x, fs, opt)
    @test all(isfinite(f0))
    @test all(f0 .>= 0.0)
    @test all(isfinite(timeaxis))
end

function test_stonemask(x, fs::Int=44100, period::Float64=5.0)
    println("test_stonemask: fs=$(fs), period=$(period)")
    opt = DioOption(71.0, 800.0, 2, period, 1)
    f0, timeaxis = dio(x, fs, opt)
    f0 = stonemask(x, fs, timeaxis, f0)
    @test all(isfinite(f0))
    @test all(f0 .>= 0.0)
end

function test_cheaptrick(x, fs::Int=44100, period::Float64=5.0)
    println("test_cheaptrick: fs=$(fs), period=$(period)")
    opt = DioOption(71.0, 800.0, 2, period, 1)
    f0, timeaxis = dio(x, fs, opt)
    f0 = stonemask(x, fs, timeaxis, f0)
    spectrogram = cheaptrick(x, fs, timeaxis, f0)
    @test all(isfinite(spectrogram))
end

function test_d4c(x, fs::Int=44100, period::Float64=5.0)
    println("test_d4c: fs=$(fs), period=$(period)")
    opt = DioOption(71.0, 800.0, 2, period, 1)
    f0, timeaxis = dio(x, fs, opt)
    f0 = stonemask(x, fs, timeaxis, f0)
    aperiodicity = d4c(x, fs, timeaxis, f0)
    @test all(isfinite(aperiodicity))
end

# speech -> {f0, envelope, aperiodicity} -> speech
function test_synthesis(x::AbstractArray, fs::Int=44100,
                        period::Float64=5.0,
                        tol::Float64=0.1)
    println("test_synthesis: fs=$(fs), period=$(period), tol=$(tol)")

    opt = DioOption(71.0, 800.0, 2, period, 1)

    # Fundamental frequency (f0) estimation by DIO
    f0, timeaxis = dio(x, fs, opt)
    @test all(isfinite(f0))
    @test all(isfinite(timeaxis))

    # F0 re-estimation by StoneMask
    f0 = stonemask(x, fs, timeaxis, f0)
    @test all(isfinite(f0))

    # Spectral envelope estimation
    spectrogram = cheaptrick(x, fs, timeaxis, f0)
    @test all(isfinite(spectrogram))

    # Aperiodicity ratio estimation by D4C
    aperiodicity = d4c(x, fs, timeaxis, f0)
    @test all(isfinite(aperiodicity))

    # Sysnthesis from f0, spectral envelope and aperiodicity ratio.
    y_length = convert(Int, (length(f0)-1)*period/1000 * fs + 1)
    y = synthesis(f0, spectrogram, aperiodicity, period, fs, length(x))
    @test all(isfinite(y))

    minlen = min(length(x), length(y))
    errorrate = mean(abs(y[1:minlen]-x[1:minlen])) / maxabs(x[1:minlen])

    println("errorrate=$(errorrate)")

    @test errorrate < tol
end

# Read test data (male)
x = vec(readdlm(joinpath(dirname(@__FILE__), "data", "x.txt")))
fs = 16000

for period in [5.0, 7.0, 10.0]
    test_dio(x, fs, period)
    test_stonemask(x, fs, period)
    test_cheaptrick(x, fs, period)
    test_d4c(x, fs, period)
end

# Test WORLD speech decomposition and re-synthesis
for (period, tol) in ([5.0, 0.15],[10.0, 0.16])
    test_synthesis(x, fs, period, tol)
end

println("WORLD decomposition and re-synthesis tests passed.")

# spectrum envelop <-> mel-cepstrum conversion

let
    @assert fs == 16000
    f0, timeaxis = dio(x, fs)
    f0 = stonemask(x, fs, timeaxis, f0)
    spectrogram = cheaptrick(x, fs, timeaxis, f0)
    spec = spectrogram[:,30]

    α = 0.41
    fftlen = get_fftsize_for_cheaptrick(fs)

    # check normalized mean squared error
    approximate_spec = mc2sp(sp2mc(spec, 25, α), α, fftlen)
    nmse25 = norm(log(spec) - log(approximate_spec))/norm(log(spec))
    @test nmse25 <= 0.06

    approximate_spec = mc2sp(sp2mc(spec, 30, α), α, fftlen)
    nmse30 = norm(log(spec) - log(approximate_spec))/norm(log(spec))
    @test nmse30 <= 0.05

    approximate_spec = mc2sp(sp2mc(spec, 40, α), α, fftlen)
    nmse40 = norm(log(spec) - log(approximate_spec))/norm(log(spec))
    @test nmse40 <= 0.03

    @test nmse25 > nmse30 > nmse40

    # For matrix input
    spec_mat = spectrogram[:,30:31]
    approximate_spec_mat = mc2sp(sp2mc(spec_mat, 25, α), α, fftlen)
    @test_approx_eq approximate_spec_mat[:,1] mc2sp(sp2mc(spec_mat[:,1], 25, α), α, fftlen)
    @test_approx_eq approximate_spec_mat[:,2] mc2sp(sp2mc(spec_mat[:,2], 25, α), α, fftlen)
end

# get_fftsize

for fs in [16000, 20000]
    c = get_fftsize_for_cheaptrick(fs)
    @test c == 1024
end

for fs in [44100, 48000]
    c = get_fftsize_for_cheaptrick(fs)
    @test c == 2048
end

# DioOption

let
    try DioOption(); catch @test false; end

    # keyword arguments
    try DioOption(f0floor=60.0); catch @test false; end
    try DioOption(f0ceil=600.0); catch @test false; end
    try DioOption(channels_in_octave=2.0); catch @test false; end
    try DioOption(period=10.0); catch @test false; end
    try DioOption(speed=5); catch @test false; end

    # f0loor
    @test_throws ArgumentError DioOption(-1,  100, 2.0, 5.0, 11)
    @test_throws ArgumentError DioOption(200, 100, 2.0, 5.0, 11)
    # channels in octave
    @test_throws ArgumentError DioOption(80.0, 640.0, -1.0, 5.0, 11)
    # period
    @test_throws ArgumentError DioOption(80.0, 640.0, 2.0, -1.0, 11)
    @test_throws ArgumentError DioOption(80.0, 640.0, 2.0, 0.0, 11)
    # speed
    try DioOption(80.0, 640.0, 2.0, 5.0, 1); catch @test false; end
    try DioOption(80.0, 640.0, 2.0, 5.0, 12); catch @test false; end
    @test_throws ArgumentError DioOption(80.0, 640.0, 2.0, 5.0, 0)
    @test_throws ArgumentError DioOption(80.0, 640.0, 2.0, 5.0, 13)
    # allowed_range
    if WORLD.version >= v"0.2.1-2"
        try DioOption(allowed_range=0.1); catch @test false; end
        try DioOption(71.0, 800.0, 2.0, 5.0, 1, 0.0); catch @test false; end
        @test_throws ArgumentError DioOption(80.0, 640.0, 2.0, 5.0, 0, -1.0)
    end
end

# matlabfunctions

let
    opt = DioOption(71.0, 800.0, 2, 5.0, 1)
    f0, timeaxis = dio(x, fs, opt)
    spectrogram = cheaptrick(x, fs, timeaxis, f0)

    logspec = log(spectrogram[:,100])
    fftsize = get_fftsize_for_cheaptrick(fs)

    freqaxis_src = collect(1:(fftsize>>1+1)) / fftsize * fs
    freqaxis_dst = freqaxis_src .* 2/3

    # spectral stretching
    interpolated_logspec = interp1(freqaxis_src, logspec, freqaxis_dst)
    @test all(isfinite(interpolated_logspec))

    interpolated_spec = exp(interpolated_logspec)

    interpolated_logspec2 = similar(logspec)
    interp1!(freqaxis_src, logspec, freqaxis_dst, interpolated_logspec2)
    @test interpolated_logspec == interpolated_logspec2
end


text:
using ChemicalIdentifiers
using Test

@testset "ChemicalIdentifiers.jl" begin
    res1 = search_chemical("water",nothing)
    res2 =  search_chemical("SMILES=O",nothing)
    @test res1.formula == res2.formula
end


text:
using Test

# eg
#
# $ julia --project
#
# Then
#
# julia> using Pkg; Pkg.test()

icy_antics() = "Is a sea ice modeling framework"
@test icy_antics() == "Is a sea ice modeling framework"


text:
"""
    Grid

Gride layout.
"""
struct Grid end

(::Grid)(g) = Point.(5 .* (vertices(g) .- 1), 0)

is_self_loop(e::AbstractEdge) = src(e) == dst(e)
any_self_loops(g::AbstractGraph) = any(is_self_loop, edges(g))


text:
function camelcase(cur_string::AbstractString)
  lcfirst(pascalcase(cur_string))
end


text:
"""
necessary functions for outer approximation framework
"""

"""
function init_oad(op::OriginalProblem, oad::OAdata)
    initiates OAdata structure using information of the OriginalProblem.
"""
function init_oad(op::OriginalProblem, oad::OAdata)
    oad.num_l_constr              = op.num_l_constr
    oad.ref_l_var                 = Float64[]
    oad.ref_u_var                 = Float64[]
    oad.ref_num_var               = 0
    oad.ref_num_nl_constr         = 0

    oad.obj_sense                 = op.obj_sense

    oad.oa_status                 = :Unknown
    oad.oa_started                = false
    oad.incumbent                 = Float64[]
    oad.new_incumbent             = false
    oad.total_time                = 0
    oad.obj_val                   = Inf
    oad.obj_bound                 = -Inf
    oad.obj_gap                   = Inf
    oad.oa_iter                   = 0

    oad.milp_sol_available        = false
    oad.ref_mip_solution          = Float64[]
    oad.mip_infeasible            = false
    oad.ref_nlp_solution          = Float64[]
    oad.nlp_infeasible            = false
    oad.ref_feasibility_solution  = Float64[]

    oad.int_idx                   = filter(i -> (op.var_type[i] in (:Int, :Bin)), 1:op.num_var)
    oad.prev_ref_mip_solution     = Float64[]

    oad.ref_linear_le_constraints = []
    oad.ref_linear_ge_constraints = []
    oad.ref_linear_eq_constraints = []
end

"""
function construct_linear_model(optimizer::MOI.AbstractOptimizer, op::OriginalProblem, oad::OAdata)
    this function generates a new mip_model with the linear constraints and linear objective function
    of the original model OriginalProblem.
"""
function construct_linear_model(optimizer::MOI.AbstractOptimizer, op::OriginalProblem, oad::OAdata)
    oad.mip_model = Model()
    lb = op.l_var
    ub = op.u_var
    push!(lb, -Inf)
    push!(ub, Inf)
    #vartype = op.var_type
    #push!(vartype, :Cont)
    # all continuous we solve relaxation first
    @variable(oad.mip_model, lb[i] <= x[i=1:(op.num_var+1)] <= ub[i])
    for i in 1:length(op.var_type)
        op.var_type[i] == :Int && JuMP.set_integer(x[i])
        op.var_type[i] == :Bin && JuMP.set_binary(x[i])
    end
    # register function into the mip model
    register_functions!(oad.mip_model, op.options.registered_functions)
    # attach linear constraints to the mip model
    backend1 = JuMP.backend(oad.mip_model);
    llc = optimizer.linear_le_constraints
    lgc = optimizer.linear_ge_constraints
    lec = optimizer.linear_eq_constraints
    #add linear constraints to the mip model
    for constr_type in [llc, lgc, lec]
        for constr in constr_type
            MOI.add_constraint(backend1, constr[1], constr[2])
        end
    end
    #set objective
    if op.has_nl_objective
        @info "the objective function is nonlinear.\n"
        op.obj_sense == :Min && @objective(oad.mip_model, Min, x[op.num_var+1])
        op.obj_sense == :Max && @objective(oad.mip_model, Max, x[op.num_var+1])
    else
        @info "the objective function is linear or quadratic.\n"
        MOI.set(oad.mip_model, MOI.ObjectiveFunction{typeof(optimizer.objective)}(), optimizer.objective)
        MOI.set(oad.mip_model, MOI.ObjectiveSense(), optimizer.sense)
    end
    oad.mip_x = x
end

"""
function add_oa_cut(model::MOI.AbstractOptimizer, op::OriginalProblem, oad::OAdata)
    this function generates oa cuts for nonlinear constraints and objective function and adds the cuts
    to the mip_model which is a milp model.
"""
function add_oa_cut(model::MOI.AbstractOptimizer, op::OriginalProblem, oad::OAdata)
    jac_IJ = MOI.jacobian_structure(model.nlp_data.evaluator)
    nlp_solution=ones(op.num_var)

    g_val = zeros(op.num_nl_constr)
    g_jac = zeros(length(jac_IJ))
    MOI.eval_constraint(model.nlp_data.evaluator, g_val, nlp_solution)
    MOI.eval_constraint_jacobian(model.nlp_data.evaluator, g_jac, nlp_solution)

    #constraints cut
    varidx_new = [zeros(Int, 0) for i in 1:op.num_nl_constr]
    coef_new = [zeros(0) for i in 1:op.num_nl_constr]

    for k in 1:length(jac_IJ)
        row = jac_IJ[k][1]
        push!(varidx_new[row], jac_IJ[k][2])
        push!(coef_new[row], g_jac[k])
    end

    for i in 1:op.num_nl_constr
        new_rhs = -g_val[i]
        for j in 1:length(varidx_new[i])
            new_rhs += coef_new[i][j] * nlp_solution[Int(varidx_new[i][j])]
        end
        if op.nlp_constr_type[i] == :(<=)
            @constraint(oad.mip_model, dot(coef_new[i], oad.mip_x[varidx_new[i]]) <= new_rhs)
        else
            @constraint(oad.mip_model, dot(coef_new[i], oad.mip_x[varidx_new[i]]) >= new_rhs)
        end
    end

    #objective cut
    if op.has_nl_objective
        f_nab = zeros(op.num_var+1)
        f_val = MOI.eval_objective(model.nlp_data.evaluator, nlp_solution)
        MOI.eval_objective_gradient(model.nlp_data.evaluator, f_nab, nlp_solution)

        new_rhs = -f_val
        varidx = zeros(Int, op.num_var+1)
        for j in 1:op.num_var
            varidx[j] = j
            new_rhs += f_nab[j] * nlp_solution[j]
        end

        varidx[op.num_var+1] = op.num_var + 1
        f_nab[op.num_var+1] = -1.0

        if op.obj_sense == :Max
            @constraint(oad.mip_model, dot(f_nab, oad.mip_x[varidx]) >= new_rhs)
        else
            @constraint(oad.mip_model, dot(f_nab, oad.mip_x[varidx]) <= new_rhs)
        end
    end
end

"""
function construct_nlp_model(model::MOI.AbstractOptimizer, op::OriginalProblem, oad::OAdata)
    this function generates a nlp_model from the original model provided as OriginalProblem.
"""
function construct_nlp_model(model::MOI.AbstractOptimizer, op::OriginalProblem, oad::OAdata)
    oad.nlp_model = Model()
    lb = op.l_var
    ub = op.u_var
    @variable(oad.nlp_model, lb[i] <= x[i=1:op.num_var] <= ub[i])
    register_functions!(oad.nlp_model, op.options.registered_functions)

    backend2 = JuMP.backend(oad.nlp_model);
    llc = model.linear_le_constraints
    lgc = model.linear_ge_constraints
    lec = model.linear_eq_constraints

    for constr_type in [llc, lgc, lec]
        for constr in constr_type
            MOI.add_constraint(backend2, constr[1], constr[2])
        end
    end

    for i in 1:op.num_nl_constr
        constr_expr = MOI.constraint_expr(model.nlp_data.evaluator, i)
        constr_expr = expr_dereferencing(constr_expr, oad.nlp_model)
        try
            JuMP.add_NL_constraint(oad.nlp_model, constr_expr)
        catch
            error("Have you registered a function? Then please register the function.\n")
        end
    end

    if model.nlp_data.has_objective
        obj_expr = MOI.objective_expr(model.nlp_data.evaluator)
        obj_expr = expr_dereferencing(obj_expr, oad.nlp_model)
        try
            JuMP.set_NL_objective(oad.nlp_model, model.sense, obj_expr)
        catch
            error("Have you registered a function? Then please register the function.\n")
        end
    elseif model.objective !== nothing
        MOI.set(oad.nlp_model, MOI.ObjectiveFunction{typeof(model.objective)}(), model.objective)
        MOI.set(oad.nlp_model, MOI.ObjectiveSense(), model.sense)
    end

    oad.nlp_x = x

    return nothing
end

"""
function solve_mip_model(model::MOI.AbstractOptimizer, op::OriginalProblem, oad::OAdata)
    this function solves mip_model with mixed interger solver.
"""

function solve_mip_model(model::MOI.AbstractOptimizer, op::OriginalProblem, oad::OAdata)
    JuMP.set_optimizer(oad.mip_model, model.options.mip_solver)
    JuMP.optimize!(oad.mip_model)

    if termination_status(oad.mip_model) == :OPTIMAL || termination_status(oad.mip_model) == :LOCALLY_SOLVED
        @info "There is a feasible solution in the mip_model.\n"
        value = JuMP.value.(oad.mip_x)
        return value
    else
        @error "There is no feasible solution from mip_model. \n"
    end
    #status = JuMP.primal_status(oad.solve_mip_model)
end

function solve_nlp_model(model::MOI.AbstractOptimizer, op::OriginalProblem, oad::OAdata)
    JuMP.set_optimizer(oad.nlp_model, model.options.nlp_solver)
    JuMP.optimize!(oad.nlp_model)

    if termination_status(oad.nlp_model) == :OPTIMAL || termination_status(oad.nlp_model) == :LOCALLY_SOLVED
        @info "There is a feasible solution in the mip_model.\n"
        value = JuMP.value.(oad.nlp_x)
        return value
    else
        @error "There is no feasible solution from mip_model. \n"
    end
    #status = JuMP.primal_status(oad.solve_mip_model)
end

"""
function construct_ref_nlp_model(model::MOI.AbstractOptimizer, op::OriginalProblem, oad::OAdata)
    this function generates a ref_nlp_model from the original model.
"""
function construct_ref_nlp_model(model::MOI.AbstractOptimizer, op::OriginalProblem, oad::OAdata)
    oad.ref_nlp_model = Model()
    lb = op.l_var
    ub = op.u_var
    vartype = op.var_type

    for i in 1:(op.nbinvars)
        push!(lb, -Inf)
        push!(ub, Inf)
        push!(vartype, :Cont)
    end

    oad.ref_l_var = lb
    oad.ref_u_var = ub
    oad.ref_var_type = vartype

    @variable(oad.ref_nlp_model, lb[i] <= x[i=1:(op.num_var+op.nbinvars)] <= ub[i])
    register_functions!(oad.ref_nlp_model, op.options.registered_functions)

    oad.ref_num_var = length(x)

    backend = JuMP.backend(oad.ref_nlp_model);
    llc = model.linear_le_constraints
    lgc = model.linear_ge_constraints
    lec = model.linear_eq_constraints

    for constr_type in [llc, lgc, lec]
        for constr in constr_type
            MOI.add_constraint(backend, constr[1], constr[2])
        end
    end

    for i in 1:op.num_nl_constr
        constr_expr = MOI.constraint_expr(model.nlp_data.evaluator, i)
        constr_expr = expr_dereferencing(constr_expr, oad.ref_nlp_model)
        try
            JuMP.add_NL_constraint(oad.ref_nlp_model, constr_expr)
        catch
            error("Have you registered a function? Then please register the function.\n")
        end
    end
    num_nl_constr = op.num_nl_constr
    #additional constraints from the nlobj
    if model.nlp_data.has_objective
        if op.obj_sense == :Max
            @info "The objective sense is Max. \n"
            @objective(oad.ref_nlp_model, Max, sum(x[i] for i in (op.num_var+1):(op.num_var+op.nbinvars)))
        elseif op.obj_sense == :Min
            @info "The objective sense is Min. \n"
            @objective(oad.ref_nlp_model, Min, sum(x[i] for i in (op.num_var+1):(op.num_var+op.nbinvars)))
        end

        obj_expr = MOI.objective_expr(model.nlp_data.evaluator)
        moi_bin_vars, moi_obj_terms = get_binvar_and_objterms(obj_expr, oad.ref_nlp_model, op)

        siz = length(moi_bin_vars)
        jump_bin_vars = []
        jump_obj_terms = []
        for i in 1:siz
            jump_bin = JuMP.VariableRef(oad.ref_nlp_model, moi_bin_vars[i].args[2])
            push!(jump_bin_vars, jump_bin)
            jump_exp = expr_dereferencing(moi_obj_terms[i], oad.ref_nlp_model)
            push!(jump_obj_terms, jump_exp)
        end
        num_nl_constr = num_nl_constr+siz
        ref_nlpcontr_type = op.nlp_constr_type
        for i in 1:siz
            L=-100000
            U=+100000
            JuMP.@constraint(oad.ref_nlp_model, x[i+op.num_var]-jump_bin_vars[i]*L >= 0)
            exp = :($(x[i+op.num_var])-($(jump_bin_vars[i])-$(1))*$(U)-$(jump_obj_terms[i]) >= 0)
            JuMP.add_NL_constraint(oad.ref_nlp_model, exp)
            push!(ref_nlpcontr_type, :>=)
        end
        oad.ref_nlp_constr_type = ref_nlpcontr_type

    elseif model.objective !== nothing
        MOI.set(oad.ref_nlp_model, MOI.ObjectiveFunction{typeof(model.objective)}(), model.objective)
        MOI.set(oad.ref_nlp_model, MOI.ObjectiveSense(), model.sense)
    end
    oad.ref_num_nl_constr = num_nl_constr
    oad.ref_nlp_x = x

    return nothing
end

function _expr_deref!(expr, exprd, m, op, dict, bin_vars, bin_obj_terms)
    for i in 2:length(expr.args)
        if isa(expr.args[i], Union{Float64,Int64,JuMP.VariableRef})
            k = 0
        elseif expr.args[i].head == :ref
            @assert isa(expr.args[i].args[2], MOI.VariableIndex)
            exprd.args[i] = JuMP.VariableRef(m, expr.args[i].args[2])
            if op.var_type[dict[exprd.args[i]]] == :Bin
                push!(bin_vars, expr.args[i])
                push!(bin_obj_terms, expr.args[i+1])   #assuming bin variable appears first
            end
        elseif expr.args[i].head == :call
            _expr_deref!(expr.args[i], exprd.args[i], m, op, dict, bin_vars, bin_obj_terms)
        else
            error("expr_dereferencing :: Unexpected term in expression tree.")
        end
    end
    return bin_vars, bin_obj_terms
end

function get_binvar_and_objterms(obj_expr, m, op)
    list = JuMP.all_variables(m)
    dict = Dict{}(list[i] => i for i in 1:op.num_var)
    bin_vars=[]
    bin_obj_terms=[]
    c_expr = copy(obj_expr)
    d_expr = copy(obj_expr)
    return _expr_deref!(c_expr, d_expr, m, op, dict, bin_vars, bin_obj_terms)
end

"""
function construct_ref_mip_model(model::MOI.AbstractOptimizer, op::OriginalProblem, oad::OAdata)
    This function genearates reformulated mip model.
"""
function construct_ref_mip_model(model::MOI.AbstractOptimizer, op::OriginalProblem, oad::OAdata)
    oad.ref_mip_model = Model()
    lb = oad.ref_l_var
    ub = oad.ref_u_var
    # all continuous we solve relaxation first
    @variable(oad.ref_mip_model, lb[i] <= x[i=1:oad.ref_num_var] <= ub[i], binary =oad.ref_var_type[i] == :Bin, integer = oad.ref_var_type[i] == :Int)
    # register function into the mip model
    register_functions!(oad.ref_mip_model, op.options.registered_functions)
    # attach linear constraints to the mip model
    backend = JuMP.backend(oad.ref_mip_model);
    llc = oad.ref_linear_le_constraints
    lgc = oad.ref_linear_ge_constraints
    lec = oad.ref_linear_eq_constraints
    for constr_type in [llc, lgc, lec]
        for constr in constr_type
            MOI.add_constraint(backend, constr[1], constr[2])
        end
    end
    #set var type
    for i in 1:length(oad.ref_var_type)
        oad.ref_var_type[i] == :Int && JuMP.set_integer(x[i])
        oad.ref_var_type[i] == :Bin && JuMP.set_binary(x[i])
    end
    #set objective
    func = JuMP.objective_function(oad.ref_nlp_model)
    func = JuMP.moi_function(func)
    MOI.set(oad.ref_mip_model, MOI.ObjectiveFunction{typeof(func)}(), func)
    MOI.set(oad.ref_mip_model, MOI.ObjectiveSense(), model.sense)
    oad.ref_mip_x = x
    oad.prev_ref_mip_solution = rand(oad.ref_num_var)
end

"""
function add_ref_oa_cut(model::MOI.AbstractOptimizer, op::OriginalProblem, oad::OAdata)
    This function facilates to add linear approximation cut to the reformulated mip model.
"""
function add_ref_oa_cut(model::MOI.AbstractOptimizer, op::OriginalProblem, oad::OAdata)
    d=JuMP.NLPEvaluator(oad.ref_nlp_model)
    features = MOI.features_available(d)
    MOI.initialize(d, features)

    jac_IJ = MOI.jacobian_structure(d)

    !(oad.nlp_infeasible) && (nlp_solution = oad.ref_nlp_solution)
    oad.nlp_infeasible && (nlp_solution = oad.ref_feasibility_solution[1:oad.ref_num_var])

    g_val = zeros(oad.ref_num_nl_constr)
    g_jac = zeros(length(jac_IJ))
    MOI.eval_constraint(d, g_val, nlp_solution)
    MOI.eval_constraint_jacobian(d, g_jac, nlp_solution)

    #constraints cut
    varidx_new = [zeros(Int, 0) for i in 1:oad.ref_num_nl_constr]
    coef_new = [zeros(0) for i in 1:oad.ref_num_nl_constr]

    for k in 1:length(jac_IJ)
        row = jac_IJ[k][1]
        push!(varidx_new[row], jac_IJ[k][2])
        push!(coef_new[row], g_jac[k])
    end

    for i in 1:oad.ref_num_nl_constr
        new_rhs = -g_val[i]
        for j in 1:length(varidx_new[i])
            new_rhs += coef_new[i][j] * nlp_solution[Int(varidx_new[i][j])]
        end
        if oad.ref_nlp_constr_type[i] == :(<=)
            @constraint(oad.ref_mip_model, dot(coef_new[i], oad.ref_mip_x[varidx_new[i]]) <= new_rhs)
        else
            @constraint(oad.ref_mip_model, dot(coef_new[i], oad.ref_mip_x[varidx_new[i]]) >= new_rhs)
        end
    end
    #no need to derive obj cut as the ref_nlp obj is linear.
end

"""
function construct_ref_feasibility_model(oad::OAdata)
    This function generates a feasibility model in case the nlp problem is infeasible.
"""
function construct_ref_feasibility_model(oad::OAdata)
    # in case of infeasible nlp problem, we need to solve a feasibility problem
    oad.ref_feasibility_model = Model()
    @variable(oad.ref_feasibility_model,  x[1:(oad.ref_num_var+oad.ref_num_nl_constr)] >= 0.0)
    #linear constraints are automatically satisfied
    #we only need to focus on nonlinear constraints
    d = JuMP.NLPEvaluator(oad.ref_nlp_model)
    features = MOI.features_available(d)
    MOI.initialize(d, features)
    #define constraints
    for i in 1:jp.ref_num_nl_constr
        constr_expr = MOI.constraint_expr(d, i)
        con_expr = expr_dereferencing(constr_expr.args[2], oad.ref_feasibility_model)
        #print($con_expr, "\n")
        if constr_expr.args[1] == :(<=)
            con_ex = :($(con_expr) <= $(x[i+oad.ref_num_var]))
            JuMP.add_NL_constraint(oad.ref_feasibility_model, con_ex)
        elseif constr_expr.args[1] == :(>=)
            con_ex = :($(con_expr) >= $(x[i+oad.ref_num_var]))
            JuMP.add_NL_constraint(oad.ref_feasibility_model, con_ex)
        elseif constr_expr.args[1] == :(==)
            @error "\n Nonlinear constarint should not be equality constarint. \n"
        end
    end
    #define obj function
    oad.obj_sense == :Min && JuMP.@objective(oad.ref_feasibility_model, Min, sum(x[i] for i in (oad.ref_num_var+1):(oad.ref_num_var+oad.ref_num_nl_constr)))
    oad.obj_sense == :Max && JuMP.@objective(oad.ref_feasibility_model, Max, sum(x[i] for i in (oad.ref_num_var+1):(oad.ref_num_var+oad.ref_num_nl_constr)))

    oad.ref_feasibility_x = x
end

"""
function reserve_ref_linear_constraints(oad::OAdata)
    This function saves linear constraints of the reformulated nlp model to the outer approximation data structure.
    This constraints are in MOI format and later will be used to generate reformulated mip model.
"""
function reserve_ref_linear_constraints(oad::OAdata)
    #list of different types of constraints
    constr_types = JuMP.list_of_constraint_types(oad.ref_nlp_model)
    #number of different constraint types
    siz = length(constr_types)
    #reserve constraints for each type
    for i in 1:siz
        num_constr = JuMP.num_constraints(oad.ref_nlp_model, constr_types[i][1], constr_types[i][2])
        constrs    = JuMP.all_constraints(oad.ref_nlp_model, constr_types[i][1], constr_types[i][2])
        for j in 1:num_constr
            con_obj = JuMP.constraint_object(constrs[j])
            func    = JuMP.moi_function(con_obj.func)
            set     = con_obj.set
            (constr_types[i][2] == MOI.EqualTo{Float64}) && push!(oad.ref_linear_eq_constraints, (func, set))
            (constr_types[i][2] == MOI.LessThan{Float64}) && push!(oad.ref_linear_le_constraints, (func, set))
            (constr_types[i][2] == MOI.GreaterThan{Float64}) && push!(oad.ref_linear_ge_constraints, (func, set))
        end
    end
end

max_stars_repo_name: MdMahbubarRahman/Duran

max_issues_repo_name: MdMahbubarRahman/Duran

max_forks_repo_name: MdMahbubarRahman/Duran


text:
import .ValueHistories: MVHistory

function Base.convert(::Type{MVHistory}, tbl::TBReadable; kwargs...)
    hist = MVHistory()

    map_summaries(tbl; kwargs...) do tag, iter, val
        push!(hist, Symbol(tag), iter, val)
    end

    return hist
end

function Base.getindex(mvhist::MVHistory, tag::AbstractString)
    return getindex(mvhist, Symbol(tag))
end

ValueHistories.MVHistory(tbl::TBReadable; kwargs...) = convert(MVHistory, tbl; kwargs...)


text:
# Filter
b = [1., 2., 3., 4.]
x = [1., 1., 0., 1., 1., 0., 0., 0.]
@test filt(b, 1., x)  == [1., 3., 5., 8., 7., 5., 7., 4.]
@test filt(b, [1., -0.5], x)  == [1., 3.5, 6.75, 11.375, 12.6875, 11.34375, 12.671875, 10.3359375]

# Convolution
a = [1., 2., 1., 2.]
b = [1., 2., 3.]
@test_approx_eq conv(a, b) [1., 4., 8., 10., 7., 6.]
@test_approx_eq conv(complex(a, ones(4)), complex(b)) complex([1., 4., 8., 10., 7., 6.], [1., 3., 6., 6., 5., 3.])

# Discrete cosine transform (DCT) tests

a = rand(8,11) + im*rand(8,11)
@test norm(idct(dct(a)) - a) < 1e-8

X = reshape([1,2,7,2,1,5,9,-1,3,4,6,9],3,4)
Y = rand(17,14)
Y[3:5,9:12] = X
sX = slice(Y,3:5,9:12)

true_Xdct = [  13.856406460551018  -3.863239728836245   2.886751345948129  -0.274551994240164; -2.828427124746190  -2.184015211898548  -4.949747468305834   3.966116180118245; 4.898979485566356  -0.194137576915510  -2.857738033247041   2.731723009609389 ]

true_Xdct_1 = [    5.773502691896258   4.618802153517007   6.350852961085884  10.969655114602890; -4.242640687119286  -2.121320343559643   4.242640687119286  -3.535533905932738; 1.632993161855452   2.041241452319315   5.715476066494083   0.408248290463863 ]

true_Xdct_2 = [    8.  -3.854030797826254  -3.0  3.761176226848022;
   4.0  -2.071929829606556   4.0  -2.388955165168770; 12.  -0.765366864730179   4.0  -1.847759065022573 ]

Xdct = dct(X)
Xdct! = float(X); dct!(Xdct!)
Xdct_1 = dct(X,1)
Xdct!_1 = float(X); dct!(Xdct!_1,1)
Xdct_2 = dct(X,2)
Xdct!_2 = float(X); dct!(Xdct!_2,2)

Xidct = idct(true_Xdct)
Xidct! = copy(true_Xdct); idct!(Xidct!)
Xidct_1 = idct(true_Xdct_1,1)
Xidct!_1 = copy(true_Xdct_1); idct!(Xidct!_1,1)
Xidct_2 = idct(true_Xdct_2,2)
Xidct!_2 = copy(true_Xdct_2); idct!(Xidct!_2,2)

pXdct = plan_dct(X)(X)
pXdct! = float(X); plan_dct!(pXdct!)(pXdct!)
pXdct_1 = plan_dct(X,1)(X)
pXdct!_1 = float(X); plan_dct!(pXdct!_1,1)(pXdct!_1)
pXdct_2 = plan_dct(X,2)(X)
pXdct!_2 = float(X); plan_dct!(pXdct!_2,2)(pXdct!_2)

pXidct = plan_idct(true_Xdct)(true_Xdct)
pXidct! = copy(true_Xdct); plan_idct!(pXidct!)(pXidct!)
pXidct_1 = plan_idct(true_Xdct_1,1)(true_Xdct_1)
pXidct!_1 = copy(true_Xdct_1); plan_idct!(pXidct!_1,1)(pXidct!_1)
pXidct_2 = plan_idct(true_Xdct_2,2)(true_Xdct_2)
pXidct!_2 = copy(true_Xdct_2); plan_idct!(pXidct!_2,2)(pXidct!_2)

sXdct = dct(sX)
psXdct = plan_dct(sX)(sX)
sYdct! = copy(Y); sXdct! = slice(sYdct!,3:5,9:12); dct!(sXdct!)
psYdct! = copy(Y); psXdct! = slice(psYdct!,3:5,9:12); plan_dct!(psXdct!)(psXdct!)

for i = 1:length(X)
    @test_approx_eq Xdct[i] true_Xdct[i]
    @test_approx_eq Xdct![i] true_Xdct[i]
    @test_approx_eq Xdct_1[i] true_Xdct_1[i]
    @test_approx_eq Xdct!_1[i] true_Xdct_1[i]
    @test_approx_eq Xdct_2[i] true_Xdct_2[i]
    @test_approx_eq Xdct!_2[i] true_Xdct_2[i]

    @test_approx_eq pXdct[i] true_Xdct[i]
    @test_approx_eq pXdct![i] true_Xdct[i]
    @test_approx_eq pXdct_1[i] true_Xdct_1[i]
    @test_approx_eq pXdct!_1[i] true_Xdct_1[i]
    @test_approx_eq pXdct_2[i] true_Xdct_2[i]
    @test_approx_eq pXdct!_2[i] true_Xdct_2[i]

    @test_approx_eq Xidct[i] X[i]
    @test_approx_eq Xidct![i] X[i]
    @test_approx_eq Xidct_1[i] X[i]
    @test_approx_eq Xidct!_1[i] X[i]
    @test_approx_eq Xidct_2[i] X[i]
    @test_approx_eq Xidct!_2[i] X[i]

    @test_approx_eq pXidct[i] X[i]
    @test_approx_eq pXidct![i] X[i]
    @test_approx_eq pXidct_1[i] X[i]
    @test_approx_eq pXidct!_1[i] X[i]
    @test_approx_eq pXidct_2[i] X[i]
    @test_approx_eq pXidct!_2[i] X[i]

    @test_approx_eq sXdct[i] true_Xdct[i]
    @test_approx_eq psXdct[i] true_Xdct[i]
    @test_approx_eq sXdct![i] true_Xdct[i]
    @test_approx_eq psXdct![i] true_Xdct[i]
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
using Pkg
Pkg.instantiate()
#Pkg.activate(".")
using Weave

#runweave= "weave" ∈ ARGS
#runnotebook= "notebook" ∈ ARGS
runweave= true
runnotebook= false

for arg in ["jmd/rt_onlym.jmd"] #ARGS
  if !(arg ∈ ["weave", "notebook"])
    src= arg

    if runweave
      println("weaving markdown for $src")
      using Weave
      weave(src,out_path="md", cache=:user, cache_path="weavecache",  doctype="github", args=Dict("md" => true), mod=Main)
    end

    if runnotebook
      println("weaving notebook for $src")
      using Weave
      notebook(src, out_path=joinpath(pwd(),"build"), nbconvert_options="--allow-errors")
    end
  end
end


text:
using HttpServer
using Mustache
using MySQL
using JSON
using DataFrames

conn = mysql_connect("127.0.0.1", "anychart_user", "anychart_pass", "anychart_db")

function loadData()
  df = mysql_execute(conn, "SELECT name, value FROM fruits ORDER BY value DESC LIMIT 5")
  data = []
  for col in eachrow(df)
    push!(data, Dict("name" => col[1], "value" => col[2]))
  end
  data
end

http = HttpHandler() do req::Request, res::Response
  chartData = JSON.json(loadData())
  data = Dict("title" => "Anychart Julia template",
   "chartData" => chartData,
   "chartTitle" => "Top 5 fruits")
  html = render_from_file("index.html", data)
  Response(html)
end

server = Server(http)
run(server, 8000)


text:
# ---
# title: 661. Image Smoother
# id: problem661
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/image-smoother/description/>
# hidden: true
# ---
#
# Given a 2D integer matrix M representing the gray scale of an image, you need
# to design a smoother to make the gray scale of each cell becomes the average
# gray scale (rounding down) of all the 8 surrounding cells and itself. If a
# cell has less than 8 surrounding cells, then use as many as you can.
#
# **Example 1:**
#
#
#
#     Input:
#     [[1,1,1],
#      [1,0,1],
#      [1,1,1]]
#     Output:
#     [[0, 0, 0],
#      [0, 0, 0],
#      [0, 0, 0]]
#     Explanation:
#     For the point (0,0), (0,2), (2,0), (2,2): floor(3/4) = floor(0.75) = 0
#     For the point (0,1), (1,0), (1,2), (2,1): floor(5/6) = floor(0.83333333) = 0
#     For the point (1,1): floor(8/9) = floor(0.88888889) = 0
#
#
# **Note:**
#
#   1. The value in the given matrix is in the range of [0, 255].
#   2. The length and width of the given matrix are in the range of [1, 150].
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
@testset "Complex scaling" begin
    B = FEDVR(range(0,stop=20,length=71), 10)
    C = FEDVR(range(0,stop=20,length=71), 10, t₀=10.0, ϕ=π/3)
    @test C.t₀ ≥ 10.0
    @test_throws ArgumentError FEDVR(range(0,stop=20,length=71), 10, t₀=30.0, ϕ=π/3)
    @test complex_rotate(5, B) == 5
    @test complex_rotate(15, B) == 15
    @test complex_rotate(5, C) == 5
    @test complex_rotate(15, C) ≈ 10 + 5*exp(im*π/3)
end

@testset "Real locations" begin
    rₘₐₓ = 20
    R = FEDVR(range(0,stop=rₘₐₓ,length=11), 10, t₀=rₘₐₓ/2, ϕ=π/3)
    R′ = FEDVR(range(0,stop=rₘₐₓ,length=11), 10)
    @test norm(real_locs(R)-real_locs(R′)) == 0
end


text:
@testset "Selections" begin

    atoms = readPDB("./structure.pdb")

    @test length(select(atoms, "name CA")) == 104
    sel = select(atoms, "index = 13")
    @test length(sel) == 1
    @test sel[1].index == 13
    @test sel[1].index_pdb == 13

    @test length(select(atoms, "index > 1 and index < 13")) == 11

    @test length(select(atoms, "protein")) == 1463

    @test length(select(atoms, "water")) == 58014

    @test length(select(atoms, "resname GLY")) == 84

    @test length(select(atoms, "segname PROT")) == 1463

    @test length(select(atoms, "residue = 2")) == 11

    @test length(select(atoms, "neutral")) == 1182

    @test length(select(atoms, "charged")) == 281

    @test length(select(atoms, "sidechain")) == 854

    @test length(select(atoms, "acidic")) == 162

    @test length(select(atoms, "basic")) == 68

    @test length(select(atoms, "hydrophobic")) == 327

    @test length(select(atoms, "hydrophobic")) == 327

    @test length(select(atoms, "aliphatic")) == 379

    @test length(select(atoms, "aromatic")) == 344

    @test length(select(atoms, "polar")) == 880

    @test length(select(atoms, "nonpolar")) == 583

    @test maxmin(atoms, "chain A").xlength ≈ [83.083, 83.028, 82.7]

    # Test editing a field
    atoms[1].index = 0
    @test atoms[1].index == 0

    # Test residue iterator
    someresidues = select(atoms, "residue < 15")
    n = 0
    m = 0.0
    for res in eachresidue(someresidues)
        if name(res) == "SER"
            n += 1
            for atom in res
                m += mass(atom)
            end
        end
    end
    @test n == 4
    @test m ≈ 348.31340000000006

    # Residue properties (discontinous set)
    lessresidues = select(someresidues, "residue < 3 or residue > 12")
    @test residue.(eachresidue(lessresidues)) == [1, 2, 13, 14]
    @test resnum.(eachresidue(lessresidues)) == [1, 2, 13, 14]
    @test name.(eachresidue(lessresidues)) == ["ALA", "CYS", "SER", "SER"]
    @test resname.(eachresidue(lessresidues)) == ["ALA", "CYS", "SER", "SER"]
    @test chain.(eachresidue(lessresidues)) == ["A", "A", "A", "A"]
    @test model.(eachresidue(lessresidues)) == [1, 1, 1, 1]

end

max_stars_repo_name: mcubeg/PDBTools

max_issues_repo_name: mcubeg/PDBTools

max_forks_repo_name: mcubeg/PDBTools