text:
using SIAN

println("Setting up the problem")

ode = @ODEmodel(
    x1'(t) = -k01 * x1(t) + k12 * x2(t) + k13 * x3(t) + k14 * x4(t) - k21 * x1(t) - k31 * x1(t) - k41 * x1(t) + u(t),
    x2'(t) = -k12 * x2(t) + k21 * x1(t),
    x3'(t) = -k13 * x3(t) + k31 * x1(t),
    x4'(t) = -k14 * x4(t) + k41 * x1(t),
    y(t) = x1(t)
)

res = identifiability_ode(ode, get_parameters(ode); p=0.99, p_mod=2^29 - 3, nthrds=1)

println(res)

max_stars_repo_name: iliailmer/SIAN-Julia

max_issues_repo_name: iliailmer/SIAN-Julia

max_forks_repo_name: iliailmer/SIAN-Julia


text:
"""
    bounded_aux(x::Array{JuMP.VariableRef}, binary_var::JuMP.VariableRef)
    bounded_aux(x::Array{JuMP.VariableRef}, y::JuMP.VariableRef, binary_var::JuMP.VariableRef)

Generates binary-bounded auxiliary variables and their bounding constraints of the same size as x + y.
"""
function bounded_aux(x::Array{JuMP.VariableRef}, binary_var::JuMP.VariableRef)
    var_bounds = get_bounds(x)
    aux_vars = @variable(x[1].model, [1:length(x)])
    bound_cons = []
    for i = 1:length(x) # Bounding auxiliary variables
        bds = var_bounds[x[i]]
        push!(bound_cons, @constraint(x[i].model, aux_vars[i] >= minimum(bds)*binary_var))
        push!(bound_cons, @constraint(x[i].model, aux_vars[i] <= maximum(bds)*binary_var))
    end
    aux_vars = (binary_var, aux_vars)
    return aux_vars, bound_cons
end

function bounded_aux(x::Array{JuMP.VariableRef}, y::JuMP.VariableRef, binary_var::JuMP.VariableRef)
    aux_vars, bound_cons = bounded_aux(x, binary_var)
    aux_dep = @variable(y.model)
    push!(bound_cons, @constraint(y.model, aux_dep >= binary_var * JuMP.lower_bound(y)))
    push!(bound_cons, @constraint(y.model, aux_dep <= binary_var * JuMP.upper_bound(y)))
    aux_vars = (aux_vars[1], aux_vars[2], aux_dep)
    return aux_vars, bound_cons
end

"""
    $(TYPEDSIGNATURES)

Generates MI constraints from GlobalModel and its BlackBoxLearners(gm.bbls), and adds them to gm.model.
"""
function add_tree_constraints!(gm::GlobalModel, bbc::BlackBoxClassifier, idx = length(bbc.learners))
    isempty(bbc.mi_constraints) || throw(OCTHaGOnException("BBC $(bbc.name) already has associated MI approximation."))
    isempty(bbc.leaf_variables) || throw(OCTHaGOnException("BBC $(bbc.name) already has associated MI variables."))
    if bbc.feas_ratio == 1.0 # Just a placeholder to show that the tree is "trained".
        z_feas = @variable(gm.model, binary = true)
        bbc.mi_constraints = Dict(1 => [@constraint(gm.model, z_feas == 1)])
        bbc.leaf_variables = Dict(1 => (z_feas, []))
        if !isempty(bbc.lls)
            for ll in bbc.lls
                z_feas = @variable(gm.model, binary = true)
                ll.mi_constraints = Dict(1 => [@constraint(gm.model, z_feas == 1)])
                ll.leaf_variables = Dict(1 => (z_feas, []))
            end
        end
    elseif get_param(bbc, :reloaded)
        bbc.mi_constraints, bbc.leaf_variables = add_feas_constraints!(gm.model, bbc.vars, bbc.learners[idx]; equality = bbc.equality, relax_var = bbc.relax_var, bigM = bbc.bigM)
        if !isempty(bbc.lls)
            for ll in bbc.lls
                ll.mi_constraints, ll.leaf_variables = add_feas_constraints!(gm.model, ll.vars, bbc.learners[idx]; equality = bbc.equality, relax_var = ll.relax_var, bigM = bbc.bigM, symbs = Symbol.(bbc.vars))
            end
        end
    elseif size(bbc.X, 1) == 0
        throw(OCTHaGOnException("Constraint " * string(bbc.name) * " has not been sampled yet, and is thus untrained."))
    elseif isempty(bbc.learners)
        throw(OCTHaGOnException("Constraint " * string(bbc.name) * " must be learned before tree constraints
                            can be generated."))
    elseif bbc.feas_ratio == 0.0
        throw(OCTHaGOnException("Constraint " * string(bbc.name) * " is INFEASIBLE but you tried to include it in
            your global problem. Find at least one feasible solution, train and try again."))
    elseif !get_param(gm, :ignore_accuracy) && !check_accuracy(bbc)
        throw(OCTHaGOnException("Constraint " * string(bbc.name) * " is inaccurately approximated. "))
    else
        bbc.mi_constraints, bbc.leaf_variables = add_feas_constraints!(gm.model, bbc.vars, bbc.learners[idx]; equality = bbc.equality, relax_var = bbc.relax_var, bigM = bbc.bigM)
        if !isempty(bbc.lls)
            for ll in bbc.lls
                ll.mi_constraints, ll.leaf_variables = add_feas_constraints!(gm.model, ll.vars, bbc.learners[idx]; equality = bbc.equality, relax_var = ll.relax_var, bigM = bbc.bigM, symbs = Symbol.(bbc.vars))
            end
        end
    end
    bbc.active_leaves = []
    return
end

function add_tree_constraints!(gm::GlobalModel, bbr::BlackBoxRegressor, idx = length(bbr.learners))
    if size(bbr.X, 1) == 0 && !get_param(bbr, :reloaded)
        throw(OCTHaGOnException("Constraint " * string(bbr.name) * " has not been sampled yet, and is thus untrained."))
    elseif bbr.convex && !bbr.equality
        clear_tree_constraints!(gm, bbr)
        mi_constraints = Dict(1 => [])
        for i = Int64.(ceil.(size(bbr.X,1) .* rand(10)))
            push!(mi_constraints[1], @constraint(gm.model, bbr.dependent_var >= sum(Array(bbr.gradients[i,:]) .* (bbr.vars .- Array(bbr.X[i, :]))) + bbr.Y[i]))
        end
        merge!(append!, bbr.mi_constraints, mi_constraints)
        if !isempty(bbr.lls)
            for lr in bbr.lls
                clear_tree_constraints!(gm, lr)
                mc = Dict(1 => [])
                # Number of initial cuts depends on the dimension of the constraint
                pt_idxs = Int64.(ceil.(size(bbr.X,1) .* rand(maximum([10, Int64(10*ceil(log(length(bbr.vars))))]))))
                update_gradients(bbr, pt_idxs)
                for i in pt_idxs
                    push!(mc[1], @constraint(gm.model, lr.dependent_var >= sum(Array(bbr.gradients[i,:]) .* (lr.vars .- Array(bbr.X[i, :]))) + bbr.Y[i]))
                end
                merge!(append!, lr.mi_constraints, mc)
                lr.active_leaves = [1]
            end
        end
        bbr.active_leaves = [1]
        return
    elseif isempty(bbr.learners)
        throw(OCTHaGOnException("Constraint " * string(bbr.name) * " must be learned before tree constraints
                            can be generated."))
    elseif isempty(bbr.ul_data)
        throw(OCTHaGOnException("Constraint " * string(bbr.name) * " is a Regressor,
        but doesn't have a ORT and/or OCT with upper/lower bounding approximators!"))
    end
    if !isempty(bbr.lls)
        bbr.thresholds[idx].first == "reg" ||
            throw(OCTHaGOnException("BBRs with LinkedRegressors are only allowed ORT approximators."))
        length(bbr.active_trees) == 0 ||
            throw(OCTHaGOnException("BBRs with LinkedRegressors cannot have more than one active tree."))
    elseif bbr.thresholds[idx].first == "reg"
        (isempty(bbr.leaf_variables) || !isnothing(bbr.thresholds[idx].second)) ||
            throw(OCTHaGOnException("Please clear previous tree constraints from $(gm.name) " *
                "before adding an unthresholded regression constraints."))
        if !isnothing(bbr.thresholds[idx].second) && !isempty(bbr.active_trees)
            bbr.thresholds[active_upper_tree(bbr)].second == bbr.thresholds[idx].second ||
                throw(OCTHaGOnException("Upper-thresholded regressors must be preceeded by upper classifiers " *
                                    "of the same threshold for $(bbr.name)."))
        end
    elseif bbr.thresholds[idx].first == "rfreg"
        isempty(bbr.leaf_variables) ||
            throw(OCTHaGOnException("Please clear previous tree constraints from $(gm.name) " *
            "before adding an random forest regressor constraints."))
            !isempty(bbr.lls) && throw(OCTHaGOnException("Random Forests cannot be used to approximate " *
                        "linked BBR $(bbr.name)."))
        isnothing(bbr.thresholds[idx].second) ||
            throw(OCTHaGOnException("RandomForestRegressors are not allowed upper bounds."))
        bbr.equality &&
            throw(OCTHaGOnException("RandomForest cannot approximate BBR $(bbr.name) since it is an equality constraint."))
    elseif bbr.thresholds[idx].first == "upper"
        all(collect(keys(bbr.mi_constraints)) .>= 0)  || throw(OCTHaGOnException("Please clear previous upper tree constraints from $(gm.name) " *
                                                          "before adding new constraints."))
    elseif bbr.thresholds[idx].first == "lower"
        all(collect(keys(bbr.mi_constraints)) .<= 0) || throw(OCTHaGOnException("Please clear previous lower tree constraints from $(gm.name) " *
                                                            "before adding new constraints."))
    end
    if bbr.thresholds[idx].first == "rfreg"
        isempty(bbr.lls) || throw(ErrorException("RandomForestRegressors augment LinkedLearners when MI approximations are made. To avoid bugs, BBRs approximated via RFs cannot have LinkedLearners."))
        trees = get_random_trees(bbr.learners[idx])
        bbr.mi_constraints, bbr.leaf_variables = add_regr_constraints!(gm.model, bbr.vars, bbr.dependent_var,
                    trees[1], bbr.ul_data[idx][1]; equality = bbr.equality, bigM = bbr.bigM, relax_var = bbr.relax_var,)
                    # mic = bbr.mi_constraints, lv = bbr.leaf_variables)
        for i=2:length(trees)
            nll = LinkedRegressor(vars = bbr.vars, dependent_var = bbr.dependent_var)
            nll.mi_constraints, nll.leaf_variables = add_regr_constraints!(gm.model, bbr.vars, bbr.dependent_var,
                                                trees[i], bbr.ul_data[idx][i];
                                                equality = bbr.equality, bigM = bbr.bigM,)
                                                # mic = nll.mi_constraints,
                                                # lv = nll.leaf_variables)
            push!(bbr.lls, nll)
        end
        bbr.active_trees[idx] = bbr.thresholds[idx]
    else
       mi_constraints, leaf_variables = add_regr_constraints!(gm.model, bbr.vars, bbr.dependent_var,
                                                bbr.learners[idx], bbr.ul_data[idx];
                                                equality = bbr.equality, bigM = bbr.bigM, relax_var = bbr.relax_var,)
                                                # mic = bbr.mi_constraints,
                                                # lv = bbr.leaf_variables)

        merge!(append!, bbr.mi_constraints, mi_constraints)
        merge!(bbr.leaf_variables, leaf_variables)
        if !isempty(bbr.lls)
            for lr in bbr.lls
                lr_mi, lr_lv = add_regr_constraints!(gm.model, lr.vars, lr.dependent_var,
                                                bbr.learners[idx], bbr.ul_data[idx];
                                                equality = bbr.equality, bigM = bbr.bigM, relax_var = lr.relax_var,
                                                symbs = Symbol.(bbr.vars),)
                                                # mic = lr.mi_constraints,
                                                # lv = lr.leaf_variables)
                merge!(append!, lr.mi_constraints, lr_mi)
                merge!(lr.leaf_variables, lr_lv)
            end
        end
        if bbr.thresholds[idx].first == "upper"
            push!(bbr.mi_constraints[-1], @constraint(gm.model, bbr.dependent_var <= bbr.thresholds[idx].second))
            [push!(lr.mi_constraints[-1], @constraint(gm.model, lr.dependent_var <= bbr.thresholds[idx])) for lr in bbr.lls]
        elseif bbr.thresholds[idx].first == "lower"
            push!(bbr.mi_constraints[1], @constraint(gm.model, bbr.dependent_var >= bbr.thresholds[idx].second))
            [push!(lr.mi_constraints[1], @constraint(gm.model, lr.dependent_var >= bbr.thresholds[idx])) for lr in bbr.lls]
        end
        bbr.active_trees[idx] = bbr.thresholds[idx]
    end
    bbr.active_leaves = []
    [lr.active_leaves = [] for lr in bbr.lls]
    return
end

function add_tree_constraints!(gm::GlobalModel, bbls::Vector{BlackBoxLearner})
    for bbl in bbls
        add_tree_constraints!(gm, bbl)
    end
    return
end

add_tree_constraints!(gm::GlobalModel) = add_tree_constraints!(gm, gm.bbls)

""" Computes the smallest possible big-M for a α'x ≥ threshold split in a tree. """
function compute_hyperplane_bigM_lower(threshold::Real, α::Vector, var_bounds::Vector)
    return threshold - sum(minimum(α[i] .* var_bounds[i]) for i=1:length(α))
end

""" Computes the smallest possible big-M for a α'x ≤ threshold split in a tree. """
function compute_hyperplane_bigM_upper(threshold::Real, α::Vector, var_bounds::Vector)
    return -threshold + sum(maximum(α[i] .* var_bounds[i]) for i=1:length(α))
end

""" Computes the smallest possible big-M for lower bounding regression in a tree leaf. """
function compute_regression_bigM_lower(β0::Real, β::Vector, var_bounds::Vector, y_bounds::Vector)
    return β0 + sum(maximum(β[i] .* var_bounds[i]) for i=1:length(β)) - minimum(y_bounds)
end

""" Computes the smallest possible big-M for upper bounding regression in a tree leaf. """
function compute_regression_bigM_upper(β0::Real, β::Vector, var_bounds::Vector, y_bounds::Vector)
    return -β0 - sum(minimum(β[i] .* var_bounds[i]) for i=1:length(β)) + maximum(y_bounds)
end

"""
Creates a set of binary feasibility constraints from a binary classification tree.
Arguments:
- m: JuMP Model
- x: independent JuMP.Variables (features in lnr)
- lnr: A fitted OptimalTreeClassifier
- equality: whether the constraint is an equality.
NOTE: mic and lv are only nonempty if we are adding an OCT approximation on top of an existing set of MI constraints and variables respectively. Leave defaults empty for basic usage.
symbs helps add constraints to linked learners, since the varkeys in lnr are different to the variables in the linked learner.
"""
function add_feas_constraints!(m::JuMP.Model, x::Array{JuMP.VariableRef}, lnr::IAI.OptimalTreeLearner;
                               equality::Bool = false, bigM::Bool = false,
                               relax_var::Union{Real, JuMP.VariableRef} = 0,
                               mic::Dict = Dict(), lv::Dict = Dict(),
                               symbs = Symbol.(x))
    check_if_trained(lnr);
    all_leaves = find_leaves(lnr)
    # Add a binary variable for each leaf
    feas_leaves = []
    if lnr isa IAI.OptimalTreeClassifier
        feas_leaves =
        [i for i in all_leaves if Bool(IAI.get_classification_label(lnr, i))];
    else
        feas_leaves = all_leaves
    end
    mi_constraints = Dict(leaf => [] for leaf in feas_leaves)
    mi_constraints[1] = []
    leaf_variables = Dict{Int64, Tuple{JuMP.VariableRef, Array}}()
    var_bounds = get_bounds(x)
    var_bounds = [var_bounds[var] for var in x] # ordered vector form
    if isempty(lv) # if the variables haven't already been created and passed!
        if bigM
            # big-M formulation variable generation and coupling
            for leaf in feas_leaves
                leaf_variables[leaf] = (@variable(m, binary = true), [])
            end
            push!(mi_constraints[1], @constraint(m, sum(leaf_variables[leaf][1] for leaf in feas_leaves) == 1))
        else
            # big-M free formulation variable generation and coupling
            for leaf in feas_leaves
                leaf_variables[leaf], bd_mi = bounded_aux(x, @variable(m, binary=true))
                append!(mi_constraints[leaf], bd_mi)
            end
            push!(mi_constraints[1], @constraint(m, sum(leaf_variables[leaf][1] for leaf in feas_leaves) == 1))
            [push!(mi_constraints[1], constr)
                for constr in @constraint(m, sum(leaf_variables[leaf][2] for leaf in feas_leaves) .== x)]
        end
    else
        leaf_variables = lv
    end
    # Getting lnr data
    upperDict, lowerDict = trust_region_data(lnr, symbs)
    for leaf in feas_leaves
        # ADDING TRUST REGIONS
        for (threshold, α) in upperDict[leaf] # upperDict is lower-bounded.
            if bigM
                M = compute_hyperplane_bigM_lower(threshold, α, var_bounds)
                push!(mi_constraints[leaf], @constraint(m, threshold ≤ sum(α .* x) + M * (1 - leaf_variables[leaf][1] + relax_var)))
            else
                push!(mi_constraints[leaf], @constraint(m, threshold * leaf_variables[leaf][1] <= sum(α .* leaf_variables[leaf][2]) + relax_var))
            end
        end
        for (threshold, α) in lowerDict[leaf] # lowerDict is upper-bounded.
            if bigM
                M = compute_hyperplane_bigM_upper(threshold, α, var_bounds)
                push!(mi_constraints[leaf], @constraint(m, M * (1 - leaf_variables[leaf][1] + relax_var) + threshold ≥ sum(α .* x)))
            else
                push!(mi_constraints[leaf], @constraint(m, threshold * leaf_variables[leaf][1] + relax_var >= sum(α .* leaf_variables[leaf][2])))
            end
        end
    end
    if equality
        infeas_leaves = [i for i in all_leaves if !Bool(IAI.get_classification_label(lnr, i))];
        [mi_constraints[leaf] = [] for leaf in infeas_leaves]
        if isempty(lv)
            if bigM
                # big-M formulation variable generation and coupling
                for leaf in infeas_leaves
                    leaf_variables[leaf] = (@variable(m, binary = true), [])
                end
                push!(mi_constraints[1], @constraint(m, sum(leaf_variables[leaf][1] for leaf in infeas_leaves) == 1))
            else
                for leaf in infeas_leaves
                    leaf_variables[leaf], bd_mi = bounded_aux(x, @variable(m, binary=true))
                    append!(mi_constraints[leaf], bd_mi)
                end
                push!(mi_constraints[1], @constraint(m, sum(leaf_variables[leaf][1] for leaf in infeas_leaves) == 1))
                [push!(mi_constraints[1], constr)
                    for constr in @constraint(m, sum(leaf_variables[leaf][2] for leaf in infeas_leaves) .== x)]
            end
        end
        for leaf in infeas_leaves
            # ADDING TRUST REGIONS
            for (threshold, α) in upperDict[leaf] # upperDict is lower-bounded.
                if bigM
                    M = compute_hyperplane_bigM_lower(threshold, α, var_bounds)
                    push!(mi_constraints[leaf], @constraint(m, threshold ≤ sum(α .* x) + M * (1 - leaf_variables[leaf][1] + relax_var)))
                else
                    push!(mi_constraints[leaf], @constraint(m, threshold * leaf_variables[leaf][1] <= sum(α .* leaf_variables[leaf][2]) + relax_var))
                end
            end
            for (threshold, α) in lowerDict[leaf] # lowerDict is upper-bounded.
                if bigM
                    M = compute_hyperplane_bigM_upper(threshold, α, var_bounds)
                    push!(mi_constraints[leaf], @constraint(m, threshold + M * (1 - leaf_variables[leaf][1] + relax_var) >= sum(α .* x)))
                else
                    push!(mi_constraints[leaf], @constraint(m, threshold * leaf_variables[leaf][1] + relax_var >= sum(α .* leaf_variables[leaf][2])))
                end
            end
        end
    end
    merge!(append!, mi_constraints, mic)
    return mi_constraints, leaf_variables
end

"""
    $(TYPEDSIGNATURES)

Creates a set of MIO constraints from a OptimalTreeClassifier that thresholds a BlackBoxRegressor.

Arguments:
- m: JuMP.Model
- x:: independent JuMP.Variables (features of learner)
- y:: dependent JuMP.Variables
- lnr:: A fitted OptimalTreeLearner
- ul_data:: Upper and lower bounding hyperplanes for data in leaves of lnr (empty by default)
- symbs:: The varkeys corresponding to the lnr.
"""
function add_regr_constraints!(m::JuMP.Model, x::Array{JuMP.VariableRef}, y::JuMP.VariableRef,
                               lnr::IAI.OptimalTreeLearner,
                               ul_data::Dict; equality::Bool = false, bigM::Bool = false,
                               relax_var::Union{Real, JuMP.VariableRef} = 0, symbs = Symbol.(x))
    # TODO: Determine whether I should relax approximator or trust regions or both.
    if lnr isa OptimalTreeRegressor
        check_if_trained(lnr)
        all_leaves = find_leaves(lnr)
        # Add a binary variable for each leaf
        mi_constraints = Dict(leaf => [] for leaf in all_leaves)
        mi_constraints[1] = []
        leaf_variables = Dict{Int64, Tuple{JuMP.VariableRef, Array, JuMP.VariableRef}}()
        var_bounds = get_bounds(x)
        var_bounds = [var_bounds[var] for var in x]
        y_bounds = get_bound(y)
        if bigM
            for leaf in all_leaves
                leaf_variables[leaf] = (@variable(m, binary = true), [], @variable(m)) # third item in tuple is placeholder. TODO: find a more elegant bookkeeping method.
            end
            push!(mi_constraints[1], @constraint(m, sum(leaf_variables[leaf][1] for leaf in all_leaves) == 1))
        else
            for leaf in all_leaves
                leaf_variables[leaf], bd_mi = bounded_aux(x, y, @variable(m, binary=true))
                append!(mi_constraints[leaf], bd_mi)
            end
            push!(mi_constraints[1], @constraint(m, sum(leaf_variables[leaf][1] for leaf in all_leaves) == 1))
            [push!(mi_constraints[1], constr)
                for constr in @constraint(m, sum(leaf_variables[leaf][2] for leaf in all_leaves) .== x)]
            push!(mi_constraints[1], @constraint(m, sum(leaf_variables[leaf][3] for leaf in all_leaves) == y))
        end
        # Getting lnr data
        pwlDict = pwl_constraint_data(lnr, symbs)
        upperDict, lowerDict = trust_region_data(lnr, symbs)
        for leaf in all_leaves
            # ADDING PWL APPROXIMATIONS
            β0, β = pwlDict[leaf]
            if bigM
                if equality
                    M_upper = compute_regression_bigM_upper(β0, β, var_bounds, y_bounds)
                    M_lower = compute_regression_bigM_lower(β0, β, var_bounds, y_bounds)
                    # Use the ridge regressor!
                    push!(mi_constraints[leaf], @constraint(m, sum(β .* x) + β0 + M_upper * (1 - leaf_variables[leaf][1] + relax_var) >= y))
                    push!(mi_constraints[leaf], @constraint(m, sum(β .* x) + β0 <= y + M_lower * (1 - leaf_variables[leaf][1] + relax_var)))
                elseif !isempty(ul_data)
                    # Use only the lower approximator
                    β0, β = ul_data[leaf] # update the lower approximator
                    M_lower = compute_regression_bigM_lower(β0, β, var_bounds, y_bounds)
                    push!(mi_constraints[leaf], @constraint(m, sum(β .* x) + β0 <= y + M_lower * (1 - leaf_variables[leaf][1] + relax_var)))
                else
                    # Use only the ridge regressor as the lower approximator.
                    M_lower = compute_regression_bigM_lower(β0, β, var_bounds, y_bounds)
                    push!(mi_constraints[leaf], @constraint(m, sum(β .* x) + β0 <= y + M_lower * (1 - leaf_variables[leaf][1] + relax_var)))
                end
            else
                if equality
                    # Use the ridge regressor!
                    push!(mi_constraints[leaf], @constraint(m, sum(β .* leaf_variables[leaf][2]) + β0 * leaf_variables[leaf][1] + relax_var >= leaf_variables[leaf][3]))
                    push!(mi_constraints[leaf], @constraint(m, sum(β .* leaf_variables[leaf][2]) + β0 * leaf_variables[leaf][1] <= leaf_variables[leaf][3] + relax_var))
                elseif !isempty(ul_data)
                    # Use only the lower approximator
                    β0, β = ul_data[leaf] # update the lower approximator
                    push!(mi_constraints[leaf], @constraint(m, sum(β .* leaf_variables[leaf][2]) + β0 * leaf_variables[leaf][1] <= leaf_variables[leaf][3] + relax_var))
                else
                    # Use only the ridge regressor as the lower approximator.
                    push!(mi_constraints[leaf], @constraint(m, sum(β .* leaf_variables[leaf][2]) + β0 * leaf_variables[leaf][1] <= leaf_variables[leaf][3] + relax_var))
                end
            end
            # ADDING TRUST REGIONS
            for (threshold, α) in upperDict[leaf] # upperDict is lower-bounded.
                if bigM
                    M = compute_hyperplane_bigM_lower(threshold, α, var_bounds)
                    push!(mi_constraints[leaf], @constraint(m, threshold ≤ sum(α .* x) + M * (1 - leaf_variables[leaf][1] + relax_var)))
                else
                    push!(mi_constraints[leaf], @constraint(m, threshold * leaf_variables[leaf][1] <= sum(α .* leaf_variables[leaf][2]) + relax_var))
                end
            end
            for (threshold, α) in lowerDict[leaf] # lowerDict is upper-bounded.
                if bigM
                    M = compute_hyperplane_bigM_upper(threshold, α, var_bounds)
                    push!(mi_constraints[leaf], @constraint(m, M * (1 - leaf_variables[leaf][1] + relax_var) + threshold ≥ sum(α .* x)))
                else
                    push!(mi_constraints[leaf], @constraint(m, threshold * leaf_variables[leaf][1] + relax_var >=
                    sum(α .* leaf_variables[leaf][2])))
                end
            end
        end
        return mi_constraints, leaf_variables
    elseif lnr isa OptimalTreeClassifier
        # Add a binary variable for each leaf
        all_leaves = find_leaves(lnr)
        # Add a binary variable for each leaf
        feas_leaves = [i for i in all_leaves if Bool(IAI.get_classification_label(lnr, i))]
        mi_constraints = Dict(leaf => [] for leaf in feas_leaves)
        mi_constraints[1] = []
        leaf_variables = Dict{Int64, Tuple{JuMP.VariableRef, Array, JuMP.VariableRef}}()
        for leaf in feas_leaves
            leaf_variables[leaf], bd_mi = bounded_aux(x, y, @variable(m, binary=true))
            append!(mi_constraints[leaf], bd_mi)
        end
        push!(mi_constraints[1], @constraint(m, sum(leaf_variables[leaf][1] for leaf in feas_leaves) == 1))
        [push!(mi_constraints[1], constr)
            for constr in @constraint(m, sum(leaf_variables[leaf][2] for leaf in feas_leaves) .== x)]
        push!(mi_constraints[1], @constraint(m, sum(leaf_variables[leaf][3] for leaf in feas_leaves) == y))
        mi_constraints, leaf_variables = add_feas_constraints!(m, x, lnr, equality = equality, bigM = bigM, mic = mi_constraints, lv = leaf_variables)
        if !isempty(ul_data)
            if all(keys(ul_data) .<= 0) || !all(keys(ul_data) .>= 0) # means an upper or upperlower bounding tree
                mi_constraints = Dict(-key => value for (key, value) in mi_constraints) # hacky sign flipping for upkeep.
                leaf_variables = Dict(-key => value for (key, value) in leaf_variables)
            end
            for leaf in collect(keys(ul_data))
                γ0, γ = ul_data[leaf]
                if !haskey(mi_constraints, leaf) # occurs with ORTS with bounding hyperplanes
                    mi_constraints[leaf] = [@constraint(m, leaf_variables[-leaf][3] <= γ0 * leaf_variables[-leaf][1] +
                        sum(γ .* leaf_variables[-leaf][2]) + relax_var)]
                elseif leaf <= 0
                    push!(mi_constraints[leaf], @constraint(m, leaf_variables[leaf][3] <= γ0 * leaf_variables[leaf][1] +
                        sum(γ .* leaf_variables[leaf][2]) + relax_var))
                else
                    push!(mi_constraints[leaf], @constraint(m, leaf_variables[leaf][3] + relax_var >=
                    γ0 * leaf_variables[leaf][1] + sum(γ .* leaf_variables[leaf][2])))
                end
            end
        end
        return mi_constraints, leaf_variables
    end
end

""" Clears upper-bounding constraints from a BBR and its associated GM. """
function clear_upper_constraints!(gm, bbr::Union{BlackBoxRegressor, LinkedRegressor})
    for (leaf_key, leaf_constrs) in bbr.mi_constraints
        if leaf_key <= 0
            while !isempty(leaf_constrs)
                constr = pop!(leaf_constrs)
                if is_valid(gm.model, constr)
                    delete(gm.model, constr)
                else
                    push!(leaf_constrs, constr) # make sure to put the constraint back.
                    throw(OCTHaGOnException("Bug: Constraints could not be removed."))
                end
            end
            delete!(bbr.mi_constraints, leaf_key)
        end
    end
    for (leaf_key, (bin_var, leaf_vars, aux_dep)) in bbr.leaf_variables
        if leaf_key <= 0
            for leaf_var in leaf_vars
                if is_valid(gm.model, leaf_var)
                    delete(gm.model, leaf_var)
                else
                    throw(OCTHaGOnException("Bug: Variables could not be removed."))
                end
            end
            if is_valid(gm.model, bin_var)
                delete(gm.model, bin_var)
            else
                throw(OCTHaGOnException("Bug: Binary variable could not be removed."))
            end
            if is_valid(gm.model, aux_dep)
                delete(gm.model, aux_dep)
            else
                throw(OCTHaGOnException("Bug: Aux dependent variable could not be removed."))
            end
            delete!(bbr.leaf_variables, leaf_key)
        end
    end
    if bbr isa BlackBoxRegressor
        idx = active_upper_tree(bbr)
        delete!(bbr.active_trees, idx)
        for lr in bbr.lls
            clear_upper_constraints!(gm, lr)
        end
    end
    return
end

""" Clears lower-bounding constraints from a BBR and its associated GM. """
function clear_lower_constraints!(gm, bbr::Union{BlackBoxRegressor, LinkedRegressor})
    for (leaf_key, leaf_constrs) in bbr.mi_constraints
        if leaf_key >= 0
            while !isempty(leaf_constrs)
                constr = pop!(leaf_constrs)
                if is_valid(gm.model, constr)
                    delete(gm.model, constr)
                else
                    push!(leaf_constrs, constr) # make sure to put the constraint back.
                    throw(OCTHaGOnException("Bug: Constraints could not be removed."))
                end
            end
            delete!(bbr.mi_constraints, leaf_key)
        end
    end
    for (leaf_key, (bin_var, leaf_vars, aux_dep)) in bbr.leaf_variables
        if leaf_key >= 0
            for leaf_var in leaf_vars
                if is_valid(gm.model, leaf_var)
                    delete(gm.model, leaf_var)
                else
                    throw(OCTHaGOnException("Bug: Variables could not be removed."))
                end
            end
            if is_valid(gm.model, bin_var)
                delete(gm.model, bin_var)
            else
                throw(OCTHaGOnException("Bug: Binary variable could not be removed."))
            end
            if is_valid(gm.model, aux_dep)
                delete(gm.model, aux_dep)
            else
                throw(OCTHaGOnException("Bug: Aux dependent variable could not be removed."))
            end
            delete!(bbr.leaf_variables, leaf_key)
        end
    end
    bbr.active_leaves = []
    if bbr isa BlackBoxRegressor
        idx = active_lower_tree(bbr)
        delete!(bbr.active_trees, idx)
        for lr in bbr.lls
            clear_lower_constraints!(gm, lr)
        end
    end
    return
end

"""
    $(TYPEDSIGNATURES)

Clears the MI-approximating constraints and variables in GlobalModel and its sub-structs.
"""
function clear_tree_constraints!(gm::GlobalModel, bbc::Union{BlackBoxClassifier, LinkedClassifier})
    for (leaf_key, leaf_constrs) in bbc.mi_constraints
        while !isempty(leaf_constrs)
            constr = pop!(leaf_constrs)
            if is_valid(gm.model, constr)
                delete(gm.model, constr)
            else
                push!(leaf_constrs, constr) # make sure to put the constraint back.
                throw(OCTHaGOnException("Bug: Constraints could not be removed."))
            end
        end
        delete!(bbc.mi_constraints, leaf_key)
    end
    for (leaf_key, (bin_var, leaf_vars)) in bbc.leaf_variables
        for leaf_var in leaf_vars
            if is_valid(gm.model, leaf_var)
                delete(gm.model, leaf_var)
            else
                throw(OCTHaGOnException("Bug: Variables could not be removed."))
            end
        end
        if is_valid(gm.model, bin_var)
            delete(gm.model, bin_var)
        else
            throw(OCTHaGOnException("Bug: Binary variable could not be removed."))
        end
        delete!(bbc.leaf_variables, leaf_key)
    end
    bbc.active_leaves = []
    if bbc isa BlackBoxClassifier
        for lc in bbc.lls
            clear_tree_constraints!(gm, lc)
        end
    end
    return
end

function clear_tree_constraints!(gm::GlobalModel, bbr::Union{BlackBoxRegressor, LinkedRegressor})
    clear_lower_constraints!(gm, bbr)
    clear_upper_constraints!(gm, bbr)
    return
end

function clear_tree_constraints!(gm::GlobalModel, bbls::Array{BlackBoxLearner})
    for bbl in bbls
        clear_tree_constraints!(gm, bbl)
    end
    return
end

clear_tree_constraints!(gm::GlobalModel) = clear_tree_constraints!(gm, gm.bbls)

"""
    update_tree_constraints!(gm::GlobalModel, bbr::BlackBoxRegressor, idx = length(bbr.learners))
    update_tree_constraints!(gm::GlobalModel, bbc::BlackBoxClassifier, idx = length(bbc.learners))

Updates the MI constraints associated with a BBL.
For BBRs, makes sure to replace the appropriate lower/upper/regressor approximations.
"""
function update_tree_constraints!(gm::GlobalModel, bbr::BlackBoxRegressor, idx = length(bbr.learners))
    if isempty(bbr.active_trees) # no mi constraints yet.
        add_tree_constraints!(gm, bbr, idx)
    elseif length(bbr.active_trees) == 1 # one set of mi_constraints
        new_threshold = bbr.thresholds[idx]
        active_idx = collect(keys(bbr.active_trees))[1]
        last_threshold = bbr.active_trees[active_idx]
        if new_threshold.first == last_threshold.first
            clear_tree_constraints!(gm, bbr)
            add_tree_constraints!(gm, bbr, idx)
            return
        elseif new_threshold == Pair("reg", nothing) || last_threshold == Pair("reg", nothing) ||
            last_threshold.first == "upperlower" || new_threshold.first == "upperlower"
            clear_tree_constraints!(gm, bbr)
            add_tree_constraints!(gm, bbr, idx)
            return
        elseif new_threshold.first == "reg" && last_threshold.first == "upper" ||
                last_threshold.first == "reg" && new_threshold.first == "upper"
            if new_threshold.second == last_threshold.second # if corresponding upper thresholds.
                add_tree_constraints!(gm, bbr, idx)
                return
            else
                clear_tree_constraints!(gm, bbr)
                add_tree_constraints!(gm, bbr, idx)
                return
            end
        elseif last_threshold.first in valid_uppers && new_threshold.first in valid_lowers ||
            last_threshold.first in valid_lowers && new_threshold.first in valid_uppers
            add_tree_constraints!(gm, bbr, idx)
            return
        else # updating a lower or upper threshold
            clear_tree_constraints!(gm, bbr)
            add_tree_constraints!(gm, bbr, idx)
            return
        end
    elseif length(bbr.active_trees) == 2 # two sets of mi_constraints
        if bbr.thresholds[idx] == Pair("reg", nothing) || bbr.thresholds[idx].first == "upperlower"
            clear_tree_constraints!(gm, bbr)
            add_tree_constraints!(gm, bbr, idx)
            return
        end
        hypertype = bbr.thresholds[idx].first # otherwise, check approximation type
        if hypertype in valid_lowers
            clear_lower_constraints!(gm, bbr)
        elseif hypertype in valid_uppers
            clear_upper_constraints!(gm, bbr)
        else
            throw(OCTHaGOnException("Hypertype $(hypertype) not recognized."))
        end
        add_tree_constraints!(gm, bbr, idx)
        return
    else
        throw(OCTHaGOnException("$(gm.name) has too many active lower/upper bounding regressors."))
    end
end

function update_tree_constraints!(gm::GlobalModel, bbc::BlackBoxClassifier, idx = length(bbc.learners))
    clear_tree_constraints!(gm, bbc)
    add_tree_constraints!(gm, bbc, idx)
    return
end

ext: jl


text:
module Octavian

using Requires: @require

using VectorizationBase, ArrayInterface, LoopVectorization

using VectorizationBase: align, AbstractStridedPointer, zstridedpointer, vsub_nsw, assume,
  static_sizeof, StridedPointer, gesp, pause, pick_vector_width, has_feature
using CPUSummary: cache_size, num_cores, cache_inclusive, cache_linesize
using LoopVectorization: preserve_buffer, CloseOpen, UpperBoundedInteger
using ArrayInterface: size, strides, offsets, indices, axes, StrideIndex
using IfElse: ifelse
using PolyesterWeave
using Static: StaticInt, Zero, One, StaticBool, True, False, gt, eq, StaticFloat64,
    roundtostaticint, floortostaticint
using ManualMemory: MemoryBuffer, load, store!

using ThreadingUtilities: _atomic_add!, _atomic_load, _atomic_store!, launch, wait, SPIN

export StaticInt
export matmul!
export matmul
export matmul_serial!
export matmul_serial

debug() = false

include("global_constants.jl")
include("types.jl")
include("integerdivision.jl")
include("memory_buffer.jl")
include("block_sizes.jl")
include("funcptrs.jl")
include("macrokernels.jl")
include("utils.jl")
include("matmul.jl")
include("complex_matmul.jl")

include("init.jl") # `Octavian.__init__()` is defined in this file

@static if VERSION >= v"1.8.0-beta1"
  let
    __init__()
    A64 = rand(100,100)
    matmul(A64,A64)
    matmul(A64',A64)
    matmul(A64,A64')
    matmul(A64',A64')
    A32 = rand(Float32,100,100)
    matmul(A32,A32)
    matmul(A32',A32)
    matmul(A32,A32')
    matmul(A32',A32')
  end
end

end # module Octavian


text:
# Causal Intervention (using TaggedΩ) #
"Intervened Random Variable x | do(theta)"
struct ReplaceRandVar{R1 <: RandVar, R2 <: RandVar} <: RandVar
  x::R1                     # Intervened Random Variable
  replmap::Dict{ID, R2}     # Intervention map
  id::ID
  ReplaceRandVar(rv::R22, replmap::Dict{ID, T2}, id = uid()) where {R22 <: RandVar, T2 <: RandVar} =
    new{R22, T2}(rv, replmap, id)
end
@inline (rv::ReplaceRandVar)(ω::Ω) = apl(rv, ω)
id(x::ReplaceRandVar) = x.id
ppapl(rv::ReplaceRandVar, ω::Ω) = apl(rv.x, tag(ω, (replmap = rv.replmap,)))

maybereplace(x::RandVar, replmap) = replace(x, replmap)

# If x is a constant, do nothing
maybereplace(x, replmap) = x

# Params of intervened RandVar should themselves be intervened (if not const)
params(rv::ReplaceRandVar) = map(p -> maybereplace(p, rv.replmap), params(rv.x))

@inline function replaceapl(rv::RandVar, tω::TaggedΩ{I, Tags{K, V}, ΩT}) where {I, K, V, ΩT <: ΩBase}
  if haskey(tω.tags.replmap, rv.id) # FIXME, double look up
    # T = (Core.Compiler).return_type(rv, typeof((tω.taggedω,)))
    # @show T
    # @show tω.taggedω
    apl(tω.tags.replmap[rv.id], tω)#::T
  else
    ppapl(rv, proj(tω, rv))
  end
end

# Conversions #
mcv(x::RandVar) = x
mcv(x) = constant(x)
upconv(x::Dict{RV}) where RV = Dict(k.id => mcv(v) for (k, v) in x)
upconv(pairs::Pair...) = Dict(k.id => mcv(v) for (k, v) in pairs)
upconv(pair::Pair) = Dict(pair.first.id => mcv(pair.second))

# """`replace(x::RandVar, replmap)`
# Causal intervention.

# Returns a randvar that is like `x` but where random variables it depend on
# (its parents) are changed to those specified in replmap

# ```julia
# θ = normal(0, 1)
# x = normal(θ, 1)
# xdo1 = replace(x, θ => 100.0)
# xdo2 = replace(x, θ => uniform(1000.0, 2000.0))
# rand((x, xdo1, xdo2))
# ```

# """
# function Base.replace end

"Causal Intervention: Set `θold` to `θnew` in `x`"
Base.replace(x::RandVar, replmap::Dict{Int, <: RandVar}) = ReplaceRandVar(x, replmap)
@spec :nocheck all([isparent(theta, x) for theta in values(tochange)])

"Causal Intervention: Set `θold` to `θnew` in `x`"
Base.replace(x::RandVar, tochange::Union{Dict, Pair}...) = replace(x, upconv(tochange...))

## Show
Base.show(io::IO, rv::ReplaceRandVar) =
  (print(io, rv.x); print(io, " | intervened "))


text:
struct Quadratic
    P
end
Flux.@functor Quadratic

(q::Quadratic)(x) = sum((x .+ q.P) .* (x .+ q.P))

@testset "Testing Simple Training" begin
    # We'd expect P to be forced to zero to minimize this quadratic objective.
    P = Float32[1.0, 1.0, 1.0]
    f = Quadratic(P)

    x = Float32[1.0, 2.0, 3.0]
    backend = nGraph.Backend()

    learning_rate = Float32(0.2)
    ex = nGraph.compile(backend, f, x; optimizer = nGraph.SGD(learning_rate))

    # We should have one implicit input/output: P
    @test length(ex.optimizer.inputs) == 1
    @test length(ex.optimizer.outputs) == 1

    @show parent(first(ex.optimizer.inputs))
    @show parent(first(ex.optimizer.outputs))

    Z = ex()
    @show parent(Z)

    @show parent(first(ex.optimizer.inputs))
    @show parent(first(ex.optimizer.outputs))

    @test all(parent(first(ex.optimizer.inputs)) .< parent(first(ex.optimizer.outputs)))

    # We can compute what the derivatives should be. Note that the inputs/outputs should have
    # swapped, so now we're looking at the implicit inputs
    dfdx(i) = 2 * (f.P[i] + x[i])

    @test P[1] - learning_rate * dfdx(1) == parent(ex.optimizer.inputs[1])[1]
    @test P[2] - learning_rate * dfdx(2) == parent(ex.optimizer.inputs[1])[2]
    @test P[3] - learning_rate * dfdx(3) == parent(ex.optimizer.inputs[1])[3]

    # Run this for a few more iterations - see if we can get the implicit parameters to be
    # negative the input
    for i in 1:20
        println(parent(ex()))
    end

    @test parent(ex())[] < 1E-3
end


text:
# using Base.Test
# using MIToS.MSA

let alphabet = "ARNDCQEGHILKMFPSTWYV-",
  msa =  ["DAWAEF",
          "DAWAED",
          "DAYCMD" ],
  badmsa=["DAWAEF",
          "DAWAED",
          "DAM"    ]

  print("""

Test Residue
============
""")

  print("""
Convert
""")
  @test Int[ Residue(char) for char in alphabet] == Int[ i for i in 1:21 ]
  @test Int[ Residue(char) for char in lowercase("ARNDCQEGHILKMFPSTWYV.")] == Int[ 21 for i in 1:21 ]

  print("""
GAP
""")
  @test Int(GAP) == 21
  @test Int8(GAP) == Int8(21)
  @test Residue('-') == GAP
  @test Char(GAP) == '-'
  @test UInt8(GAP) == UInt8('-')

  print("""
Convert to and from ASCIIString
""")
  @test res"ARNDCQEGHILKMFPSTWYV-" == Residue[ Residue(char) for char in alphabet]
  @test Residue(alphabet) == Residue[ Residue(char) for char in alphabet]
  @test ascii(Residue(alphabet)) == alphabet
  @test string(Residue(alphabet)) == alphabet # "AR..." instead of the standar "[A,R,..."

  @test convert(Matrix{Residue}, msa) == Residue['D' 'A' 'W' 'A' 'E' 'F'
                                                 'D' 'A' 'W' 'A' 'E' 'D'
                                                 'D' 'A' 'Y' 'C' 'M' 'D']

  @test_throws ErrorException convert(Matrix{Residue}, badmsa)

  print("""
Ambiguous or not standard residues (are gaps on MIToS)
""")

  for res in Char['X', # Any amino acid
                  'B', # Aspartic acid or Asparagine
                  'Z', # Glutamine or Glutamic acid
                  'O', # Pyrrolysine
                  'U', # Selenocysteine
                  'J'] # Leucine or Isoleucine

    @test Residue(res) == GAP

  end

  print("""
Comparisons
""")

  for res in Residue(alphabet)
    @test res == res
  end

  for i in 1:21, j in 1:21
    if j != i
      @test Residue(i) != Residue(j)
    end
  end

  print("""
Random
""")
  for i in 1:8000
    @test rand(Residue) != GAP
    @test typeof( rand(Residue) ) == Residue
  end
  @test size(rand(Residue, 20)) == (20,)
  @test typeof(rand(Residue, 20)) == Array{Residue, 1}
  @test size(rand(Residue, 20,30)) == (20,30)
  @test typeof(rand(Residue, 20, 30)) == Array{Residue, 2}

end

print("""

Three letters <=> One letter (Residue)
--------------------------------------
""")

@test three2residue("ALA") == Residue('A')
@test three2residue("Ala") == Residue('A')
@test three2residue("ala") == Residue('A')

@test_throws ErrorException three2residue("AL")
@test_throws ErrorException three2residue("Alanine")

@test residue2three(Residue('A')) == "ALA"


text:
# This file was generated, do not modify it. # hide
using StatsBase
describe(boston, :min, :max, :mean, :median, :std)


text:
module _TestIndirectImportsDownstream2

using IndirectImports

@indirect import _TestIndirectImportsUpstream

end # module


text:
using ModelingToolkitStandardLibrary.Electrical, ModelingToolkit, OrdinaryDiffEq, Test
using ModelingToolkitStandardLibrary.Electrical: _step, _square_wave, _triangular_wave, _cos_wave, _damped_sine_wave, _ramp

# using Plots

@parameters t

@testset "sensors" begin
    @named source = ConstantVoltage(V=10)
    @named resistor = Resistor(R=1)
    @named capacitor = Capacitor(C=1)
    @named ground = Ground()

    @named voltage_sensor = VoltageSensor()
    @named current_sensor = CurrentSensor()
    @named power_sensor = PowerSensor()

    connections = [
        connect(source.p, resistor.p)
        connect(resistor.n, current_sensor.p)
        connect(current_sensor.n, power_sensor.pc)
        connect(power_sensor.nc, capacitor.p)
        connect(capacitor.n, source.n, ground.g)
        connect(capacitor.p, voltage_sensor.p)
        connect(capacitor.n, voltage_sensor.n)
        connect(capacitor.p, power_sensor.pv)
        connect(capacitor.n, power_sensor.nv)
    ]

    @named model = ODESystem(connections, t; systems=[resistor, capacitor, source, ground, voltage_sensor, current_sensor, power_sensor])
    sys = structural_simplify(model)
    prob = ODAEProblem(sys, Pair[], (0.0, 10.0))
    sol = solve(prob, Tsit5())

    # Plots.plot(sol; vars=[capacitor.v, voltage_sensor.v])
    # Plots.plot(sol; vars=[power_sensor.power, capacitor.i * capacitor.v])
    # Plots.plot(sol; vars=[resistor.i, current_sensor.i])
    @test sol[capacitor.v] ≈ sol[voltage_sensor.v] atol=1e-3
    @test sol[power_sensor.power] ≈ sol[capacitor.i * capacitor.v] atol=1e-3
    @test sol[resistor.i] ≈ sol[current_sensor.i] atol=1e-3
end

# simple voltage divider
@testset "voltage divider" begin
    @named source = ConstantVoltage(V=10)
    @named R1 = Resistor(R=1e3)
    @named R2 = Resistor(R=1e3)
    @named ground = Ground()

    connections = [
        connect(source.p, R1.p)
        connect(R1.n, R2.p)
        connect(R2.n, source.n, ground.g)
    ]

    @named model = ODESystem(connections, t, systems=[R1, R2, source, ground])
    sys = structural_simplify(model)
    prob = ODEProblem(sys, Pair[], (0, 2.0))
    sol = solve(prob, Rodas4()) # has no state; does not work with Tsit5
    @test sol[R1.p.v][end] ≈ 10 atol=1e-3
    @test sol[R1.n.v][end] ≈ 5 atol=1e-3
    @test sol[R2.n.v][end] ≈ 0 atol=1e-3
end

# simple RC
@testset "RC" begin
    @named source = ConstantVoltage(V=10)
    @named resistor = Resistor(R=1)
    @named capacitor = Capacitor(C=1)
    @named ground = Ground()

    connections = [
        connect(source.p, resistor.p)
        connect(resistor.n, capacitor.p)
        connect(capacitor.n, source.n, ground.g)
    ]

    @named model = ODESystem(connections, t; systems=[resistor, capacitor, source, ground])
    sys = structural_simplify(model)
    prob = ODAEProblem(sys, [capacitor.v => 0.0], (0.0, 10.0))
    sol = solve(prob, Tsit5())

    # Plots.plot(sol; vars=[source.v, capacitor.v])
    @test sol[capacitor.v][end] ≈ 10 atol=1e-3
end

# simple RL
@testset "RL" begin
    @named source = ConstantVoltage(V=10)
    @named resistor = Resistor(R=1)
    @named inductor = Inductor(L=1.0)
    @named ground = Ground()

    connections = [
        connect(source.p, resistor.p)
        connect(resistor.n, inductor.p)
        connect(inductor.n, source.n, ground.g)
    ]

    @named model = ODESystem(connections, t; systems=[resistor, inductor, source, ground])
    sys = structural_simplify(model)
    prob = ODAEProblem(sys, [inductor.i => 0.0], (0.0, 10.0))
    sol = solve(prob, Tsit5())

    # Plots.plot(sol; vars=[inductor.i, inductor.i])
    @test sol[inductor.i][end] ≈ 10 atol=1e-3
end

# RC with different voltage sources
@testset "RC with voltage sources" begin
    @named source_const = ConstantVoltage(V=10)
    @named source_sin = SineVoltage(offset=1, amplitude=10, frequency=2, start_time=0.5, phase=0)
    @named source_step = StepVoltage(offset=1, height=10, start_time=0.5)
    @named source_tri = TriangularVoltage(offset=1, start_time=0.5, amplitude=10, frequency=2)
    @named source_dsin = ExpSineVoltage(offset=1, amplitude=10, frequency=2, start_time=0.5, phase=0, damping=0.5)
    @named source_ramp = RampVoltage(offset=1, height=10, start_time=0.5, duration=1)
    sources = [source_const, source_sin, source_step, source_tri, source_dsin, source_ramp]

    @named resistor = Resistor(R=1)
    @named capacitor = Capacitor(C=1)
    @named ground = Ground()

    for source in sources
        connections = [
            connect(source.p, resistor.p)
            connect(resistor.n, capacitor.p)
            connect(capacitor.n, source.n, ground.g)
        ]

        @named model = ODESystem(connections, t; systems=[resistor, capacitor, source, ground])
        sys = structural_simplify(model)
        prob = ODAEProblem(sys, [capacitor.v => 0.0], (0.0, 10.0))
        @test_nowarn sol = solve(prob, Tsit5())
        @test_nowarn sol = solve(prob, Rodas4())

        # Plots.plot(sol; vars=[source.v, capacitor.v])
    end
end

# RL with different voltage sources
@testset "RL with voltage sources" begin
    @named source_const = ConstantVoltage(V=10)
    @named source_sin = SineVoltage(offset=1, amplitude=10, frequency=2, start_time=0.5, phase=0)
    @named source_step = StepVoltage(offset=1, height=10, start_time=0.5)
    @named source_tri = TriangularVoltage(offset=1, start_time=0.5, amplitude=10, frequency=2)
    @named source_dsin = ExpSineVoltage(offset=1, amplitude=10, frequency=2, start_time=0.5, phase=0, damping=0.5)
    @named source_ramp = RampVoltage(offset=1, height=10, start_time=0.5, duration=1)
    sources = [source_const, source_sin, source_step, source_tri, source_dsin, source_ramp]

    @named resistor = Resistor(R=1.0)
    @named inductor = Inductor(L=1.0)
    @named ground = Ground()

    for source in sources
        connections = [
            connect(source.p, resistor.p)
            connect(resistor.n, inductor.p)
            connect(inductor.n, source.n, ground.g)
        ]

        @named model = ODESystem(connections, t; systems=[resistor, inductor, source, ground])
        sys = structural_simplify(model)
        prob = ODAEProblem(sys, [inductor.i => 0.0], (0.0, 10.0))
        @test_nowarn sol = solve(prob, Tsit5())
        @test_nowarn sol = solve(prob, Rodas4())

        # Plots.plot(sol; vars=[source.i, inductor.i])
    end
end

# RC with different current sources
@testset "RC with current sources" begin
    @named source_const = ConstantCurrent(I=10)
    @named source_sin = SineCurrent(offset=1, amplitude=10, frequency=2, start_time=0.5, phase=0)
    @named source_step = StepCurrent(offset=1, height=10, start_time=0.5)
    @named source_tri = TriangularCurrent(offset=1, start_time=0.5, amplitude=10, frequency=2)
    @named source_dsin = ExpSineCurrent(offset=1, amplitude=10, frequency=2, start_time=0.5, phase=0, damping=0.5)
    @named source_ramp = RampCurrent(offset=1, height=10, start_time=0.5, duration=1)
    sources = [source_const, source_sin, source_step, source_tri, source_dsin, source_ramp]

    @named resistor = Resistor(R=1)
    @named capacitor = Capacitor(C=1)
    @named ground = Ground()

    for source in sources
        connections = [
            connect(source.p, resistor.p)
            connect(resistor.n, capacitor.p)
            connect(capacitor.n, source.n, ground.g)
        ]

        @named model = ODESystem(connections, t; systems=[resistor, capacitor, source, ground])
        sys = structural_simplify(model)
        prob = ODAEProblem(sys, [capacitor.v => 0.0], (0.0, 10.0))
        @test_nowarn sol = solve(prob, Tsit5())
        @test_broken sol = solve(prob, Rodas4())

        # Plots.plot(sol; vars=[source.v, capacitor.v])
    end
end

@testset "Integrator" begin
    R=1e3
    f=1
    Vin=5
    @named ground = Ground()
    @named R1 = Resistor(R=R)
    @named R2 = Resistor(R=100*R)
    @named C1 = Capacitor(C=1/(2 * pi * f * R))
    @named opamp = IdealOpAmp()
    @named square = SquareVoltage(amplitude=Vin)
    @named sensor = VoltageSensor()

    connections = [
        connect(square.p, R1.p)
        connect(R1.n, C1.n, R2.p, opamp.n1)
        connect(opamp.p2, C1.p, R2.n)
        connect(opamp.p1, ground.g, opamp.n2, square.n)
        connect(opamp.p2, sensor.p)
        connect(sensor.n, ground.g)
    ]
    @named model = ODESystem(connections, t, systems = [R1, R2, opamp, square, C1, ground, sensor])
    sys = structural_simplify(model)
    u0 = [
        C1.v => 0.0
        R1.v => 0.0
    ]
    prob = ODEProblem(sys, u0, (0, 100.0))
    sol = solve(prob, Rodas4())
    @test sol[opamp.v2] == sol[C1.v] # Not a great one however. Rely on the plot
    @test sol[opamp.p2.v] == sol[sensor.v]

    # plot(sol, vars=[sensor.v, square.v, C1.v])
end

@testset "Voltage function generators" begin
    st, o, h, f, A, et, ϕ, d, δ = 0.7, 1.25, 3, 2, 2.5, 2.5, π/4, 0.1, 0.0001

    @named res = Resistor(R=1)
    @named cap = Capacitor(C=1)
    @named ground = Ground()
    @named voltage_sensor = VoltageSensor()
    @named vstep = StepVoltage(start_time=st, offset=o, height=h)
    @named vsquare = SquareVoltage(offset=o, start_time=st, amplitude=A, frequency=f)
    @named vtri = TriangularVoltage(offset=o, start_time=st, amplitude=A, frequency=f)
    # @named vsawtooth = SawToothVoltage(amplitude=A, start_time=st, frequency=f, offset=o)
    @named vcosine = CosineVoltage(offset=o, amplitude=A, frequency=f, start_time=st, phase=ϕ)
    @named vdamped_sine = ExpSineVoltage(offset=o, amplitude=A, frequency=f, start_time=st, phase=ϕ, damping=d)
    @named vramp = RampVoltage(offset=o, start_time=st, duration=et-st, height=h)

    vsources = [vtri, vsquare, vstep, vcosine, vdamped_sine, vramp]
    waveforms(i, x) = getindex([o .+ (x .> st) .* _triangular_wave.(x, δ, f, A, st),
                                o .+ (x .> st) .* _square_wave.(x, δ, f, A, st),
                                o .+  _step.(x, δ, h, st),
                                # o .+ (x .> st). * _sawtooth_wave.(x, δ, f, A, st),
                                o .+ (x .> st) .* _cos_wave.(x, f, A, st, ϕ),
                                o .+ (x .> st) .* _damped_sine_wave.(x, f, A, st, ϕ, d),
                                o .+ _ramp.(x, δ, st, et, h)], i)
    for i in 1:length(vsources)
        vsource = vsources[i]
        # @info Symbolics.getname(vsource)
        eqs = [
            connect(vsource.p, voltage_sensor.p, res.p)
            connect(res.n, cap.p)
            connect(ground.g, voltage_sensor.n, vsource.n, cap.n)
        ]
        @named vmodel = ODESystem(eqs, t, systems = [voltage_sensor, res, cap, vsource, ground])
        vsys = structural_simplify(vmodel)

        u0 = [
            vsource.v => 1
            res.v => 1
        ]

        prob = ODAEProblem(vsys, u0, (0, 10.0))
        sol = solve(prob, dt=0.1, Tsit5())

        @test sol[vsource.v][1150:end] ≈ waveforms(i, sol.t)[1150:end] atol=1e-1
        # For visual inspection
        # plt = plot(sol; vars=[vsource.v])
        # savefig(plt, "test_voltage_$(Symbolics.getname(vsource))")
    end
end

@testset "Current function generators" begin
    st, o, h, f, A, et, ϕ, d, δ = 0.7, 1.25, 3, 2, 2.5, 2.5, π/4, 0.1, 0.0001

    @named ground = Ground()
    @named res = Resistor(R=1.0)
    @named cap = Capacitor(C=1)
    @named current_sensor = CurrentSensor()
    @named istep = StepCurrent(start_time=st, offset=o, height=h)
    @named isquare = SquareCurrent(offset=o, start_time=st, amplitude=A, frequency=f)
    @named itri = TriangularCurrent(offset=o, start_time=st, amplitude=A, frequency=f)
    # @named isawtooth = SawToothCurrent(amplitude=A, start_time=st, frequency=f, offset=o)
    @named icosine = CosineCurrent(offset=o, amplitude=A, frequency=f, start_time=st, phase=ϕ)
    @named idamped_sine = ExpSineCurrent(offset=o, amplitude=A, frequency=f, start_time=st, phase=ϕ, damping=d)
    @named iramp = RampCurrent(offset=o, start_time=st, duration=et-st, height=h)

    isources = [itri, isquare, istep, icosine, idamped_sine, iramp]
    waveforms(i, x) = getindex([o .+ (x .> st) .* _triangular_wave.(x, δ, f, A, st),
                                o .+ (x .> st) .* _square_wave.(x, δ, f, A, st),
                                o .+  _step.(x, δ, h, st),
                                # o .+ (x .> st). * _sawtooth_wave.(x, δ, f, A, st),
                                o .+ (x .> st) .* _cos_wave.(x, f, A, st, ϕ),
                                o .+ (x .> st) .* _damped_sine_wave.(x, f, A, st, ϕ, d),
                                o .+ _ramp.(x, δ, st, et, h)], i)

    for i in 1:length(isources)
        isource = isources[i]
        eqs = [
            connect(isource.p, current_sensor.n)
            connect(current_sensor.p, res.p)
            connect(res.n, cap.p)
            connect(isource.n, ground.g, cap.n)
        ]
        @named model = ODESystem(eqs, t, systems = [current_sensor, isource, res, cap, ground])
        isys = structural_simplify(model)

        u0 = [
            isource.i => 1.0
            res.v => 1.0
            cap.v => 0.0
        ]
        prob = ODAEProblem(isys, u0, (0, 10.0))
        sol = solve(prob, Tsit5())

        @test sol[isource.i][1150:end] ≈ waveforms(i, sol.t)[1150:end] atol=1e-1
        # For visual inspection
        # plt = plot(sol)
        # savefig(plt, "test_current_$(Symbolics.getname(isource))")
    end
end


text:
function mk_t!(C::GphysChannel, nx::Integer, ts::Int64)
  t = Array{Int64, 2}(undef, 2, 2)
  setindex!(t, one(Int64), 1)
  setindex!(t, nx, 2)
  setindex!(t, ts, 3)
  setindex!(t, zero(Int64), 4)
  setfield!(C, :t, t)
  return nothing
end

function mk_t(nx::Integer, ts::Int64)
  t = Array{Int64, 2}(undef, 2, 2)
  setindex!(t, one(Int64), 1)
  setindex!(t, nx, 2)
  setindex!(t, ts, 3)
  setindex!(t, zero(Int64), 4)
  return t
end


text:
# Rotation around the origin
rotate(d::EuclideanDomain{2}, θ) = rotation_map(θ).(d)

rotate(d::EuclideanDomain{3}, phi, theta, psi) = rotation_map(phi,theta,psi).(d)
# Rotation around a fixed center.
rotate(d::EuclideanDomain{2}, θ, center::SVector{T}) where {T} = (Translation(center) ∘ rotation_map(θ) ∘ Translation(-center)).(d)

rotate(d::EuclideanDomain{3}, phi, theta, psi, center::SVector{T}) where {T} = (Translation(center) ∘ rotation_map(phi,theta,psi) ∘ Translation(-center)).(d)

# Maps having to do with coordinate transforms.

"""
A Cartesion to Polar map. First dimension is interpreted as radial distance,
second as an angle.
The unit circle is mapped to the square `[-1,1]x[-1,1]`.
"""
struct CartToPolarMap{T} <: Map{SVector{2,T}}
end

CartToPolarMap() = CartToPolarMap{Float64}()

mapsize(m::CartToPolarMap) = (2,2)

applymap(map::CartToPolarMap{T}, x) where {T} =
    SVector{2,T}(sqrt(x[1]^2+x[2]^2)*2-1, atan(x[2],x[1])/pi)

function jacobian(m::CartToPolarMap{T}, x) where {T}
    d = sqrt(x[1]^2+x[2]^2)
    r = 1/(1+(x[2]/x[1])^2)
    SMatrix{2,2,T}(2*x[1]/d, 1/pi*r*(-x[2]/x[1]^2), 2*x[2]/d, 1/pi*r*1/x[1])
end

inverse(m::CartToPolarMap{T}) where {T} = PolarToCartMap{T}()
inverse(m::CartToPolarMap, x) = inverse(m)(x)

isreal(m::CartToPolarMap) = true

convert(::Type{Map{SVector{2,T}}}, ::CartToPolarMap) where {T} = CartToPolarMap{T}()

==(m1::CartToPolarMap, m2::CartToPolarMap) = true

"""
A Polar to Cartesian map. The angle is mapped to the second dimension,
radius to the first.
The square `[-1,1]x[-1,1]` is mapped to the unit circle.
"""
struct PolarToCartMap{T} <: Map{SVector{2,T}}
end

PolarToCartMap() = PolarToCartMap{Float64}()

mapsize(m::PolarToCartMap) = (2,2)

applymap(map::PolarToCartMap{T}, x) where {T} = SVector{2,T}((x[1]+1)/2*cos(pi*x[2]), (x[1]+1)/2*sin(pi*x[2]))

jacobian(m::PolarToCartMap{T}, x) where {T} =
    SMatrix{2,2,T}(cos(pi*x[2])/2, sin(pi*x[2])/2, -pi*(x[1]+1)/2*sin(pi*x[2]), pi*(x[1]+1)/2*cos(pi*x[2]))

inverse(m::PolarToCartMap{T}) where {T} = CartToPolarMap{T}()
inverse(m::PolarToCartMap, x) = inverse(m)(x)

isreal(m::PolarToCartMap) = true

convert(::Type{Map{SVector{2,T}}}, ::PolarToCartMap) where {T} = PolarToCartMap{T}()

==(m1::PolarToCartMap, m2::PolarToCartMap) = true

#############################
# Rotations around the origin
#############################

# Rotation in positive (counterclockwise) direction
# (Note: the SMatrix constructor expects the arguments column-first)
rotationmatrix(theta) = SMatrix{2,2}(cos(theta), sin(theta), -sin(theta), cos(theta))

# Rotation about X-axis (phi), Y-axis (theta) and Z-axis (psi)
# As above, the matrix is given column-by-column
rotationmatrix(phi,theta,psi) =
    SMatrix{3,3}(cos(theta)*cos(psi), cos(theta)*sin(psi), -sin(theta),
        -cos(phi)*sin(psi)+sin(phi)*sin(theta)*cos(psi), cos(phi)*cos(psi)+sin(phi)*sin(theta)*sin(psi), sin(phi)*cos(theta),
        sin(phi)*sin(psi)+cos(phi)*sin(theta)*cos(psi), -sin(phi)*cos(psi)+cos(phi)*sin(theta)*sin(psi), cos(phi)*cos(theta))

rotation_map(theta) = LinearMap(rotationmatrix(theta))

rotation_map(phi, theta, psi) = LinearMap(rotationmatrix(phi,theta,psi))


text:
using Test
using AcceleratedArrays
using SplitApplyCombine
using Dictionaries

@test isempty(setdiff(detect_ambiguities(Base, AcceleratedArrays, Dictionaries), detect_ambiguities(Base, Dictionaries) ))

@testset "AcceleratedArrays" begin
    include("Interval.jl")
    include("MaybeVector.jl")
    include("SingleVector.jl")

    include("HashIndex.jl")
    include("UniqueHashIndex.jl")
    include("SortIndex.jl")
    include("UniqueSortIndex.jl")
end


text:
# TODO should this have AbstractArray interface?
struct Memory
    v::Vector{UInt8}
end

Memory(v::AbstractVector{UInt8}) = Memory(Vector(v))

Memory(ℓ::Integer=65536) = Memory(zeros(UInt8, ℓ))

abstract type AbstractΠ end  # uses capital Π

# note that, somewhat confusingly, T is the type of the pointer, not the value being pointed to
struct Π{T<:Unsigned} <: AbstractΠ
    addr::T

    Π{T}(addr::Unsigned) where T<:Unsigned = new(addr)
end

Π(addr::T) where {T<:Unsigned} = Π{T}(addr)

const Π8 = Π{UInt8}  # zero page pointer
const Π16 = Π{UInt16}  # standard 6502 memory pointer

Π8(x::AbstractVector{UInt8}) = Π{UInt8}(x[1])

Π16(x::UInt8, y::UInt8) = Π{UInt16}(bincat(x, y))
Π16(x::AbstractVector{UInt8}) = Π16(x[1], x[2])

reset!(m::Memory) = (m.v .= zeros(UInt8, length(m.v)))

Base.fetch(m::Memory, idx::T) where T<:Unsigned = m.v[idx+one(T)]
Base.fetch(m::Memory, idx::AbstractVector{T}) where T<:Unsigned = m.v[idx+one(T)]

Base.getindex(m::Memory, idx) = m.v[idx + 0x0001]
Base.getindex(m::Memory, ptr::Π) = fetch(m, ptr.addr)

Base.setindex!(m::Memory, val::UInt8, idx::T) where T<:Unsigned = (m.v[idx+one(T)] = val)
Base.setindex!(m::Memory, val::UInt8, ptr::Π) = setindex!(m, val, ptr.addr)

function Base.setindex!(m::Memory, v::AbstractVector{UInt8},
                        idx::AbstractVector{T}) where T <: Integer
    m.v[idx .+ one(T)] = v
end

Base.view(m::Memory, idx) = view(m.v, idx + 0x0001)
Base.view(m::Memory, ptr::Π) = view(m, ptr.addr)

#==============================================================================================
    <referencing>
==============================================================================================#
dereference(ptr::Π, m::Memory) = m[ptr]
deref(ptr::Π, m::Memory) = dereference(ptr, m)
(↦)(ptr::Π, m::Memory) = dereference(ptr, m)  # this symbol is \mapsto

# note that the views always return SubArray, even for a single element
derefview(ptr::Π, m::Memory) = view(m, ptr)

dereference(::Type{UInt8}, ptr::Π, m::Memory)::UInt8 = dereference(ptr, Π)
deref(::Type{UInt8}, ptr::Π, m::Memory) = dereference(ptr, Π)

# these are for dereferencing length 2 arrays into single UInt16's
function dereference(::Type{UInt16}, ptr::Π{T}, m::Memory)::UInt16 where {T}
    least_sig = deref(ptr, m)
    most_sig = deref(ptr + one(T), m)
    0x0100 * most_sig + least_sig
end
deref(::Type{UInt16}, ptr::Π, m::Memory) = dereference(UInt16, ptr, m)
(⇾)(ptr::Π, m::Memory) = dereference(UInt16, ptr, m)  # this symbol is \rightarrowtriangle

# dereference arrays
function dereference(ptr::Π{T}, m::Memory, ℓ::Integer) where {T}
    start_idx = ptr.addr + one(T)
    end_idx = start_idx + ℓ - 1
    m.v[start_idx:end_idx]
end
deref(ptr::Π, m::Memory, ℓ::Integer) = dereference(ptr, m, ℓ)

# returns SubArray
function derefview(ptr::Π, m::Memory, ℓ::Integer)
    start_idx = ptr.addr + one(T)
    end_idx = start_idx + ℓ - 1
    view(m.v, start_idx:end_idx)
end

store!(m::Memory, ptr::Π, val::UInt8) = (m[ptr.addr] = val)
# this stores backwards from the supplied memory address
function storeback!(m::Memory, ptr::Π, val::UInt16)
    m[ptr] = UInt8(val >> 8)
    m[ptr - 0x01] = UInt8(0x00ff & val)
    val
end

Base.:+(ptr1::Π{T}, ptr2::Π{T}) where {T} = Π{T}(ptr1.addr + ptr2.addr)
Base.:-(ptr1::Π{T}, ptr2::Π{T}) where {T} = Π{T}(ptr1.addr - ptr2.addr)

Base.:+(ptr::Π{T}, val::Unsigned) where {T} = Π{T}(ptr.addr + convert(T, val))
Base.:-(ptr::Π{T}, val::Unsigned) where {T} = Π{T}(ptr.addr - convert(T, val))

#==============================================================================================
    </referencing>
==============================================================================================#


text:
__precompile__(true)

module DistributedArrays

using Primes
using Primes: factor

importall Base
import Base.Callable
import Base.BLAS: axpy!

export (.+), (.-), (.*), (./), (.%), (.<<), (.>>), div, mod, rem, (&), (|), ($)
export DArray, SubDArray, SubOrDArray, @DArray
export dzeros, dones, dfill, drand, drandn, distribute, localpart, localindexes, ppeval, samedist
export close, darray_closeall

include("core.jl")
include("serialize.jl")
include("mapreduce.jl")
include("linalg.jl")
include("sort.jl")

end # module


text:
###############################################################################
#
# The Stopping structure eases the implementation of algorithms and the
# stopping criterion.
#
# The following examples illustre the various possibilities offered by Stopping
#
###############################################################################

using Test, NLPModels, Stopping

#printstyled("How to State ")
include("howtostate.jl")
#printstyled("passed ✓ \n", color = :green)
#printstyled("How to State for NLP ")
include("howtostate-nlp.jl")
#printstyled("passed ✓ \n", color = :green)
#printstyled("How to Stop ")
include("howtostop.jl")
#printstyled("passed ✓ \n", color = :green)
#printstyled("How to Stop II ")
include("howtostop-2.jl")
#printstyled("passed ✓ \n", color = :green)
#printstyled("How to Stop for NLP ")
include("howtostop-nlp.jl")
#printstyled("passed ✓ \n", color = :green)


text:
module ProximalBundleMethod

import LinearAlgebra
import Printf
import Convex
import SCS
# import Gurobi

const dot = LinearAlgebra.dot
const norm = LinearAlgebra.norm
const minimize = Convex.minimize

include("solver.jl")

end # module


text:
using HDF5
using Random
include("./karger.jl")
using .Karger

# USAGE: julia src/julia/calculate_potential.jl filename N
# where N is the number of samples to use and filename just the name (without path or .h5).
# The file has to be an HDF5 file with fields n, edges, weights and seeds. The output
# will be written in a file with the same name in the "potentials" directory

Random.seed!(0)

graph_path = "results/graphs/" * ARGS[1] * ".h5"
result_path = "results/karger_potentials/" * ARGS[1] * ".h5"

n = h5read(graph_path, "n")
edges = h5read(graph_path, "edges")
weights = h5read(graph_path, "weights")
seeds = h5read(graph_path, "seeds")
g = Graph(n, transpose(edges) .+ 1, weights)
mkpath(dirname(result_path))
cp(graph_path, result_path, force=true)
pots = potential(g, parse(Int, ARGS[2]), seeds)
segmentation = Vector{Int}(undef, n)
max_pots = zeros(n)
for (key, pot) in pots
    h5write(result_path, "potential/" * string(key), pot)
    for i in 1:n
        if pot[i] > max_pots[i]
            segmentation[i] = key
            max_pots[i] = pot[i]
        end
    end
end
h5write(result_path, "segmentation", segmentation)


text:
struct BoundingBox
    xmin::Real
    xmax::Real
    ymin::Real
    ymax::Real
end

const colorplaces = [:r,:g,:b]

iterate(c::ColorTypes.RGB) = (getfield(c,first(colorplaces)),2)
iterate(c::ColorTypes.RGB,i) = i > 3 ? nothing : (getfield(c,colorplaces[i]),i+1)

iterate(c::NamedColor) = iterate(c.color)
iterate(c::NamedColor,i) = iterate(c.color,i)

+(P::Vec2,Q::Vec2) = Vec2(P.x+Q.x,P.y+Q.y)
-(Q::Vec2) = Vec2(-Q.x,-Q.y)
-(P::Vec2,Q::Vec2) = P + (-Q)
*(c::Real,P::Vec2) = Vec2(c*P.x,c*P.y)

function boundingbox(path::Union{Path2D,Polygon2D})
    # Coordinates are returned in the order xmin, xmax, ymax, ymin
    return BoundingBox(minimum(pt.x for pt in path.points),
                       maximum(pt.x for pt in path.points),
                       minimum(pt.y for pt in path.points),
                       maximum(pt.y for pt in path.points))
end

function boundingbox(c::Circle2D)
    return BoundingBox(c.center.x - c.radius,
                       c.center.x + c.radius,
                       c.center.y - c.radius,
                       c.center.y + c.radius)
end

function boundingbox(p::Point2D)
    return BoundingBox(p.P.x,p.P.x,p.P.y,p.P.y)
end

function boundingbox(A::Array{BoundingBox,1};border=0.05)
    a = minimum(bb.xmin for bb in A)
    b = maximum(bb.xmax for bb in A)
    c = maximum(bb.ymax for bb in A)
    d = minimum(bb.ymin for bb in A)
    if a == b || c == d
        error("Bounding box has empty interior")
    end
    margin = border*min(b-a,c-d)
    return BoundingBox(a - margin,
                       b + margin,
                       d - margin,
                       c + margin)
end

function isclip(G::GraphicElement)
    if isa(G,Circle2D) || isa(G,Polygon2D)
        G.clip
    else
        false
    end
end

function boundingbox(P::Plot2D;kwargs...)
    k = something(findlast(isclip,P.elements),1)
    boundingbox(map(boundingbox,P.elements[k:end]);kwargs...)
end

function boundingbox(P::PixelMap)
    BoundingBox(P.lowerleft[1],P.upperright[1],
                P.lowerleft[2],P.upperright[2])
end

function aspectratio(bb::BoundingBox)
    return (bb.xmax - bb.xmin)/(bb.ymax - bb.ymin)
end

function boundingbox(L::Label2D)
    return BoundingBox(L.location.x,
                       L.location.x,
                       L.location.y,
                       L.location.y)
end

lwcorrect(x) = x == 0 ? 0.7 : x

function arrowhead(P::Path2D,bb::BoundingBox)
    A,B = P.points[end-1:end]
    z = (B-A).x + im*((B-A).y)
    θ = angle(z)
    ψ = 16π/180
    l = 5e-3*P.arrow.size*min(bb.xmax-bb.xmin,bb.ymax-bb.ymin)
    arrowloc = P.arrow.position
    tip = A + arrowloc*(B-A)
    C = tip-l*Vec2(cos(θ+ψ),sin(θ+ψ))
    D = tip-l*Vec2(cos(θ-ψ),sin(θ-ψ))
    return Polygon2D([tip,C,D,tip],
                      pen=NoPen(),
                      fillpen=P.pen)
end

# Note: more Cairo functionality moved to __init__ function in init.jl


text:
using TimeZones: TimeZones, ZonedDateTime, @tz_str
export ZonedDateTime, @tz_str
using Dates: Dates, Date
export Dates, Date

@system Calendar begin
    init ~ preserve::datetime(extern, parameter)
    last => nothing ~ preserve::datetime(extern, parameter, optional)
    time(t0=init, t=context.clock.time) => t0 + convert(Cropbox.Dates.Second, t) ~ track::datetime
    date(time) => Cropbox.Dates.Date(time) ~ track::date
    step(context.clock.step) ~ preserve(u"hr")
    stop(time, last) => begin
        isnothing(last) ? false : (time >= last)
    end ~ flag
    count(init, last, step) => begin
        if isnothing(last)
            nothing
        else
            # number of update!() required to reach `last` time
            (last - init) / step
        end
    end ~ preserve::int(round, optional)
end

export Calendar


text:
int main() {
  int b = 3;
  int c = b.length;
  return 0;
}

max_stars_repo_name: fredefox/javalette

max_issues_repo_name: fredefox/javalette

max_forks_repo_name: fredefox/javalette


text:
# Read a set of vectors stored in the bvec format (int + n * float)
# The function returns a set of output uint8 vector (one vector per column)
#
# Syntax:
#   v = fvecs_read (filename)        -> read all vectors
#   v = fvecs_read (filename, n)      -> read n vectors
#   v = fvecs_read (filename, [a b]) -> read the vectors from a to b (indices starts from 1)

function fvecs_read(
  bounds::UnitRange,
  filename::AbstractString="./data/deep/deep10M.fvecs")

  @assert bounds.start >= 1

  # open the file and count the number of descriptors
  open(filename, "r") do fid

    # Read the vector size
    d = read(fid, Int32, 1)
    @assert length(d) == 1
    @show vecsizeof = 1 * 4 + d[1] * 4

    # Get the number of vectrors
    seekend(fid)
    vecnum = position(fid) / vecsizeof

    # compute the number of vectors that are really read and go in starting positions
    @show n = bounds.stop - bounds.start + 1
    seekstart(fid)
    skip(fid, (bounds.start - 1) * vecsizeof)

    # read n vectors
    v = read(fid, Float32, (d[1] + 1) * n)
    v = reshape(v, d[1] + 1, n)

    # Check if the first column (dimension of the vectors) is correct
    @assert sum( v[1,2:end] .== v[1, 1]) == n - 1
    v = v[2:end, :]

    return v

  end
end

function fvecs_read(n::Integer, filename::AbstractString)
  return fvecs_read(1:n,filename)
end

function fvecs_read(n::Integer)
  return fvecs_read(1:n)
end

max_issues_repo_name: xiaosu-zhu/local-search-quantization

max_forks_repo_name: xiaosu-zhu/local-search-quantization


text:
using ODEInterfaceDiffEq, Test

callback_f = function (du,u,p,t)
  du[1] = u[2]
  du[2] = -9.81
end

condtion= function (u,t,integrator) # Event when event_f(u,t,k) == 0
  u[1]
end

affect! = nothing
affect_neg! = function (integrator)
  integrator.u[2] = -integrator.u[2]
end

callback = ContinuousCallback(condtion,affect!,affect_neg!)

u0 = [50.0,0.0]
tspan = (0.0,25.0)
prob = ODEProblem(callback_f,u0,tspan)

sol = solve(prob,dopri5(),callback=callback,dtmax=0.5)
@test sol(4.0)[1] > 0
sol = solve(prob,dopri5(),callback=callback,save_everystep=true)
@test sol(4.0)[1] > -1e-12


text:
module interaction

using ..system
using LinearAlgebra

export GreenTensor,G_ij,Gamma_ij,Omega_ij

"""
    interaction.F(ri::Vector, rj::Vector, µi::Vector, µj::Vector)

General F function for arbitrary positions and dipole orientations.

Arguments:
* ri: Position of first spin
* rj: Position of second spin
* µi: Dipole orientation of first spin.
* µj: Dipole orientation of second spin.
"""
function F(ri::Vector, rj::Vector, µi::Vector, µj::Vector)
    rij = ri - rj
    rij_norm = norm(rij)
    rijn = rij./rij_norm
    μi_ = normalize(μi)
    μj_ = normalize(μj)
    T = float(promote_type(eltype(rij),eltype(μi_),eltype(μj_)))
    if rij_norm == 0
        T(2/3.)
    else
        ξ = 2π*rij_norm
        T(dot(µi_, µj_)*(sin(ξ)/ξ + cos(ξ)/ξ^2 - sin(ξ)/ξ^3) + dot(µi_, rijn)*dot(rijn, µj_)*(-sin(ξ)/ξ - 3*cos(ξ)/ξ^2 + 3*sin(ξ)/ξ^3))
    end
end

"""
    interaction.G(ri::Vector, rj::Vector, µi::Vector, µj::Vector)

General G function for arbitrary positions and dipole orientations.

Arguments:
* ri: Position of first spin
* rj: Position of second spin
* µi: Dipole orientation of first spin.
* µj: Dipole orientation of second spin.
"""
function G(ri::Vector, rj::Vector, µi::Vector, µj::Vector)
    rij = ri - rj
    rij_norm = norm(rij)
    rijn = rij./rij_norm
    μi_ = normalize(μi)
    μj_ = normalize(μj)
    T = float(promote_type(eltype(rij),eltype(μi_),eltype(μj_)))
    if rij_norm == 0
        zero(T)
    else
        ξ = 2π*rij_norm
        T(dot(µi_, µj_)*(-cos(ξ)/ξ + sin(ξ)/ξ^2 + cos(ξ)/ξ^3) + dot(µi_, rijn)*dot(rijn, µj_)*(cos(ξ)/ξ - 3*sin(ξ)/ξ^2 - 3*cos(ξ)/ξ^3))
    end
end

"""
    interaction.Omega(ri::Vector, rj::Vector, µi::Vector, µj::Vector, γi::Real=1, γj::Real=1)

Arguments:
* ri: Position of first spin
* rj: Position of second spin
* µi: Dipole orientation of first spin.
* µj: Dipole orientation of second spin.
* γi: Decay rate of first spin.
* γj: Decay rate of second spin.

Note that the dipole moments `μi` and `μj` are normalized internally. To account
for dipole moments with different lengths you need to scale the decay rates
`γi` and `γj`, respectively.
"""
function Omega(ri::Vector, rj::Vector, µi::Vector, µj::Vector, γi::Real=1, γj::Real=1)
    return 0.75*sqrt(γi*γj)*G(ri, rj, µi, µj)
end

"""
    interaction.Gamma(ri::Vector, rj::Vector, µi::Vector, µj::Vector, γi::Real=1, γj::Real=1)

Arguments:
* ri: Position of first spin
* rj: Position of second spin
* µi: Dipole orientation of first spin.
* µj: Dipole orientation of second spin.
* γi: Decay rate of first spin.
* γj: Decay rate of second spin.

Note that the dipole moments `μi` and `μj` are normalized internally. To account
for dipole moments with different lengths you need to scale the decay rates
`γi` and `γj`, respectively.
"""
function Gamma(ri::Vector, rj::Vector, µi::Vector, µj::Vector, γi::Real=1, γj::Real=1)
    return 1.5*sqrt(γi*γj)*F(ri, rj, µi, µj)
end

"""
    interaction.OmegaMatrix(S::SpinCollection)

Matrix of the dipole-dipole interaction for a given SpinCollection.
"""
function OmegaMatrix(S::system.SpinCollection)
    spins = S.spins
    mu = S.polarizations
    gamma = S.gammas
    N = length(spins)
    Ω = zeros(Float64, N, N)
    for i=1:N, j=1:N
        if i==j
            continue
        end
        Ω[i,j] = interaction.Omega(spins[i].position, spins[j].position, mu[i], mu[j], gamma[i], gamma[j])
    end
    return Ω
end

"""
    interaction.GammaMatrix(S::SpinCollection)

Matrix of the collective decay rate for a given SpinCollection.
"""
function GammaMatrix(S::system.SpinCollection)
    spins = S.spins
    mu = S.polarizations
    gamma = S.gammas
    N = length(spins)
    return [
        interaction.Gamma(spins[i].position, spins[j].position, mu[i], mu[j], gamma[i], gamma[j])
        for i=1:N, j=1:N
    ]
end

"""
    GreenTensor(r::Vector, k::Number=2π)

Calculate the Green's Tensor at position r for wave number k defined by

```math
G = e^{ikr}\\Big[\\left(\\frac{1}{kr} + \\frac{i}{(kr)^2} - \\frac{1}{(kr)^3}\\right)*I -
    \\textbf{r}\\textbrf{r}^T\\left(\\frac{1}{kr} + \\frac{3i}{(kr)^2} - \\frac{3}{(kr)^3}\\right)\\Big]
```

Choosing `k=2π` corresponds to the position `r` being given in units of the
wavelength associated with the dipole transition.

Returns a 3×3 complex Matrix.
"""
function GreenTensor(r::Vector{<:Number},k::Real=2π)
    n = norm(r)
    rn = r./n
    return exp(im*k*n)*(
        (1/(k*n) + im/(k*n)^2 - 1/(k*n)^3).*Matrix(I,3,3) +
        -(1/(k*n) + 3im/(k*n)^2 - 3/(k*n)^3).*(rn*rn')
    )
end

"""
    Gamma_ij(r1::Vector,r2::Vector,μ1::Vector,μ2::Vector;gamma1=1,gamma2=1,k0=2π)

From the field overlap with the Green's tensor, compute the mutual decay rate as

```math
Γ_{ij} = \\frac{3}{2} μᵢ* ⋅ ℑ(G) ⋅ μⱼ.
```

"""
function Gamma_ij(r1::Vector{<:Real},r2::Vector{<:Real},μ1::Vector{<:Number},μ2::Vector{<:Number};k0=2π,kwargs...)
    T = float(promote_type(eltype(r1),eltype(r2),eltype(μ1),eltype(μ2),typeof(k0)))
    if (r1 == r2) && (μ1 == μ2)
        return one(T)
    elseif (r1 == r2) && (μ1 != μ2)
        return zero(T)
    else
        G_im = imag(GreenTensor(r1 - r2, k0))
        return 1.5*dot(μ1, G_im, μ2)::T
    end
end

"""
    Omega_ij(r1::Vector,r2::Vector,μ1::Vector,μ2::Vector;gamma1=1,gamma2=1,k0=2π)

From the field overlap with the Green's tensor, compute the energy shifts as

```math
Ω_{ij} = -\\frac{3}{4} μᵢ* ⋅ ℜ(G) ⋅ μⱼ.
```

"""
function Omega_ij(r1::Vector{<:Real},r2::Vector{<:Real},μ1::Vector{<:Number},μ2::Vector{<:Number};k0=2π,kwargs...)
    T = float(promote_type(eltype(r1),eltype(r2),eltype(μ1),eltype(μ2),typeof(k0)))
    if r1 == r2
        return zero(T)
    else
        G_re = real(GreenTensor(r1 - r2, k0))
        return -0.75*dot(μ1,G_re,μ2)::T
    end
end

end # module


text:
using Test
using JSExpr
using JSExpr: @crawl, deparse, JSTerminal, JSString

@testset "call" begin
    @testset "call builtin functions" begin
        @test @js(console.log("foo")) == JSString("""console.log("foo")""")
    end
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    CoDaDistance

Distance between compositions in Aitchison geometry.
"""
struct CoDaDistance <: Metric end

(::CoDaDistance)(x, y) = distance(x, y)

result_type(::CoDaDistance, x, y) = Float64


text:
G = Graph(10)
A = 100

for i in 1:A
    new = sample(collect(nodes(G)), 2, replace=true)
    add_arc!(G, new[1], new[2])
end

u = 1
@test in_neighbor(G, u) == keys(badj(G)[u])
@test out_neighbor(G, u) == keys(fadj(G)[u])

@test nodes(G) == keys(fadj(G))
@test length(nodes(G)) == G.nn

tes_out_neighbor = out_neighbor(G, u)
out_neighbor_arc = Arc[]
for v in tes_out_neighbor
    auv = collect(arcs(G, u, v))
    append!(out_neighbor_arc, auv)
end
@test out_neighbor_arc == collect(arcs(G, u, :))


text:
export ruleVBLogitOut, ruleVBLogitIn1, ruleVBLogitXi

function ruleVBLogitOut(marg_out::Any,
                        marg_in1::ProbabilityDistribution{Univariate},
                        marg_xi::ProbabilityDistribution{Univariate})

    return Message(Univariate, Bernoulli, p=logisticSigmoid(unsafeMean(marg_in1)))
end

function ruleVBLogitIn1(marg_out::ProbabilityDistribution{Univariate},
                        marg_in1::Any,
                        marg_xi::ProbabilityDistribution{Univariate})

    xi_hat = unsafeMode(marg_xi)

    return Message(Univariate, GaussianWeightedMeanPrecision, xi=unsafeMean(marg_out) - 0.5, w=2*logisticLambda(xi_hat))
end

function ruleVBLogitXi(marg_out::ProbabilityDistribution{Univariate},
                       marg_in1::ProbabilityDistribution{Univariate},
                       marg_xi::Any)

    return Message(Univariate, Function, mode=sqrt(unsafeMean(marg_in1)^2 + unsafeCov(marg_in1)))
end


text:
module Shufflers

function Independent(popSize)

    function shuffle()
    end
    (
     popSize=popSize,
     stateType=nothing,
     init=nothing,
     shuffle=shuffle,
     callEvery=0
    )
end

end


text:
module StackPointers

using VectorizationBase

export StackPointer, stack_pointer_call

struct StackPointer # default to Float64
    ptr::Ptr{Float64}
end
@inline Base.pointer(s::StackPointer) = s.ptr
@inline Base.pointer(s::StackPointer, ::Type{T}) where {T} = Base.unsafe_convert(Ptr{T}, s.ptr)
@inline Base.pointer(s::StackPointer, ::Type{Float64}) = s.ptr

@inline Base.convert(::Type{Ptr{T}}, s::StackPointer) where {T} = pointer(s, T)
@inline Base.unsafe_convert(::Type{Ptr{T}}, s::StackPointer) where {T} = pointer(s, T)

@inline Base.:+(sp::StackPointer, i::Integer) = StackPointer(gep(sp.ptr, i))
@inline Base.:+(sp::StackPointer, i::Integer...) = StackPointer(gep(sp.ptr, +(i...)))
@inline Base.:+(i::Integer, sp::StackPointer) = StackPointer(gep(sp.ptr, i))
@inline Base.:-(sp::StackPointer, i::Integer) = StackPointer(gep(sp.ptr, -i))

@inline VectorizationBase.align(s::StackPointer) = StackPointer(VectorizationBase.align(s.ptr))

"""
To support using StackPointers, overload stack_pointer_call.
For example, in PaddedMatrices (which defines an appropriate `similar` method):

@inline function StackPointers.stack_pointer_call(::typeof(*), sp::StackPointer, A::AbstractMatrix, B::AbstractVecOrMat)
    (sp, C) = similar(sp, A, B)
    sp, mul!(C, A, B)
end

"""
@inline stack_pointer_call(f::F, sp::StackPointer, args::Vararg{<:Any,N}) where {F,N} = (sp, f(args...))

# function extract_func_sym(f::Expr)::Symbol
#     if f.head === :(.)
#         f.args[2].value
#     elseif f.head === :curly
#         ff = first(f.args)
#         ff isa Symbol ? ff : ff.args[2].value
#     end
# end

# include("precompile.jl")
# _precompile_()

end # module


text:
# Classical Multidimensional Scaling

## convert Gram matrix to Distance matrix

function gram2dmat!(D::AbstractMatrix{DT}, G::AbstractMatrix) where {DT<:Real}
    # argument checking
    m = size(G, 1)
    n = size(G, 2)
    m == n || error("D should be a square matrix.")
    size(D) == (m, n) ||
        throw(DimensionMismatch("Sizes of D and G do not match."))

    # implementation
    for j = 1:n
        for i = 1:j-1
            @inbounds D[i,j] = D[j,i]
        end
        D[j,j] = zero(DT)
        for i = j+1:n
            @inbounds D[i,j] = sqrt(G[i,i] + G[j,j] - 2 * G[i,j])
        end
    end
    return D
end

gram2dmat(G::AbstractMatrix{T}) where {T<:Real} = gram2dmat!(similar(G, T), G)

## convert Distance matrix to Gram matrix

function dmat2gram!(G::AbstractMatrix{GT}, D::AbstractMatrix) where GT
    # argument checking
    n = LinearAlgebra.checksquare(D)
    size(G) == (n, n) ||
        throw(DimensionMismatch("Sizes of G and D do not match."))

    # implementation
    u = zeros(GT, n)
    s = 0.0
    for j = 1:n
        s += (u[j] = sum(abs2, view(D,:,j)) / n)
    end
    s /= n

    for j = 1:n
        for i = 1:j-1
            @inbounds G[i,j] = G[j,i]
        end
        for i = j:n
            @inbounds G[i,j] = (u[i] + u[j] - abs2(D[i,j]) - s) / 2
        end
    end
    return G
end

momenttype(T) = typeof((zero(T) * zero(T) + zero(T) * zero(T))/ 2)
dmat2gram(D::AbstractMatrix{T}) where {T<:Real} = dmat2gram!(similar(D, momenttype(T)), D)

## Classical MDS

"""Classical MDS type"""
struct MDS{T<:Real}
    d::Real                  # original dimension
    X::AbstractMatrix{T}     # fitted data, X (d x n)
    λ::AbstractVector{T}     # sqrt. eigenvalues in feature space, √λ (k x 1)
    U::AbstractMatrix{T}     # eigenvectors in feature space, U (n x k)
end

## properties

indim(M::MDS) = M.d
outdim(M::MDS) = size(M.U,2)

projection(M::MDS) = M.U
eigvals(M::MDS) = M.λ

## use

"""Calculate out-of-sample multidimensional scaling transformation"""
function transform(M::MDS{T}, x::AbstractVector{T}; distances=false) where {T<:Real}
    d = if isnan(M.d) # model has only distance matrix
        @assert distances "Cannot transform points if model was fitted with a distance matrix. Use point distances."
        size(x, 1) != size(M.X, 1) && throw(
            DimensionMismatch("Point distances should be calculated to all original points"))
        x
    else
        if distances
            size(x, 1) != size(M.X, 2) && throw(
                DimensionMismatch("Point distances should be calculated to all original points."))
            x
        else
            size(x, 1) != size(M.X, 1) && throw(
                DimensionMismatch("Points and original data must have same dimensionality."))
            pairwise((x,y)->norm(x-y), M.X, x)
        end
    end

    # get distance matrix
    D = isnan(M.d) ? M.X  : pairwise((x,y)->norm(x-y), M.X)
    d = d.^2

    # b = 0.5*(ones(n,n)*d./n - d + D*ones(n,1)./n - ones(n,n)*D*ones(n,1)./n^2)
    mD = mean(D.^2, dims=2)
    b = (d  .- mean(d, dims=1) .- mD .+ mean(mD)) ./ -2

    # sqrt(λ)⁻¹U'b
    λ = vcat(M.λ, zeros(T, outdim(M) - length(M.λ)))
    return M.U' * b ./ sqrt.(λ)
end

"""Calculate multidimensional scaling transformation"""
function transform(M::MDS{T}) where {T<:Real}
    # if there are non-positive missing eigval then padd with zeros
    λ = vcat(M.λ, zeros(T, outdim(M) - length(M.λ)))
    return diagm(0=>sqrt.(λ)) * M.U'
end

## show

function Base.show(io::IO, M::MDS)
    print(io, "Classical MDS(indim = $(indim(M)), outdim = $(outdim(M)))")
end

## interface functions

function fit(::Type{MDS}, X::AbstractMatrix{T};
             maxoutdim::Int = size(X,1)-1,
             distances=false) where T<:Real

    # get distance matrix and space dimension
    D, d = if !distances
        pairwise((x,y)->norm(x-y), X), size(X,1)
    else
        X, NaN
    end
    G = dmat2gram(D)

    n = size(D, 1)
    m = min(maxoutdim, n) #Actual number of eigenpairs wanted

    E = eigen!(Hermitian(G))

    #Sometimes dmat2gram produces a negative definite matrix, and the sign just
    #needs to be flipped. The heuristic to check for this robustly is to check
    #if there is a negative eigenvalue of magnitude larger than the largest
    #positive eigenvalue, and flip the sign of eigenvalues if necessary.
    mineig, maxeig = extrema(E.values)
    if mineig < 0 && abs(mineig) > abs(maxeig)
        #do flip
        ord = sortperm(E.values)
        λ = -E.values[ord[1:m]]
    else
        ord = sortperm(E.values; rev=true)
        λ = E.values[ord[1:m]]
    end

    for i in 1:m
        if λ[i] <= 0
            #Keeping all remaining eigenpairs would not change solution (if 0)
            #or make the answer _worse_ (if <0).
            #The least squares solution would want to throw all these away.
            @warn("Gramian has only $(i-1) positive eigenvalue(s)")
            m = i-1
            ord = ord[1:m]
            λ = λ[1:m]
            break
        end
    end

    #Check if the last considered eigenvalue is degenerate
    if m>0
        nevalsmore = sum(abs.(E.values[ord[m+1:end]] .- λ[m]) .< n*eps())
        nevals = sum(abs.(E.values .- λ[m]) .< n*eps())
        if nevalsmore > 1
            @warn("The last eigenpair is degenerate with $(nevals-1) others; $nevalsmore were ignored. Answer is not unique")
        end
    end
    U = E.vectors[:, ord[1:m]]

    #Add trailing zero coordinates if dimension of embedding space (p) exceeds
    #number of eigenpairs used (m)
    if m < maxoutdim
        U = [U zeros(T, n, maxoutdim-m)]
    end

    return MDS(d, X, λ, U)
end

@deprecate classical_mds(D::AbstractMatrix, p::Int) transform(fit(MDS, D, maxoutdim=p, distances=true))

function stress(M::MDS)
    # calculate distances if original data was stored
    DX = isnan(M.d) ? M.X : pairwise((x,y)->norm(x-y), M.X)
    DY = pairwise((x,y)->norm(x-y), transform(M))
    n = size(DX,1)
    return sqrt(2*sum((DX - DY).^2)/sum(DX.^2));
end


text:
module Electre

import GreyDecision.GreyNumbers: GreyNumber
import GreyDecision.Utility:
    unitize, normalize, weightize, apply_columns, reverseminmax, rowmaxs

function electre(
    decisionMat::Array{GreyNumber{T},2},
    weights::Array{Float64,1},
    fns::Array{Function,1},
) where {T<:Real}

    w = unitize(weights)

    nalternatives, ncriteria = size(decisionMat)

    normalizedMat = normalize(decisionMat)
    weightednormalizedMat = weightize(normalizedMat, w)

    fitnessTable = []
    nonfitnessTable = []

    for i = 1:nalternatives
        for j = 1:nalternatives
            if i != j
                betterlist = []
                worselist = []
                for h = 1:ncriteria
                    if fns[h] == maximum
                        if weightednormalizedMat[i, h] >= weightednormalizedMat[j, h]
                            push!(betterlist, h)
                        else
                            push!(worselist, h)
                        end
                    else
                        if weightednormalizedMat[i, h] <= weightednormalizedMat[j, h]
                            push!(betterlist, h)
                        else
                            push!(worselist, h)
                        end
                    end
                end
                fitnesstableelement = Dict(:i => i, :j => j, :set => betterlist)
                push!(fitnessTable, fitnesstableelement)

                nonfitnesstableelement = Dict(:i => i, :j => j, :set => worselist)
                push!(nonfitnessTable, nonfitnesstableelement)
            end
        end
    end

    fitnessmatrix = Array{GreyNumber{Float64}, 2}(undef, nalternatives, nalternatives)
    nonfitnessmatrix = Array{GreyNumber{Float64}, 2}(undef, nalternatives, nalternatives)

    for elements in fitnessTable
        i = elements[:i]
        j = elements[:j]
        elemlist = elements[:set]

        CC = sum(w[elemlist])
        fitnessmatrix[i, j] = GreyNumber(CC, CC)
    end

    nonfitnessmatrix = Array{GreyNumber{Float64}, 2}(undef, nalternatives, nalternatives)
    for elements in nonfitnessTable
        i = elements[:i]
        j = elements[:j]
        elemlist = elements[:set]

        r_ik = weightednormalizedMat[i, elemlist]
        r_jk = weightednormalizedMat[j, elemlist]
        r_ik_full = weightednormalizedMat[i, :]
        r_jk_full = weightednormalizedMat[j, :]

        nom = maximum(abs.(r_ik - r_jk))
        dom = maximum(abs.(r_ik_full - r_jk_full))
        nonfitnessmatrix[i, j] = nom / dom
    end

    C = Array{GreyNumber{Float64}, 1}(undef, nalternatives)
    D = Array{GreyNumber{Float64}, 1}(undef, nalternatives)
    for i = 1:nalternatives
        C[i] = sum(fitnessmatrix[i, :]) - sum(fitnessmatrix[:, i])
        D[i] = sum(nonfitnessmatrix[i, :]) - sum(nonfitnessmatrix[:, i])
    end

    best_C_index = sortperm(C) |> last
    best_D_index = sortperm(D) |> first
    best = nothing

    if best_C_index == best_D_index
        best = (best_C_index,)
    else
        best = (best_C_index, best_D_index)
    end

    result = Dict(
        :C => C,
        :D => D,
        :bestindex => best
    )

    return result
end
end # end of module Electre


text:
# Test conversion between representations

using Test, CIJ
import StaticArrays

@testset "Conversion" begin
    @testset "Stiffness-Compliance" begin
        let c = CIJ.ol()[1], atol=1e-12, c′ = copy(c), ca = Array(c), ca′ = copy(ca)
            # ECs
            s = C2S(c)
            @test s[1,1] ≈  1.16918e-8 atol=atol
            @test s[2,2] ≈  2.03838e-8 atol=atol
            @test s[3,3] ≈  1.70592e-8 atol=atol
            @test s[1,2] ≈ -3.04176e-9 atol=atol
            @test s[1,3] ≈ -2.58468e-9 atol=atol
            @test s[2,3] ≈ -5.77166e-9 atol=atol
            for i in 4:6
                @test s[i,i] ≈ 1/c[i,i]
            end
            @test all(isapprox.([vec(s[1:3,4:6]); s[4,5:6]; s[5,6]], 0.0, atol=atol))
            @test C2S(S2C(c)) ≈ c
            @test C2S!(c′) ≈ s
            @test typeof(c) == typeof(s)
            # Arrays
            ca = Array(c)
            sa = C2S(ca)
            @test typeof(ca) == typeof(sa)
            @test S2C(C2S(ca)) ≈ ca
            @test sa ≈ s
            @test C2S!(ca′) == sa
            @test S2C!(ca′) ≈ c
        end
    end
end

@testset "Cij/Cijkl" begin
    let vp = rand(), vs = vp/1.7, ϵ = γ = δ = max(rand(), 0.01), (A, B, Γ) = 360 .* rand(3),
            c = CIJ.rot3(CIJ.thom(vp, vs, ϵ, γ, δ), A, B, Γ), ca = Array(c)
        c4 = cijkl(c)
        @test c4 isa StaticArrays.MArray
        @test size(c4) == (3, 3, 3, 3)
        @test c4[1,1,1,1] == c[1,1]
        @test c4[2,2,2,2] == c[2,2]
        @test c4[3,3,3,3] == c[3,3]
        @test c4[2,3,2,3] == c[4,4]
        @test c4[1,3,1,3] == c[5,5]
        @test c4[1,2,1,2] == c[6,6]
        @test c4[1,1,2,2] == c[1,2]
        @test c4[1,1,3,3] == c[1,3]
        @test c4[2,2,3,3] == c[2,3]
        @test c4[1,1,1,3] == c[1,5]
        @test cij(cijkl(c)) ≈ c
        c4a = cijkl(ca)
        @test c4a isa StaticArrays.MArray
        @test cij(cijkl(ca)) ≈ ca
    end
end


text:
f_xy(x,y)=x+y;
f_loga(x)=log(x+1000);


text:
# Taken from GLVisualize example "Meshcreation"

using GLVisualize, GeometryTypes, GLAbstraction, Colors

println("Visualization packages loaded.")

using ValidatedNumerics
using IntervalConstraintProgramming

println("Constraint packages loaded.")

S1 = @constraint 3 <= x^2 + y^2 + z^2 <= 4
S2 = @constraint 2 <= (x-0.5)^2 + (y-0.4)^2 + (z-0.3)^2 <= 4
S = S1 ∩ S2

X = IntervalBox(-10..10, -10..10, -10..10)

@time inner, boundary = pave(S, X, 0.2)
@show length(inner)

println("Set inversion finished")

function cube(X::IntervalBox, inner=true)
    lo = [x.lo for x in X]  # what's the good way to do this?
    hi = [x.hi for x in X]
    #color = Float32((hi[3] + 2.) / 4)
    #c = Float32(abs(hi[3]) / 2.)
    #mycolor = RGBA(1f0, c, c/2, 0.2f0)
    if inner
        mycolor = RGBA(1f0, 0f0, 0f0, 0.2f0)
    else
        mycolor = RGBA(0f0, 1f0, 0f0, 0.2f0)
    end
    return (HyperRectangle{3, Float32}(Vec3f0(lo), Vec3f0(hi - lo)), mycolor)
end

window = glscreen()

baselen = 0.4f0
dirlen = 2f0
# create an array of differently colored boxes in the direction of the 3 axes
cubes = map(cube, inner)
boundarycubes = [cube(x, false) for x in boundary]

# convert to an array of normal meshes
# note, that the constructor is a bit weird. GLNormalMesh takes a tuple of
# a geometry and a color. This means, the geometry will be converted to a GLNormalMesh
# and the color will be added afterwards, so the resulting type is a GLNormalColorMesh
meshes = map(GLNormalMesh, vcat(cubes, boundarycubes))
# merge them into one big mesh
# the resulting type is a GLNormalAttributeMesh, since we merged meshes with different
# attributes (colors). An array of the colors will be created and each vertex in the
# mesh will be asigned to one of the colors found there.
colored_mesh = merge(meshes)
view(visualize(colored_mesh), window)

renderloop(window)


text:
## set parts of BeastXML

function set_chain_length!(bx::BEASTXMLElement, cl::Int)
    mcmc_el = get_mcmc(bx)
    set_chain_length!(mcmc_el, cl)
end

function set_screen_logEvery(bx::BEASTXMLElement, sle::Int)
    mcmc_el = get_mcmc(bx)
    set_screen_logEvery(mcmc_el, sle)
end

function set_file_logEvery(bx::BEASTXMLElement, fle::Int)
    mcmc_el = get_mcmc(bx)
    set_file_logEvery(mcmc_el, fle)
end

function set_filename(bx::BEASTXMLElement, filename::AbstractString)
    mcmc_el = get_mcmc(bx)
    if isnothing(mcmc_el)
        error("No MCMC component. Cannot set log file names.")
    end
    set_filename(mcmc_el, string(filename))
    timer_el = get_timer(bx)
    if !isnothing(timer_el)
        set_filename(timer_el, "$(filename)_timer.txt")
    end
end

function set_data_dates(bx::BEASTXMLElement, dates::AbstractVector{Float64})
    data_el = get_data(bx)
    set_dates(data_el, dates)
end

function set_full_eval!(bx::BEASTXMLElement, n_eval::Int)
    mcmc = get_mcmc(bx)
    set_full_eval!(mcmc, n_eval)
end

## Model Parameters
function set_mbd_precision(bx::BEASTXMLElement, mat::AbstractArray{Float64, 2})
    mbd = get_mbd(bx)
    set_precision(mbd, mat)
end

function set_residual_precision(bx::BEASTXMLElement, mat::AbstractArray{Float64, 2})
    ext = get_extension(bx)
    set_precision(ext, mat)
end

function set_loadings(bx::BEASTXMLElement, L::AbstractArray{Float64, 2})
    loadings = get_loadings(bx)
    set_mat!(loadings, L)
end

function set_factor_precisions(bx::BEASTXMLElement, precisions::AbstractArray{Float64, 1})
    precs = get_factor_precisions(bx)
    set_value(precs, precisions)
end

function set_muliplicative_gamma_indices(bx::BEASTXMLElement, inds::Vector{Int})
    op = get_multiplicative_gamma_op(bx)
    set_indices!(op, inds)
end

function set_factor_precision_indices(bx::BEASTXMLElement, inds::Vector{Int})
    precisions = get_factor_precisions(bx)
    op = get_op_by_parameter(bx, precisions)
    set_indices!(op, inds)
end

function set_options!(bx::BEASTXMLElement, options::MCMCOptions)
    mcmc = get_mcmc(bx)
    set_options!(mcmc, options)
end