text:
"""
Simple evaluators implement a single objective function.
"""
type SimpleEvaluator <: Evaluator
  evaluator::Function
  stage::AbstractString
  threads::Int
end

"""
Provides a definition for a simple evaluator.
"""
type SimpleEvaluatorDefinition <: EvaluatorDefinition
  evaluator::Function
  stage::AbstractString
  threads::Int
end

"""
Composes a simple evaluator from its definition.
"""
function compose!(e::SimpleEvaluatorDefinition, p::Population)
  if e.stage == ""
    e.stage = p.demes[1].species.genotype.label
  end
  SimpleEvaluator(e.evaluator, e.stage, e.threads)
end

"""
TODO: Short explanation of what a simple evaluator is.
"""
simple(f::Function; stage = "", threads = 1) =
  SimpleEvaluatorDefinition(f, stage, max(1, threads))

"""
Evaluates a provided phenome according to a simple evaluation method.
"""
evaluate{T}(ev::SimpleEvaluator, sch::FitnessScheme, phenome::T) =
  ev.evaluator(sch, phenome)


text:
# This example highlights the use of a linear minimization oracle
# using an LP solver defined in MathOptInterface
# we compare the performance of the two LMOs, in- and out of place
#
# to get accurate timings it is important to run twice so that the compile time of Julia for the first run
# is not tainting the results
using FrankWolfe
using ProgressMeter
using Arpack
using Plots
using DoubleFloats
using ReverseDiff

using LinearAlgebra
using LaTeXStrings

using JuMP
const MOI = JuMP.MOI

import GLPK

n = Int(1e3)
k = 10000

xpi = rand(n);
total = sum(xpi);
const xp = xpi ./ total;

f(x) = norm(x - xp)^2
function grad!(storage, x)
    @. storage = 2 * (x - xp)
    return nothing
end

lmo_radius = 2.5
lmo = FrankWolfe.FrankWolfe.ProbabilitySimplexOracle(lmo_radius)

x00 = FrankWolfe.compute_extreme_point(lmo, zeros(n))
gradient = collect(x00)

x_lmo, v, primal, dual_gap, trajectory_lmo = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo,
    collect(copy(x00)),
    max_iteration=k,
    line_search=FrankWolfe.Shortstep(),
    L=2,
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    verbose=true,
    trajectory=true,
)

# create a MathOptInterface Optimizer and build the same linear constraints
o = GLPK.Optimizer()
x = MOI.add_variables(o, n)

# x_i ≥ 0
for xi in x
    MOI.add_constraint(o, xi, MOI.GreaterThan(0.0))
end
# ∑ x_i == 1
MOI.add_constraint(
    o,
    MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 0.0),
    MOI.EqualTo(lmo_radius),
)

lmo_moi = FrankWolfe.MathOptLMO(o)

x, v, primal, dual_gap, trajectory_moi = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo_moi,
    collect(copy(x00)),
    max_iteration=k,
    line_search=FrankWolfe.Shortstep(),
    L=2,
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    verbose=true,
    trajectory=true,
)

# formulate the LP using JuMP
m = JuMP.Model(GLPK.Optimizer)
@variable(m, y[1:n] ≥ 0)
# ∑ x_i == 1
@constraint(m, sum(y) == lmo_radius)

lmo_jump = FrankWolfe.MathOptLMO(m.moi_backend)

x, v, primal, dual_gap, trajectory_jump = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo_jump,
    collect(copy(x00)),
    max_iteration=k,
    line_search=FrankWolfe.Shortstep(),
    L=2,
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    verbose=true,
    trajectory=true,
)

x_lmo, v, primal, dual_gap, trajectory_lmo_blas = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo,
    x00,
    max_iteration=k,
    line_search=FrankWolfe.Shortstep(),
    L=2,
    print_iter=k / 10,
    emphasis=FrankWolfe.blas,
    verbose=true,
    trajectory=true,
)

x, v, primal, dual_gap, trajectory_jump_blas = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo_jump,
    x00,
    max_iteration=k,
    line_search=FrankWolfe.Shortstep(),
    L=2,
    print_iter=k / 10,
    emphasis=FrankWolfe.blas,
    verbose=true,
    trajectory=true,
)

# Defined the x-axis for the series, when plotting in terms of iterations.
iteration_list = [[x[1] + 1 for x in trajectory_lmo], [x[1] + 1 for x in trajectory_moi]]
# Defined the x-axis for the series, when plotting in terms of time.
time_list = [[x[5] for x in trajectory_lmo], [x[5] for x in trajectory_moi]]
# Defined the y-axis for the series, when plotting the primal gap.
primal_gap_list = [[x[2] for x in trajectory_lmo], [x[2] for x in trajectory_moi]]
# Defined the y-axis for the series, when plotting the dual gap.
dual_gap_list = [[x[4] for x in trajectory_lmo], [x[4] for x in trajectory_moi]]
# Defined the labels for the series using latex rendering.
label = [L"\textrm{Closed-form LMO}", L"\textrm{GLPK LMO}"]

FrankWolfe.plot_results(
    [primal_gap_list, primal_gap_list, dual_gap_list, dual_gap_list],
    [iteration_list, time_list, iteration_list, time_list],
    label,
    ["", "", L"\textrm{Iteration}", L"\textrm{Time}"],
    [L"\textrm{Primal Gap}", "", L"\textrm{Dual Gap}", ""],
    xscalelog=[:log, :identity, :log, :identity],
    yscalelog=[:log, :log, :log, :log],
    legend_position=[:bottomleft, nothing, nothing, nothing],
    filename="moi_compare.pdf",
)

ext: jl


text:
function npqreg(y, x, tau, method=IP(); m=50, h=2, xrange=nothing)
    T = eltype(x)
    if xrange === nothing
        xrange = collect(range(extrema(x)..., length=m))
    end
    z = copy(x)
    Z = hcat(fill(1, length(x)), z)
    ghat = copy(xrange)
    dghat = copy(xrange)
    for i in 1:length(xrange)
        z .= x .- xrange[i]
        Z[:, 2] .= z

        w = pdf.(Normal(), z./h)
        widx = findall(x-> x > eps(T), w)

        if length(widx) > 1
            # more than one observation with positive
            # weights, so we can estimate value and derivative

            # first we check that there is actually variation
            # in the second column
            variation = false
            z1 = Z[widx[1], 2]
            for j = 2:length(widx)
                if Z[widx[j], 2] !== z1
                    variation = true
                    break
                end
            end
            if !variation
                ghat[i] = qreg_coef(y[widx], fill(1.0, length(widx), 1), tau, method)[1]
                dghat[i] = NaN
                continue
            end

            w = w[widx]
            wy = w.*y[widx]
            wZ = w.*Z[widx, :]
            r = qreg_coef(wy, wZ, tau, method)
            ghat[i] = r[1]
            dghat[i] = r[2]
        elseif length(widx) == 0
            # no positive weights, xrange too narrow
            ghat[i] = NaN
            dghat[i] = NaN
        else
            # only one subject with positive weights, we
            # can get the value but not the derivative
            ghat[i] = y[widx[1]]
            dghat[i] = NaN
        end
    end

    xrange, ghat, dghat
end


text:
"""
Example using NLPSaUT with SNOPT7

# HS71
Polynomial objective and constraints
min x1 * x4 * (x1 + x2 + x3) + x3
st  x1 * x2 * x3 * x4 >= 25
    x1^2 + x2^2 + x3^2 + x4^2 = 40
    1 <= x1, x2, x3, x4 <= 5

Initial x = (1,5,5,1)
 (1.000..., 4.743..., 3.821..., 1.379...)
"""

using JuMP
using SNOPT7

push!(LOAD_PATH, "../src/")
using NLPSaUT

"""
User-defined fitness function

# Returns:
	- `Vector`: vector stores [f, h, g] where f is objective, h is equality, and g is inequality constraints
"""
function f_fitness(x...)
	# objective
    f = x[1]^2 - x[2]

    # equality constraints
    h = zeros(1,)
    h = x[1]^3 + x[2] - 1

    # inequality constraints
    g = zeros(2,)
    g[1] = x[1]^2 + x[2]^2 - 1
    g[2] = x[2] - 2

    fitness = vcat(f,h,g)[:]
    return fitness
end

# ------------------------------------------------------------------------- #
## Solve model
# build snopt model
model = Model(optimizer_with_attributes(SNOPT7.Optimizer,
                                        "print_level"=>5,
                                        "system_information"=>"yes"))

# problem dimensions
nx = 2                   # number of decision vectors
nh = 1                   # number of equality constraints
ng = 2                   # number of inequality constraints
lx = -10*ones(nx,)
ux =  10*ones(nx,)
x0 = [1.2, 0.9]

# get model
order = 2
diff_f = "forward"
x = build_model!(model, f_fitness, nx, nh, ng, lx, ux, x0, order, diff_f)
println(model)

# run optimizer
optimize!(model)

# print results
println(termination_status(model))
println("Decision vector: ")
println(value.(x))
println("Objective: ")
println(objective_value(model))

max_stars_repo_name: Yuricst/NLPSaUT

max_issues_repo_name: Yuricst/NLPSaUT

max_forks_repo_name: Yuricst/NLPSaUT


text:
using ExpmV, Expokit
export AQCprocessor

abstract AQCprocessor{N}

include("hamiltonian.jl")
include("aqc.jl")
include("adiaprocessing.jl")
include("demoncooling.jl")

max_stars_repo_name: Roger-luo/openQC

max_issues_repo_name: Roger-luo/openQC

max_forks_repo_name: Roger-luo/openQC


text:
using SafeTestsets
@safetestset "Benchmark Tests" begin include("debspkgtests.jl") end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/bound.jl
# Computes interval bounds of various functions.
#############################################################################

###
### AFFINE FUNCTIONS
###

function lower_interval_bound(m::Optimizer, f::AffineFunctionIneq, y::NodeBB)

    terms = f.terms
    lo_bnds = y.lower_variable_bounds
    up_bnds = y.upper_variable_bounds

    sol_branch_map = m._sol_to_branch_map
    lower_interval_bound = f.constant
    for i = 1:f.len
        coeff, indx = @inbounds terms[i]

        if m._branch_variables[indx]
            mapped_vi = @inbounds sol_branch_map[indx]
            xL = @inbounds lo_bnds[mapped_vi]
            xU = @inbounds up_bnds[mapped_vi]
        else
            xL = @inbounds m._working_problem._variable_info[indx].lower_bound
            xU = @inbounds m._working_problem._variable_info[indx].upper_bound
        end

        lower_interval_bound += (coeff > 0.0) ? coeff*xL : coeff*xU
    end

    return lower_interval_bound
end

function interval_bound(m::Optimizer, f::AffineFunctionEq, y::NodeBB)
    terms = f.terms
    lo_bnds = y.lower_variable_bounds
    up_bnds = y.upper_variable_bounds

    sol_branch_map = m._sol_to_branch_map
    lower_interval_bound = f.constant
    upper_interval_bound = f.constant
    for i = 1:f.len
        coeff, indx = @inbounds terms[i]

        if m._branch_variables[indx]
            mapped_vi = @inbounds sol_branch_map[indx]
            xL = @inbounds lo_bnds[mapped_vi]
            xU = @inbounds up_bnds[mapped_vi]
        else
            xL = @inbounds m._working_problem._variable_info[indx].lower_bound
            xU = @inbounds m._working_problem._variable_info[indx].upper_bound
        end

        if coeff > 0.0
            lower_interval_bound += coeff*xL
            upper_interval_bound += coeff*xU
        else
            lower_interval_bound += coeff*xU
            upper_interval_bound += coeff*xL
        end
    end

    return lower_interval_bound, upper_interval_bound
end

###
### QUADRATIC FUNCTIONS
###

function lower_interval_bound(m::Optimizer, f::BufferedQuadraticIneq, n::NodeBB)

    sol_branch_map = m._sol_to_branch_map
    lo_bnds = n.lower_variable_bounds
    up_bnds = n.upper_variable_bounds
    lower_interval_bound = Interval{Float64}(f.func.constant)

    for aff_term in f.func.affine_terms

        coeff = aff_term.coefficient

        vi = aff_term.variable_index.value
        if m._branch_variables[vi]
            mapped_vi = @inbounds sol_branch_map[vi]
            xL = @inbounds lo_bnds[mapped_vi]
            xU = @inbounds up_bnds[mapped_vi]
        else
            xL = @inbounds m._working_problem._variable_info[vi].lower_bound
            xU = @inbounds m._working_problem._variable_info[vi].upper_bound
        end

        lower_interval_bound += coeff > 0.0 ? coeff*xL : coeff*xU
    end

    # assumes branching on all quadratic terms, otherwise we'd need to distinguish
    # between look ups to the node bounds and lookups to the variable info in the
    # working problem.
    for quad_term in f.func.quadratic_terms

        coeff = quad_term.coefficient

        vi1 = quad_term.variable_index_1.value
        vi2 = quad_term.variable_index_2.value
        mapped_vi1 = sol_branch_map[vi1]

        xL = @inbounds lo_bnds[mapped_vi1]
        xU = @inbounds up_bnds[mapped_vi1]

        if vi1 === vi2
            if coeff > 0.0
                lower_interval_bound += (0.0 < xL) ? 0.5*coeff*xL*xL : ((xU <= 0.0) ? 0.5*coeff*xU*xU : 0.0)
            else
                lower_interval_bound += (xL < xU) ? 0.5*coeff*xU*xU : 0.5*coeff*xL*xL
            end
        else
            mapped_vi2 = sol_branch_map[vi2]

            il2b = @inbounds lo_bnds[mapped_vi2]
            iu2b = @inbounds up_bnds[mapped_vi2]
            lower_interval_bound += coeff*Interval{Float64}(xL, xU)*Interval{Float64}(il2b, iu2b)

        end
    end

    return lower_interval_bound.lo
end

function interval_bound(m::Optimizer, f::BufferedQuadraticEq, n::NodeBB)

    sol_branch_map = m._sol_to_branch_map
    lo_bnds = n.lower_variable_bounds
    up_bnds = n.upper_variable_bounds
    val_intv = Interval(f.func.constant)

    for aff_term in f.func.affine_terms

        coeff = aff_term.coefficient

        vi = aff_term.variable_index.value
        if m._branch_variables[vi]
            mapped_vi = @inbounds sol_branch_map[vi]
            xL = @inbounds lo_bnds[mapped_vi]
            xU = @inbounds up_bnds[mapped_vi]
        else
            xL = @inbounds m._working_problem._variable_info[vi].lower_bound
            xU = @inbounds m._working_problem._variable_info[vi].upper_bound
        end

        val_intv += coeff*Interval(xL, xU)
    end

    # assumes branching on all quadratic terms, otherwise we'd need to distinguish
    # between look ups to the node bounds and lookups to the variable info in the
    # working problem.
    for quad_term in f.func.quadratic_terms

        coeff = quad_term.coefficient

        vi1 = quad_term.variable_index_1.value
        vi2 = quad_term.variable_index_2.value

        mapped_vi1 = @inbounds sol_branch_map[vi1]

        xL = @inbounds lo_bnds[mapped_vi1]
        xU = @inbounds up_bnds[mapped_vi1]

        if vi1 === vi2
            val_intv += 0.5*coeff*pow(Interval(xL, xU), 2)
        else
            mapped_vi2 = @inbounds sol_branch_map[vi2]
            @inbounds il2b = lo_bnds[mapped_vi2]
            @inbounds iu2b = up_bnds[mapped_vi2]
            val_intv += coeff*Interval(xL, xU)*Interval(il2b, iu2b)
        end
    end

    return val_intv.lo, val_intv.hi
end

###
### SECOND-ORDER CONE
###
function lower_interval_bound(m::Optimizer, d::BufferedSOC, n::NodeBB)

    sol_branch_map = m._sol_to_branch_map
    lo_bnds = n.lower_variable_bounds
    up_bnds = n.upper_variable_bounds
    vec_of_vi = d.variables.variables

    norm_bound = Interval(0.0)
    for i = 2:length(vec_of_vi)
        mapped_vi = @inbounds sol_branch_map[vec_of_vi[i].value]
        x = Interval{Float64}(lo_bnds[mapped_vi], up_bnds[mapped_vi])
        norm_bound += pow(x, 2)
    end
    norm_bound = sqrt(norm_bound)

    mapped_vi = @inbounds sol_branch_map[vec_of_vi[1].value]
    lower_bound = norm_bound.lo -(@inbounds up_bnds[mapped_vi])

    return lower_bound
end

###
### NONLINEAR FUNCTIONS
###

function lower_interval_bound(m::Optimizer, d::BufferedNonlinearFunction{V}, n::NodeBB) where V
    if !d.has_value
        forward_pass!(m._working_problem._relaxed_evaluator, d)
    end

    expr = d.expr
    if expr.isnumber[1]
        lower_value = expr.numberstorage[1]
    else
        lower_value = expr.setstorage[1].Intv.lo
    end

    return lower_value
end

function interval_bound(m::Optimizer, d::BufferedNonlinearFunction{V}, n::NodeBB) where V
    if !d.has_value
        forward_pass!(d.evaluator, d)
    end

    expr = d.expr
    if expr.isnumber[1]
        interval_value = Interval(expr.numberstorage[1])
    else
        interval_value = expr.setstorage[1].Intv
    end

    return interval_value.lo, interval_value.hi
end


text:
"""
spin24.jl - multiple sampling robustness for the δ
"""

WDIR = joinpath(@__DIR__, "../../")
include(joinpath(WDIR, "src", "spin", "spin.jl"))

using Altro
using Distributions
using HDF5
using ForwardDiff
using LinearAlgebra
using Random
using RobotDynamics
using StaticArrays
using Zygote
using TrajectoryOptimization
const RD = RobotDynamics
const TO = TrajectoryOptimization

# paths
const EXPERIMENT_META = "spin"
const EXPERIMENT_NAME = "spin24"
const SAVE_PATH = joinpath(WDIR, "out", EXPERIMENT_META, EXPERIMENT_NAME)

# problem
const CONTROL_COUNT = 1
const ASTATE_SIZE = HDIM_ISO
const ACONTROL_SIZE = CONTROL_COUNT
const INITIAL_STATE1 = [1., 0, 0, 0]
# state indices
const STATE1_IDX = 1:HDIM_ISO
# control indices
const CONTROL_IDX = 1:CONTROL_COUNT

# model
mutable struct Model <: RD.AbstractModel
end
@inline RD.state_dim(model::Model) = ASTATE_SIZE
@inline RD.control_dim(model::Model) = ACONTROL_SIZE

# dynamics
function RD.discrete_dynamics(::Type{RK3}, model::Model, astate::StaticVector{ASTATE_SIZE},
                              acontrol::StaticVector{ACONTROL_SIZE}, time::Real, dt::Real)
    negi_hc = acontrol[CONTROL_IDX[1]] * NEGI_H1_ISO
    state1 = exp(dt * (FQ_NEGI_H0_ISO + negi_hc)) * astate
    return state1
end

# cost function
struct Cost{N,M,T} <: TO.CostFunction
    hess_astate::Symmetric{T, SMatrix{N,N,T}}
    target_state1::SVector{HDIM_ISO, T}
    R::Diagonal{T, SVector{M,T}}
    q::T
end

function Cost(target_state1::SVector{HDIM_ISO, T}, q::T, r::T) where {T}
    N = ASTATE_SIZE
    M = ACONTROL_SIZE
    # For reasons unkown to the author, this cost function works properly
    # when a -1 is thrown infront of the hessian.
    hess_astate = -1 * q * hessian_gate_error_iso2(target_state1)
    hess_astate = Symmetric(SMatrix{N, N}(hess_astate))
    R = Diagonal(SVector{M}(r))
    return Cost{N,M,T}(hess_astate, target_state1, R, q)
end
@inline TO.state_dim(cost::Cost{N,M,T}) where {N,M,T} = N
@inline TO.control_dim(cost::Cost{N,M,T}) where {N,M,T} = M
@inline Base.copy(cost::Cost{N,M,T}) where {N,M,T} = Cost{N,M,T}(
    cost.hess_astate, cost.target_state1, cost.R, cost.q
)

@inline TO.stage_cost(cost::Cost{N,M,T}, astate::SVector{N}) where {N,M,T} = (
    cost.q * gate_error_iso2(astate, cost.target_state1)
)

@inline TO.stage_cost(cost::Cost{N,M,T}, astate::SVector{N}, acontrol::SVector{M}) where {N,M,T} = (
    TO.stage_cost(cost, astate) + 0.5 * acontrol' * cost.R * acontrol
)

function TO.gradient!(E::TO.QuadraticCostFunction, cost::Cost{N,M,T}, astate::SVector{N,T}) where {N,M,T}
    E.q = cost.q * jacobian_gate_error_iso2(astate, cost.target_state1)
    return false
end

function TO.gradient!(E::TO.QuadraticCostFunction, cost::Cost{N,M,T}, astate::SVector{N,T},
                      acontrol::SVector{M,T}) where {N,M,T}
    TO.gradient!(E, cost, astate)
    E.r = cost.R * acontrol
    E.c = 0
    return false
end

function TO.hessian!(E::TO.QuadraticCostFunction, cost::Cost{N,M,T}, astate::SVector{N,T}) where {N,M,T}
    E.Q = cost.hess_astate
    return true
end

function TO.hessian!(E::TO.QuadraticCostFunction, cost::Cost{N,M,T}, astate::SVector{N,T},
                     acontrol::SVector{M,T}) where {N,M,T}
    TO.hessian!(E, cost, astate)
    E.R = cost.R
    E.H .= 0
    return true
end

# main
function run_traj(;gate_type=xpiby2, evolution_time=56.8, solver_type=altro,
                  sqrtbp=false, integrator_type=rk3, q=1e0, r=1e0,
                  dt_inv=Int64(1e1), smoke_test=false, constraint_tol=1e-8, al_tol=1e-4,
                  pn_steps=2, max_penalty=1e8, verbose=true, save=true,
                  max_iterations=Int64(2e5), line_search_lower_bound=1e-8,
                  line_search_upper_bound=10, dJ_counter_limit=Int(1e3))
    model = Model()
    n = RD.state_dim(model)
    m = RD.control_dim(model)
    t0 = 0.

    # initial state
    x0 = SVector{n}([
        INITIAL_STATE1;
    ])

    # target state
    if gate_type == xpiby2
        target_state1 = XPIBY2_ISO_1
        target_state2 = XPIBY2_ISO_2
    elseif gate_type == ypiby2
        target_state1 = YPIBY2_ISO_1
        target_state2 = YPIBY2_ISO_2
    elseif gate_type == zpiby2
        target_state1 = ZPIBY2_ISO_1
        target_state2 = ZPIBY2_ISO_2
    end
    xf = SVector{n}([
        target_state1;
    ])

    # initial trajectory
    dt = dt_inv^(-1)
    N = Int(floor(evolution_time * dt_inv)) + 1
    U0 = [SVector{m}([
        fill(1e-2, CONTROL_COUNT);
    ]) for k = 1:N-1]
    X0 = [SVector{n}([
        fill(NaN, n);
    ]) for k = 1:N]
    Z = Traj(X0, U0, dt * ones(N))

    # cost function
    ck = Cost(target_state1, q, r)
    cf = Cost(target_state1, q * N, r)
    objective = Objective(ck, cf, N)

    # constraints
    constraints = ConstraintList(n, m, N)

    # solve
    prob = Problem{IT_RDI[integrator_type]}(model, objective, constraints, x0, xf, Z, N, t0, evolution_time)
    solver = ALTROSolver(prob)
    verbose_pn = verbose ? true : false
    verbose_ = verbose ? 2 : 0
    projected_newton = solver_type == altro ? true : false
    constraint_tolerance = solver_type == altro ? constraint_tol : al_tol
    iterations_inner = smoke_test ? 1 : 300
    iterations_outer = smoke_test ? 1 : 30
    n_steps = smoke_test ? 1 : pn_steps
    set_options!(
        solver, square_root=sqrtbp, constraint_tolerance=constraint_tolerance,
        projected_newton_tolerance=al_tol, n_steps=n_steps,
        penalty_max=max_penalty, verbose_pn=verbose_pn, verbose=verbose_,
        projected_newton=projected_newton, iterations_inner=iterations_inner,
        iterations_outer=iterations_outer, iterations=max_iterations,
        line_search_lower_bound=line_search_lower_bound,
        line_search_upper_bound=line_search_upper_bound,
        dJ_counter_limit=dJ_counter_limit,
    )
    Altro.solve!(solver)

    # post-process
    acontrols_raw = TO.controls(solver)
    acontrols_arr = permutedims(reduce(hcat, map(Array, acontrols_raw)), [2, 1])
    astates_raw = TO.states(solver)
    astates_arr = permutedims(reduce(hcat, map(Array, astates_raw)), [2, 1])
    cmax = TO.max_violation(solver)
    cmax_info = TO.findmax_violation(TO.get_constraints(solver))
    iterations_ = Altro.iterations(solver)

    result = Dict(
        "acontrols" => acontrols_arr,
        "evolution_time" => evolution_time,
        "astates" => astates_arr,
        "q" => q,
        "r" => r,
        "cmax" => cmax,
        "cmax_info" => cmax_info,
        "dt" => dt,
        "solver_type" => Integer(solver_type),
        "sqrtbp" => Integer(sqrtbp),
        "max_penalty" => max_penalty,
        "constraint_tol" => constraint_tol,
        "al_tol" => al_tol,
        "integrator_type" => Integer(integrator_type),
        "gate_type" => Integer(gate_type),
        "save_type" => Integer(jl),
        "iterations" => iterations_,
        "max_iterations" => max_iterations,
    )

    return result
end

max_stars_repo_name: SchusterLab/rbqoc

max_issues_repo_name: SchusterLab/rbqoc

max_forks_repo_name: SchusterLab/rbqoc


text:
TimeCost = []
D = Distributions.Normal(0.0, 1.0)
for x in 1:10000
  # sampling, and get the time cost (gc time included), then store
  push!(TimeCost, @timed( NovDist.LatinHyperCube(D, 100,100) )[2] )
end
# averaging
StatsBase.mean(TimeCost)

PyPlot.plt[:hist]( tmp, 40, density = true );
PyPlot.grid(true); PyPlot.ylabel("Density"); PyPlot.xlabel("x");
PyPlot.title("PDF of 10000 samples from N(0,1)");
PyPlot.tight_layout()

TimeCost = []
Sigma = [0.5 -0.25; -0.25 1.5];
Mu = [ 0.0, 1.0 ];
D = Distributions.MvNormal(Mu, Sigma)
for x in 1:10000
  push!(TimeCost, @timed( NovDist.LatinHyperCube(D, 100,100) )[2] )
end
StatsBase.mean(TimeCost)
StatsBase.std(TimeCost)

#

max_stars_repo_name: Clpr/NumToolLibs

max_issues_repo_name: Clpr/NumToolLibs

max_forks_repo_name: Clpr/NumToolLibs


text:
module DBdatatype # Note: correspondence here is complicated by platform issues.
# See http://julia.readthedocs.org/en/release-0.3/manual/calling-c-and-fortran-code/
using Compat

const DB_INT = 16
const DB_SHORT = 17
const DB_LONG = 18
const DB_FLOAT = 19
const DB_DOUBLE = 20
const DB_CHAR = 21
const DB_LONG_LONG = 22
const DB_NOTYPE = 25           #=unknown type =#

const JuliaSilotypemap= @compat Dict(Int64 => DB_LONG_LONG, Int32 => DB_INT, Int16 =>DB_SHORT, Int8=>DB_CHAR, Float32 => DB_FLOAT, Float64 => DB_DOUBLE)
# const CJuliaSilotypemap=[Clonglong => DB_LONG_LONG, Cint => DB_INT, Cshort =>DB_SHORT, Cchar=>DB_CHAR, Cfloat => DB_FLOAT, Cdouble => DB_DOUBLE]
const SiloJuliatypemap = @compat Dict(DB_LONG_LONG => Int64, DB_INT=> Int32, DB_SHORT=>Int16, DB_CHAR=>Int8, DB_FLOAT=>Float32, DB_DOUBLE=>Float64 )
# const SiloCJuliatypemap=[DB_LONG_LONG => Clonglong, DB_INT=> Cint, DB_SHORT=>Cshort , DB_CHAR=>Cchar, Cfloat => DB_FLOAT, DB_DOUBLE=>Cdouble]

end #module


text:
using Pkg
Pkg.activate(".")

using RadMod2D

include("./models2D.jl")
include("./plot2D.jl")

# backend
# using GLMakie
using CairoMakie
CairoMakie.activate!(type = "svg")
filetype = ".svg"

# parameter
r1x = 0.8
r1y = 0.8
r2x = 1.4
r2y = 1.4
elemsize = 0.05
n = 15

# model
m = model_rect_in_rect(r1x, r1y, r2x, r2y, elemsize)

# 2D plot of model
fig1 = Figure(resolution = (350, 350), font = "Arial", fontsize = 12)
ax1 = fig1[1,1] = Axis(fig1)
ax1.xlabel = "X in m"
ax1.ylabel = "Y in m"
ax1.aspect = DataAspect()
plot_model(fig1, ax1, m, shownvec = true, shownodes = false, showcom = false, showleg = true, legpos = "right", linewidth = 1.5)

# view factors
vfmat = zeros(Float64, m.nelem, m.nelem)
existing_vf!(m, vfmat)
dx, dy = create_tiles(m, n)
t_occ = check_tile_occupation(m, dx, dy, n)
blocking_vf_with_tiles!(m, vfmat, dx, dy, n, t_occ)
calculating_vf!(m, vfmat, normit = true)

# net radiation method
epsilon = zeros(m.nelem,1)
set_bc_part!(m, epsilon, 1, 0.3)
set_bc_part!(m, epsilon, 2, 0.6)
temp = zeros(m.nelem,1)
set_bc_part!(m, temp, 1:2, 300)
temp[79:106,1] .= 600 # bc for specific elements
@time Q, G = tempsolver(m, vfmat, temp, epsilon)
area = [m.elem[i].area for i = m.elem2par[1].first:m.elem2par[end].last]
q = Q[:] ./ area[:]

# 2D plot of heat flux density
fig2 = Figure(resolution = (350, 350), font = "Arial", fontsize = 12)
ax2 = fig2[1,1] = Axis(fig2)
ax2.xlabel = "X in m"
ax2.ylabel = "Y in m"
ax2.aspect = DataAspect()
cbar = plot_model_with_value(fig2, ax2, m, q, "Heat flux density in W/m", cmap = :jet, showcbar = true, cbarpos = "right", linewidth = 1.5)
cbar.height = Relative(2/3)

# display(fig1)
# display(fig2)
save("fig_rect_in_rect_mesh2D"*filetype, fig1)
save("fig_rect_in_rect_q2D"*filetype, fig2)

# include("./test/run_example2.jl")


text:
export ELFStrTab

"""
    ELFStrTab

ELF string table type, containing information about a `Section` that contains
a string table.  Note that an ELF file may contain multiple string tables, use
the different constructors to reference the different string tables.
"""
struct ELFStrTab{H <: ELFHandle} <: StrTab{H}
    section_ref::SectionRef{H}

    # Do a little bit of section checking
    function ELFStrTab{H}(section::ELFSectionRef{H}) where {H <: ELFHandle}
        # If this isn't actually a string table section, warn the user!
        if section_type(section) != SHT_STRTAB
            section_idx = section_number(section)
            msg = strip("""
            Loading Section $(section_idx) as a StrTab, despite it having section
            type $(section_type_string(section)), not SHT_STRTAB.
            """)
            warn(replace(msg, "\n", " "), key=(handle(section), section_idx))
        end

        return new(section)
    end
end

"""
    StrTab(section::SectionRef{ELFHandle})

Special short-cut constructor to build an `StrTab` directly from an ELF
`Section`, without attempting to automatically discover the correct section for
a string table.
"""
StrTab(section::ELFSectionRef{H}) where {H <: ELFHandle} = ELFStrTab{H}(section)

"""
    StrTab(oh::ELFHandle)

Given an ELF object, construct an `ELFStrTab` that refers to the section
header string table.
"""
function StrTab(oh::H) where {H <: ELFHandle}
    section_index = section_header_strtab_index(oh)
    strtab_section = Sections(oh)[section_index]
    return ELFStrTab{H}(strtab_section)
end

handle(strtab::ELFStrTab) = handle(strtab.section_ref)

"""
    section_header_strtab_index(oh::ELFHandle)

In an ELF object file, the section header string table is located within a
designated ELF `Section`, this method returns the index of that section.  This
is used by `StrTab(::ELFHandle)` to locate the correct `SectionRef` to
construct an `StrTab` object from.  The returned value is an index (e.g. it is
1-based, not zero-based as the value within the ELF object itself).
"""
function section_header_strtab_index(oh::ELFHandle)
    return header(oh).e_shstrndx + 1
end

function strtab_lookup(strtab::ELFStrTab, index)
    seek(strtab.section_ref, index)
    return strip(readuntil(handle(strtab), '\0'), '\0')
end


text:
# Register functions and handle literals
"""
$(SIGNATURES)

Registers a function call as a primitive for the `Operation` graph of the
ModelingToolkit IR. Example:

```julia
@register f(x,y)
```

registers `f` as a possible two-argument function.

You may also want to tell ModelingToolkit the derivative of the registered
function. Here is an example to do it

```julia
julia> using ModelingToolkit

julia> foo(x, y) = sin(x) * cos(y)
foo (generic function with 1 method)

julia> @parameters t; @variables x(t) y(t) z(t); @derivatives D'~t;

julia> @register foo(x, y)
foo (generic function with 4 methods)

julia> foo(x, y)
foo(x(t), y(t))

julia> ModelingToolkit.derivative(::typeof(foo), (x, y), ::Val{1}) = cos(x) * cos(y) # derivative w.r.t. the first argument

julia> ModelingToolkit.derivative(::typeof(foo), (x, y), ::Val{2}) = -sin(x) * sin(y) # derivative w.r.t. the second argument

julia> isequal(expand_derivatives(D(foo(x, y))), expand_derivatives(D(sin(x) * cos(y))))
true
```
"""
macro register(sig)
    splitsig = splitdef(:($sig = nothing))
    name = splitsig[:name]

    # Extract the module and function name from the signature
    if name isa Symbol
        mod = __module__  # Calling module
        funcname = name
    else
        mod = name.args[1]
        funcname = name.args[2].value
    end

    args = splitsig[:args]
    typargs = typed_args(args)
    defs = :()
    for typarg in typargs
        splitsig[:args] = typarg
        if mod == (@__MODULE__)  # If the calling module is ModelingToolkit itself...
            splitsig[:body] = :(Operation($name, Expression[$(args...)]))
        else
            # Register the function's associated model so we can inject it in later.
            splitsig[:body] = quote
                get!(ModelingToolkit.registered_external_functions, Symbol($("$funcname")), $mod)
                Operation($name, Expression[$(args...)])
            end
        end
        defs = :($defs; $(combinedef(splitsig)))
    end
    esc(defs)
end
# Create all valid combinations of Expression,Number for function signature
function typed_args(args)
    cases = Iterators.product(Tuple((:Expression, :Number) for i = 1:length(args))...)
    typargs = Vector{Any}[]
    for case in cases
        typarg = [Expr(:(::), args[i], case[i]) for i = 1:length(args)]
        push!(typargs, typarg)
    end
    typargs[1:end-1] # Last case is all Numbers but that should be specialized
end

# Binary & unary operators and functions
import DiffRules
for (M, f, arity) in DiffRules.diffrules()
    fun = :($M.$f)
    sig = arity == 1 ? :($fun(x)) :
          arity == 2 && :($fun(x,y))
    @eval @register $sig
end

for fun ∈ [:!]
    basefun = Expr(:., Base, QuoteNode(fun))
    sig = :($basefun(x))
    @eval @register $sig
end

for fun ∈ [:<, :>, :(==), :&, :|, :div]
    basefun = Expr(:., Base, QuoteNode(fun))
    sig = :($basefun(x,y))
    @eval @register $sig
end

# ifelse
#@register Base.ifelse(cond,t,f)

# special cases
Base.:^(x::Expression,y::T) where T <: Integer = Operation(Base.:^, Expression[x, y])
Base.:^(x::Expression,y::T) where T <: Rational = Operation(Base.:^, Expression[x, y])

@register Base.conj(x)
@register Base.getindex(x,i)
@register Base.binomial(n,k)
@register Base.copysign(x,y)

Base.getindex(x::Operation,i::Int64) = Operation(getindex,[x,i])
Base.one(::Operation) = 1

# Ensure that Operations that get @registered from outside the ModelingToolkit
# module can work without having to bring in the associated function into the
# ModelingToolkit namespace. We basically store information about functions
# registered at runtime in a ModelingToolkit variable,
# `registered_external_functions`. It's not pretty, but we are limited by the
# way GeneralizedGenerated builds a function (adding "ModelingToolkit" to every
# function call).
# ---
const registered_external_functions = Dict{Symbol,Module}()
function inject_registered_module_functions(expr)
    MacroTools.postwalk(expr) do x
        # Find all function calls in the expression and extract the function
        # name and calling module.
        MacroTools.@capture(x, f_module_.f_name_(xs__))
        if isnothing(f_module)
            MacroTools.@capture(x, f_name_(xs__))
        end

        if !isnothing(f_name)
            # Set the calling module to the module that registered it.
            mod = get(registered_external_functions, f_name, f_module)
            if !isnothing(mod)
                x.args[1] = :(getproperty($mod, $(Meta.quot(f_name))))
            end
        end

        return x
    end
end


text:
@testset "pawns" begin
    bb = Bobby.set_board()
    pl = Array{Bobby.Move,1}()
    pvl = Bobby.get_pawns_list(pl, bb)
    @test length(pvl) == 16

    pl = Array{Bobby.Move,1}()
    pvl = Bobby.get_pawns_list(pl, bb, "black")
    @test length(pvl) == 16

    @test test_fen("k7/8/8/8/8/1p6/P7/7K w - - 0 1", 1, [6])
    @test test_fen("k7/8/8/8/8/8/P7/7K w - - 0 1", 1, [5])
    @test test_fen("k7/8/8/8/1p6/8/P7/7K w - - 0 1", 1, [5])
    @test test_fen("k7/8/8/8/8/1p6/P7/7K b - - 0 1", 1, [5])
    @test test_fen("k7/8/8/8/Pp6/8/8/7K b - a3 0 1", 1, [5])
    @test test_fen("k7/8/8/8/P7/1p6/8/7K b - - 0 1", 1, [4])
    @test test_fen("k7/8/8/8/P7/8/1p6/7K b - - 0 1", 1, [7])
    @test test_fen("k7/8/8/8/P7/8/7p/6NK b - - 0 1", 1, [7])
    @test test_fen("k7/8/8/8/8/p7/P7/7K w - - 0 1", 1, [3])
    @test test_fen("k7/8/8/8/8/2p5/3P4/7K w - - 0 1", 1, [6])
    @test test_fen("k7/8/8/8/8/2p1p3/3P4/7K w - - 0 1", 1, [7])
    @test test_fen("k7/8/8/8/2p1p3/3P1P2/8/7K b - d3 0 1", 1, [8])
    @test test_fen("k7/3p4/8/8/8/8/8/7K b - - 0 1", 1, [5])
    @test test_fen("k7/8/8/2P5/8/8/8/7K w - - 0 1", 1, [4])
    @test test_fen("k7/8/8/2Pp4/8/8/8/7K w - - 0 1", 1, [4])
    @test test_fen("k7/8/8/2Pp4/8/8/8/7K w - d6 0 1", 1, [5])
    @test test_fen("k7/3pP3/8/8/8/8/8/7K w - - 0 1", 1, [7])
    @test test_fen("k7/8/1P6/8/8/8/8/7K b - - 0 1", 1, [2])
    @test test_fen("k7/1P6/8/8/8/8/8/7K b - - 0 1", 1, [3])
    @test test_fen("k7/8/8/8/8/2P5/1P6/7K w - - 0 1", 1, [6])
end


text:
include("../utils.jl")

@testset "Cholesky Space" begin
    M = Manifolds.CholeskySpace(3)
    @test repr(M) == "CholeskySpace(3)"

    types = [Matrix{Float64}]
    TEST_FLOAT32 && push!(types, Matrix{Float32})
    TEST_STATIC_SIZED && push!(types, MMatrix{3,3,Float64,9})

    for T in types
        A(α) = [1.0 0.0 0.0; 0.0 cos(α) sin(α); 0.0 -sin(α) cos(α)]
        ptsF = [#
            cholesky([1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1]).L,
            cholesky([2.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 1]).L,
            cholesky(
                A(π / 6) * [1.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 1] * transpose(A(π / 6)),
            ).L,
        ]
        pts = [convert(T, a) for a in ptsF]
        test_manifold(
            M,
            pts;
            test_injectivity_radius=false,
            test_default_vector_transport=true,
            vector_transport_methods=[
                ParallelTransport(),
                SchildsLadderTransport(),
                PoleLadderTransport(),
            ],
            test_vee_hat=false,
            exp_log_atol_multiplier=8.0,
            test_inplace=true,
        )
    end
    @testset "Test Error cases in is_point and is_vector" begin
        pt1f = zeros(2, 3) # wrong size
        pt2f = [1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 1.0] # nonpos diag
        pt3f = [2.0 0.0 1.0; 0.0 1.0 0.0; 0.0 0.0 4.0] # no lower and nonsym
        pt4 = [2.0 0.0 0.0; 1.0 2.0 0.0; 0.0 0.0 4.0]
        @test !is_point(M, pt1f)
        @test_throws DomainError is_point(M, pt1f, true)
        @test !is_point(M, pt2f)
        @test_throws DomainError is_point(M, pt2f, true)
        @test !is_point(M, pt3f)
        @test_throws DomainError is_point(M, pt3f, true)
        @test is_point(M, pt4)
        @test !is_vector(M, pt3f, pt1f)
        @test_throws DomainError is_vector(M, pt3f, pt1f, true)
        @test !is_vector(M, pt4, pt1f)
        @test_throws DomainError is_vector(M, pt4, pt1f, true)
        @test !is_vector(M, pt4, pt3f)
        @test_throws DomainError is_vector(M, pt4, pt3f, true)
        @test is_vector(M, pt4, pt2f)
    end
end


text:
module Laplace

export postwid, _PWcoeffs
export LT_hyperbola
export LT_hyper_fixed
#export talbot_improved

include("postwidder.jl")
include("contours.jl")

end


text:
#############################################################################
# ECOS.jl
# Wrapper around the ECOS solver https://github.com/ifa-ethz/ecos
# See http://github.com/JuliaOpt/ECOS.jl
#############################################################################
# ECOSSolverInterface.jl
# MathProgBase.jl interface for the ECOS.jl solver wrapper
#############################################################################

import MathProgBase
const MPB = MathProgBase

#############################################################################
# Define the MPB Solver and Model objects
export ECOSSolver
struct ECOSSolver <: MPB.AbstractMathProgSolver
    options
end
ECOSSolver(;kwargs...) = ECOSSolver(kwargs)

mutable struct ECOSMathProgModel <: MPB.AbstractConicModel
    nvar::Int                           # Number of variables
    nconstr::Int                        # Number of rows in the MPB-form A matrix
    nineq::Int                          # Number of inequalities Gx <=_K h
    neq::Int                            # Number of equalities Ax = b
    npos::Int                           # Number of positive orthant cones
    ncones::Int                         # Number of SO cones
    conedims::Vector{Int}               # Dimension of each SO cone
    nexp_cones::Int                     # Number of exponential cones
    G::ECOSMatrix                       # The G matrix (inequalties)
    A::ECOSMatrix                       # The A matrix (equalities)
    c::Vector{Float64}                  # The objective coeffs (always min)
    orig_sense::Symbol                  # Original objective sense
    h::Vector{Float64}                  # RHS for inequality
    b::Vector{Float64}                  # RHS for equality
    # Post-solve
    solve_stat::Symbol
    solve_time::Float64
    obj_val::Float64
    primal_sol::Vector{Float64}
    dual_sol_eq::Vector{Float64}
    dual_sol_ineq::Vector{Float64}
    # Maps x∈K to ECOS duals
    # .._ind maps a column to a row index
    # .._type  is the original cone type of each column
    col_map_ind::Vector{Int}
    col_map_type::Vector{Symbol}
    # Maps b-Ax∈K to ECOS duals
    # .._ind maps a row to an index
    # .._type  is the original cone type of each row
    row_map_ind::Vector{Int}
    row_map_type::Vector{Symbol}
    # To reorder solution to MPB form
    fwd_map::Vector{Int}
    options
end
ECOSMathProgModel(;kwargs...) = ECOSMathProgModel(0,0,0,0,0,0,
                                        Int[],0,
                                        ECOSMatrix(spzeros(0,0)),
                                        ECOSMatrix(spzeros(0,0)),
                                        Float64[], :Min,
                                        Float64[], Float64[],
                                        :NotSolved, 0.0, 0.0,
                                        Float64[],
                                        Float64[], Float64[],
                                        Int[], Symbol[],
                                        Int[], Symbol[], Int[],
                                        kwargs)

#############################################################################
# Begin implementation of the MPB low-level interface
# Implements
# - model
# - optimize!
# - status
# - numvar
# - numconstr
# http://mathprogbasejl.readthedocs.org/en/latest/lowlevel.html

MPB.ConicModel(s::ECOSSolver) = ECOSMathProgModel(;s.options...)
MPB.LinearQuadraticModel(s::ECOSSolver) = MPB.ConicToLPQPBridge(MPB.ConicModel(s))

function MPB.optimize!(m::ECOSMathProgModel)
    ecos_prob_ptr = setup(
        m.nvar, m.nineq, m.neq,
        m.npos, m.ncones, m.conedims, m.nexp_cones,
        m.G, m.A,
        m.c,
        m.h, m.b; m.options...)
    # Note: ECOS modifies problem data in setup() and restores it on cleanup()

    t = time()
    flag = solve(ecos_prob_ptr)
    m.solve_time = time() - t
    if flag == ECOS_OPTIMAL
        m.solve_stat = :Optimal
    elseif flag == ECOS_PINF
        m.solve_stat = :Infeasible
    elseif flag == ECOS_DINF  # Dual infeasible = primal unbounded, probably
        m.solve_stat = :Unbounded
    elseif flag == ECOS_MAXIT
        m.solve_stat = :UserLimit
    elseif flag == ECOS_OPTIMAL + ECOS_INACC_OFFSET
        m.solve_stat = :Suboptimal
    else
        m.solve_stat = :Error
    end
    # Extract solution
    ecos_prob = unsafe_wrap(Array, ecos_prob_ptr, 1)[1]
    m.primal_sol = unsafe_wrap(Array,ecos_prob.x, m.nvar)[:]
    m.dual_sol_eq   = unsafe_wrap(Array,ecos_prob.y, m.neq)[:]
    m.dual_sol_ineq = unsafe_wrap(Array,ecos_prob.z, m.nineq)[:]
    cleanup(ecos_prob_ptr, 0)
    # Compute this value after cleanup with restored c.
    m.obj_val = dot(m.c, m.primal_sol) * (m.orig_sense == :Max ? -1 : +1)
end

MPB.status(m::ECOSMathProgModel) = m.solve_stat
MPB.getobjval(m::ECOSMathProgModel) = m.obj_val
MPB.getsolution(m::ECOSMathProgModel) = m.primal_sol[m.fwd_map]

MPB.numvar(m::ECOSMathProgModel) = m.nvar
MPB.numconstr(m::ECOSMathProgModel) = m.nineq + m.neq

#############################################################################
# Begin implementation of the MPB conic interface
# Implements
# - supportedcones
# - loadconicproblem!
# - getconicdual
# http://mathprogbasejl.readthedocs.org/en/latest/conic.html

MPB.supportedcones(m::ECOSSolver) = [:Free,:Zero,:NonNeg,:NonPos,:SOC,:ExpPrimal]

function MPB.loadproblem!(m::ECOSMathProgModel, c, A, b, constr_cones, var_cones)
    if size(A,2) == 0
        Base.warn("Input matrix has no columns. ECOS is known to crash in this corner case.")
    end
    m.nconstr = size(A,1)
    # If A is sparse, we should use an appropriate "zeros"
    zeromat = isa(A,SparseMatrixCSC) ? spzeros : zeros

    # We don't support SOCRotated, SDP, or ExpDual
    bad_cones = [:SOCRotated, :SDP, :ExpDual]
    for cone_vars in constr_cones
        cone_vars[1] in bad_cones && error("Cone type $(cone_vars[1]) not supported")
    end
    for cone_vars in var_cones
        cone_vars[1] in bad_cones && error("Cone type $(cone_vars[1]) not supported")
    end

    # MathProgBase form             ECOS form
    # min  c'x                      min  c'x
    #  st b-Ax ∈ K_1                 st   Ax = b
    #        x ∈ K_2                    h-Gx ∈ K
    #
    # Mapping:
    # * For the constaints (K_1)
    #   * If :Zero cone, then treat as equality constraint in ECOS form
    #   * Otherewise trivially maps to h-Gx in ECOS form
    # * For the variables (K_2)
    #   * If :Free, do nothing
    #   * If :Zero, put in as equality constraint
    #   * If rest, stick in h-Gx
    #
    # Approach:
    # 0. Figure out a mapping between MPB and ECOS
    # 1. Map non-SOC/Exp variables to ECOS form
    # 2. Map non-SOC/Exp constraints to ECOS form
    # 3. Map SOC variables and constraints to ECOS form
    # 4. Map Exp variables and constraints to ECOS form

    # Allocate space for the ECOS variables
    num_vars = length(c)
    fwd_map = Array{Int}(undef, num_vars)  # Will be used for SOCs
    rev_map = Array{Int}(undef, num_vars)  # Need to restore sol. vec.
    idxcone = Array{Symbol}(undef, num_vars)  # We'll use this for non-SOC/Exp
    m.col_map_type = Array{Symbol}(undef, num_vars) # In MPB indices

    # Now build the mapping between MPB variables and ECOS variables
    pos = 1
    for (cone, idxs) in var_cones
        for i in idxs
            fwd_map[i]   = pos   # fwd_map = MPB idx -> ECOS idx
            rev_map[pos] = i     # rev_map = ECOS idx -> MPB idx
            idxcone[pos] = cone
            m.col_map_type[i] = cone
            pos += 1
        end
    end

    # Rearrange data into the internal ordering, make copy
    ecos_c = c[rev_map]
    ecos_A = zeromat(0,num_vars)
    ecos_b = Float64[]

    # Mapping for duals
    m.col_map_ind = zeros(Int,num_vars)
    m.row_map_ind = zeros(Int, length(b))
    m.row_map_type  = Array{Symbol}(undef, length(b))
    function update_row_map(cone_type, cur_ind)
        for (cone,idxs) in constr_cones
            if cone == cone_type
                for idx in idxs
                    m.row_map_ind[idx] = cur_ind
                    m.row_map_type[idx] = cone_type
                    cur_ind += 1
                end
            end
        end
        cur_ind
    end

    ###################################################################
    # PHASE ONE  -  MAP x ∈ K_2 to ECOS form, except SOC/Exp

    # If a variable is in :Zero cone, fix at 0 with equality constraint.
    for j = 1:num_vars
        idxcone[j] != :Zero && continue

        new_row    = zeromat(1,num_vars)
        new_row[j] = 1.0
        ecos_A     = vcat(ecos_A, new_row)
        ecos_b     = vcat(ecos_b, 0.0)
        m.col_map_ind[rev_map[j]] = length(ecos_b)
    end

    # G matrix:
    # * 1 row ∀ :NonNeg & :NonPos cones
    # * 1 row ∀ variable in SOC/Exp cones
    # We will only handle the first case here, the rest in phase 3.
    num_G_row_negpos = 0
    for j = 1:num_vars
        !(idxcone[j] == :NonNeg || idxcone[j] == :NonPos) && continue
        num_G_row_negpos += 1
        m.col_map_ind[rev_map[j]] = num_G_row_negpos
    end
    ecos_G = zeromat(num_G_row_negpos,num_vars)
    ecos_h =   zeros(num_G_row_negpos)

    # Handle the :NonNeg, :NonPos cases
    num_pos_orth = 0
    G_row = 1
    for j = 1:num_vars
        if idxcone[j] == :NonNeg
            ecos_G[G_row,j] = -1.0
            G_row += 1
            num_pos_orth += 1
        elseif idxcone[j] == :NonPos
            ecos_G[G_row,j] = +1.0
            G_row += 1
            num_pos_orth += 1
        end
    end
    @assert G_row == num_pos_orth + 1

    ###################################################################
    # PHASE TWO  -  MAP b-Ax ∈ K_1 to ECOS form, except SOC/Exp

    # Zero rows for Ax=b, NonNegPos rows to append to G,h
     eq_rows = Int[]
    pos_rows = Int[]
    neg_rows = Int[]
    for (cone,idxs) in constr_cones
        cone == :Free && error("Free cone constraints not handled")
        (cone == :SOC || cone == :ExpPrimal) && continue  # Handle later
        idxset = vec(collect(idxs))
        if cone == :Zero
            append!(eq_rows, idxset)
            continue
        end
        cone == :NonNeg && append!(pos_rows, idxset)
        cone == :NonPos && append!(neg_rows, idxset)
    end
    # Update mappings - eq, nonneg, nonpos
     eq_cur_ind = length(ecos_b) + 1
     eq_cur_ind = update_row_map(:Zero, eq_cur_ind)
    neq_cur_ind = length(ecos_h) + 1
    neq_cur_ind = update_row_map(:NonNeg, neq_cur_ind)
    neq_cur_ind = update_row_map(:NonPos, neq_cur_ind)
    # Equality constraints / Zero cones
    ecos_A = vcat(ecos_A,  A[eq_rows,rev_map])
    ecos_b = vcat(ecos_b,  b[eq_rows])
    ecos_G = vcat(ecos_G,  A[pos_rows,rev_map])
    ecos_h = vcat(ecos_h,  b[pos_rows])
    ecos_G = vcat(ecos_G, -A[neg_rows,rev_map])  # b-a'x <= 0 - flip sign,
    ecos_h = vcat(ecos_h, -b[neg_rows])          # then maps to a row in h-Gx
    G_row        += length(pos_rows) + length(neg_rows)
    num_pos_orth += length(pos_rows) + length(neg_rows)

    ###################################################################
    # PHASE THREE  -  MAP SOC variables and constraints to ECOS form

    # Handle the SOC variable cones
    # MPB  form: vector of var (y,x) is in the SOC ||x|| <= y
    # ECOS form: h - Gx ∈ Q  -->  0 - Ix ∈ Q
    num_G_row_soc = 0
    num_G_row_exp = 0
    for j = 1:num_vars
        if idxcone[j] == :SOC
            num_G_row_soc += 1
        elseif idxcone[j] == :ExpPrimal
            num_G_row_exp += 1
        end
    end
    ecos_G = vcat(ecos_G, zeromat(num_G_row_soc,num_vars))
    ecos_h = vcat(ecos_h,   zeros(num_G_row_soc))

    num_SOC_cones = 0
    SOC_conedims  = Int[]
    for (cone, idxs) in var_cones
        cone != :SOC && continue
        # Found a new SOC
        num_SOC_cones += 1
        push!(SOC_conedims, length(idxs))
        # Add the entries (carrying on from pos. orthant)
        for j in idxs
            ecos_G[G_row,fwd_map[j]] = -1.0
            m.col_map_ind[j] = G_row
            G_row += 1
        end
    end
    @assert G_row == num_pos_orth + num_G_row_soc + 1

    # Handle the SOC constraint cones
    # Collect all the rows we'll be appending to G,h
    all_rows = Int[]
    for (cone,idxs) in constr_cones
        if cone == :SOC
            num_SOC_cones += 1
            push!(SOC_conedims, length(idxs))
            idx_list   = collect(idxs)
            all_rows   = vcat(all_rows,   idx_list)
        end
    end
    neq_cur_ind = update_row_map(:SOC, neq_cur_ind)
    ecos_G = vcat(ecos_G, A[all_rows,rev_map])
    ecos_h = vcat(ecos_h, b[all_rows])

    ###################################################################
    # PHASE FOUR  -  MAP Exp variables and constraints to ECOS form

    # Note: The MPB definition of the ExpPrimal cone and the ECOS
    # defintions swap the 2nd and 3rd components, so we have to
    # manually permute the indices.

    num_exp_cones = 0
    ecos_G = vcat(ecos_G, zeromat(num_G_row_exp,num_vars))
    ecos_h = vcat(ecos_h,   zeros(num_G_row_exp))

    G_row += sum(all_rows)
    for (cone, idxs) in var_cones
        cone != :ExpPrimal && continue
        # Found a new Exp
        num_exp_cones += 1
        @assert length(idxs) == 3
        # Add the entries in permuted order
        ecos_G[G_row,fwd_map[idxs[1]]] = -1.0
        ecos_G[G_row+2,fwd_map[idxs[2]]] = -1.0
        ecos_G[G_row+1,fwd_map[idxs[3]]] = -1.0
        m.col_map_ind[idxs[1]] = G_row
        m.col_map_ind[idxs[2]] = G_row+2
        m.col_map_ind[idxs[3]] = G_row+1
        G_row += 3
    end
    @assert G_row == num_pos_orth + num_G_row_soc + sum(all_rows) + num_G_row_exp + 1

    exp_rows = Int[]
    for (cone,idxs) in constr_cones
        if cone == :ExpPrimal
            num_exp_cones += 1
            @assert length(idxs) == 3
            append!(exp_rows, [idxs[1],idxs[3],idxs[2]])

            # override update_row_map to handle the permutation
            m.row_map_type[idxs] .= :ExpPrimal
            m.row_map_ind[idxs[1]] = neq_cur_ind
            m.row_map_ind[idxs[3]] = neq_cur_ind+1
            m.row_map_ind[idxs[2]] = neq_cur_ind+2
            neq_cur_ind += 3
        end
    end
    ecos_G = vcat(ecos_G, A[exp_rows,rev_map])
    ecos_h = vcat(ecos_h, b[exp_rows])

    ###################################################################
    # Store in the ECOS structure
    m.nvar          = num_vars          # Num variable
    m.nineq         = size(ecos_G,1)    # Num inequality constraints
    m.neq           = length(ecos_b)    # Num equality constraints
    m.npos          = num_pos_orth      # Num ineq. constr. in +ve orthant
    m.ncones        = num_SOC_cones     # Num second-order cones
    m.conedims      = SOC_conedims      # Num contr. in each SOC
    m.nexp_cones    = num_exp_cones     # Num exponential cones
    m.G             = ECOSMatrix(ecos_G)
    m.A             = ECOSMatrix(ecos_A)
    m.c             = ecos_c
    m.orig_sense    = :Min
    m.h             = ecos_h
    m.b             = ecos_b
    m.fwd_map       = fwd_map           # Used to return solution
end

function MPB.getdual(m::ECOSMathProgModel)
    duals = zeros(length(m.row_map_ind))
    for (mpb_row,ecos_row) in enumerate(m.row_map_ind)
        cone = m.row_map_type[mpb_row]
        if cone == :Zero
            # This MPB constraint ended up in ECOS equality block
            duals[mpb_row] = m.dual_sol_eq[ecos_row]
        else
            # Ended up in ECOS inequality block
            if cone == :NonPos
                duals[mpb_row] = -m.dual_sol_ineq[ecos_row]
            else
                duals[mpb_row] = m.dual_sol_ineq[ecos_row]
            end
        end
    end
    return duals
end

function MPB.getvardual(m::ECOSMathProgModel)
    duals = zeros(length(m.col_map_ind))
    for (mpb_col,ecos_row) in enumerate(m.col_map_ind)
        cone = m.col_map_type[mpb_col]
        cone == :Free && continue # dual is zero
        if cone == :Zero
            # This MPB zero var ended up in ECOS equality block
            duals[mpb_col] = m.dual_sol_eq[ecos_row]
        else
            # Ended up in ECOS inequality block
            if cone == :NonPos
                duals[mpb_col] = -m.dual_sol_ineq[ecos_row]
            else
                duals[mpb_col] = m.dual_sol_ineq[ecos_row]
            end
        end
    end
    return duals
end

MPB.getsolvetime(m::ECOSMathProgModel) = m.solve_time


text:
using JuMP, Mosek, MosekTools, LinearAlgebra, Suppressor, StatsBase, Ipopt

function getMatrixTraceNorm_EDM_noisy(coords, entries, D, n, numEntriesSampled, ys, num_anchors, useRadiorange=false, radiorange=1.0)
    m=Model(Mosek.Optimizer)
    @variable(m, G[1:n+num_anchors+2, 1:n+num_anchors+2], PSD)
    @constraint(m, G[1:2,1:2].==[1.0 0.0; 0.0 1.0])
    # Note that G_full=[I X;X' G], so we can recover X from G_full quite simply when rank(G)=2 in the underlying relaxation.
    @variable(m, slack[1:numEntriesSampled])
    @variable(m, absslack[1:numEntriesSampled])
    @constraint(m, absslack.>=slack)
    @constraint(m, absslack.>=-slack)
    mu=100.0

    for t=1:Int(numEntriesSampled)
        if (coords[entries[t],1]<=num_anchors) && (coords[entries[t],2]<=num_anchors)
            #don't need to consider this case as doesn't yield any new information
        elseif (coords[entries[t],1]<=num_anchors) && (coords[entries[t],2]>num_anchors)
            e_k=zeros(n+num_anchors)
            e_k[coords[entries[t],2]]=1.0
            @constraint(m, Compat.dot([ys[coords[entries[t],1],:]; -e_k]*[ys[coords[entries[t],1],:]; -e_k]', G).==D[coords[entries[t],1], coords[entries[t],2]]+slack[t])
        elseif (coords[entries[t],1]>num_anchors) && (coords[entries[t],2]<=num_anchors)
            e_k=zeros(n+num_anchors)
            e_k[coords[entries[t],1]]=1.0
            @constraint(m,Compat.dot([ys[coords[entries[t],2],:]; -e_k]*[ys[coords[entries[t],2],:]; -e_k]', G).==D[coords[entries[t],1], coords[entries[t],2]]+slack[t])
        elseif (coords[entries[t],1]>num_anchors) && (coords[entries[t],2]>num_anchors)
            e_i=zeros(n+num_anchors)
            e_i[coords[entries[t],1]]=1.0
            e_k=zeros(n+num_anchors)
            e_k[coords[entries[t],2]]=1.0
            @constraint(m, Compat.dot([zeros(2); e_i.-e_k]*[zeros(2); e_i.-e_k]', G).==D[coords[entries[t],1], coords[entries[t],2]]+slack[t])
        end
    end
    if useRadiorange
        altentries=setdiff((1:(n+num_anchors)^2), entries)
        for t=1:((n+num_anchors)^2 -Int(numEntriesSampled))
            if (coords[altentries[t],1]<=num_anchors) && (coords[altentries[t],2]<=num_anchors)
                #don't need to consider this case as doesn't yield any new information
            elseif (coords[altentries[t],1]<=num_anchors) && (coords[altentries[t],2]>num_anchors)
                e_k=zeros(n+num_anchors)
                e_k[coords[altentries[t],2]]=1.0
                @constraint(m, Compat.dot([ys[coords[altentries[t],1],:]; -e_k]*[ys[coords[altentries[t],1],:]; -e_k]', G).>=radiorange^2)
            elseif (coords[altentries[t],1]>num_anchors) && (coords[altentries[t],2]<=num_anchors)
                e_k=zeros(n+num_anchors)
                e_k[coords[altentries[t],1]]=1.0
                @constraint(m,Compat.dot([ys[coords[altentries[t],2],:]; -e_k]*[ys[coords[altentries[t],2],:]; -e_k]', G).>=radiorange^2)
            elseif (coords[altentries[t],1]>num_anchors) && (coords[altentries[t],2]>num_anchors)
                e_i=zeros(n+num_anchors)
                e_i[coords[altentries[t],1]]=1.0
                e_k=zeros(n+num_anchors)
                e_k[coords[altentries[t],2]]=1.0
                @constraint(m, Compat.dot([zeros(2); e_i.-e_k]*[zeros(2); e_i.-e_k]', G).>=radiorange^2)
            end
        end
    end

    @objective(m, Min, sum(G[i,i] for i=1:n)+mu*sum(absslack))

    @suppress optimize!(m)
   return value.(G[3+num_anchors:end, 1:2]), value.(G[3+num_anchors:end, 3+num_anchors:end]), value.(G[3:end, 1:2]) #return X
end

function getMatrixTraceNorm_EDM(coords, entries, D, n, numEntriesSampled, ys, num_anchors, imposeTriangleInequalities=false)
    m=Model(Mosek.Optimizer)
    @variable(m, G[1:n+num_anchors+2, 1:n+num_anchors+2], PSD)
    @constraint(m, G[1:2,1:2].==[1.0 0.0; 0.0 1.0])
    # Note that G_full=[I X;X' G], so we can recover X from G_full quite simply when rank(G)=2 in the underlying relaxation.

    for t=1:Int(numEntriesSampled)
        if (coords[entries[t],1]<=num_anchors) && (coords[entries[t],2]<=num_anchors)
            #don't need to consider this case as doesn't yield any new information
        elseif (coords[entries[t],1]<=num_anchors) && (coords[entries[t],2]>num_anchors)
            e_k=zeros(n+num_anchors)
            e_k[coords[entries[t],2]]=1.0
            @constraint(m, Compat.dot([ys[coords[entries[t],1],:]; -e_k]*[ys[coords[entries[t],1],:]; -e_k]', G).==D[coords[entries[t],1], coords[entries[t],2]])
        elseif (coords[entries[t],1]>num_anchors) && (coords[entries[t],2]<=num_anchors)
            e_k=zeros(n+num_anchors)
            e_k[coords[entries[t],1]]=1.0
            @constraint(m,Compat.dot([ys[coords[entries[t],2],:]; -e_k]*[ys[coords[entries[t],2],:]; -e_k]', G).==D[coords[entries[t],1], coords[entries[t],2]])
        elseif (coords[entries[t],1]>num_anchors) && (coords[entries[t],2]>num_anchors)
            e_i=zeros(n+num_anchors)
            e_i[coords[entries[t],1]]=1.0
            e_k=zeros(n+num_anchors)
            e_k[coords[entries[t],2]]=1.0
            @constraint(m, Compat.dot([zeros(2); e_i.-e_k]*[zeros(2); e_i.-e_k]', G).==D[coords[entries[t],1], coords[entries[t],2]])
        end
    end
    @objective(m, Min, sum(G[i,i] for i=1:n))

    @suppress optimize!(m)
   return value.(G[3+num_anchors:end, 1:2]), value.(G[3+num_anchors:end, 3+num_anchors:end]), value.(G[3:end, 1:2]) #return X
end

function getMatrixTraceNorm_EDM_Gram(coords, entries, D, n, numEntriesSampled, d_radio, lambda, gamma, k)
    m=Model(Mosek.Optimizer)
    @variable(m, G[1:n, 1:n], PSD)
    @variable(m, G2[1:n, 1:n])
    @variable(m, theta[1:n, 1:n], Symmetric)
    @variable(m, Y[1:n, 1:n], Symmetric)
    @variable(m, xi[1:n, 1:n])
    @variable(m, absxi[1:n, 1:n])
    @constraint(m, absxi.>=xi)
    @constraint(m, absxi.>=-xi)
    @constraint(m, G.==G2) #Since Mosek can't have a matrix which is part of two seperate PSD constraints at once

    for t=1:Int(numEntriesSampled)
        it=coords[entries[t],1]
        jt=coords[entries[t],2]
        @constraint(m, G[it,it]+G[jt,jt]-2*G[it,jt]+xi[it,jt]==D[it,jt])
    end

    @constraint(m, Symmetric(Matrix(1.0I, n, n)-Y) in PSDCone());
    @constraint(m, Symmetric([Y G2; G2' theta]) in PSDCone());
    @constraint(m, sum(Y[i,i] for i=1:n)<=k)

    @objective(m, Min, sum(G[i,i] for i=1:n)+lambda*sum(absxi)+(1/(2*gamma))*sum(theta[i,i] for i=1:n))#sum(G[i,j]^2 for i=1:n for j=1:n))
    @show d_radio, lambda, gamma
    @suppress optimize!(m)
    @show objective_value(m)
    # @show value.(G)
    # @show value.(Y)
   return value.(G), value.(Y)
end

max_stars_repo_path: Euclidean Distance Matrix Completion/convex_penalties_edm.jl

max_stars_repo_name: ryancorywright/MixedProjectionSoftware

max_issues_repo_name: ryancorywright/MixedProjectionSoftware

max_forks_repo_name: ryancorywright/MixedProjectionSoftware


text:
module LeaveOneOut

    using ParetoSmooth
    using Statistics
    using StructArrays
    using Tullio
    using VectorizationBase

    export psis_loo, PsisLoo, Loo, AbstractLoo

    const LOO_METHODS = ["psis"]

    function psis_loo(log_likelihood::ArrayType,
        score::Function=identity;  # Default to log_likelihood (elpd) score function
        rel_eff=similar(loo, 0)
        ) where {F<:AbstractFloat, ArrayType<:AbstractArray{F, 3}}

        dims = size(log_likelihood)
        data_size = dims[1]
        mcmc_count = dims[2] * dims[3]  # total number of samples from posterior
        # score_array::ArrayType = similar(log_likelihood)
        # score_array .= score(log_likelihood)
        score_array = log_likelihood

        psis_object = psis(log_likelihood, rel_eff)
        weights = psis_object.weights
        ξ = psis_object.pareto_k
        ess = psis_object.ess

        @tullio pointwise_ev[i] := weights[i,j,k] * score_array[i,j,k]
        # Replace with quantile mcse estimate from R LOO package?
        @tullio pointwise_mcse[i] := sqrt <|
            (weights[i,j,k] * score_array[i,j,k] .- pointwise_ev[i])^2 ./ (ess[i] - 1)
        @tullio pointwise_naive[i] := log_likelihood[i,j,k] |> _ / mcmc_count
        pointwise_p_eff = pointwise_ev - pointwise_naive
        points = (
            estimate = pointwise_ev,
            mcse = pointwise_mcse,
            p_eff = pointwise_p_eff,
            pareto_k = ξ
            )
        pointwise = StructArray{LooPoint{F}}(points)

        @tullio ev := points.estimate[i]
        @tullio ev_naive := pointwise_naive[i]
        p_eff = ev - ev_naive

        ev_se = std(pointwise_ev; mean=ev) / sqrt(data_size)
        p_eff_se = std(pointwise_p_eff; mean=p_eff) / sqrt(data_size)

        vals = Dict("Score Est"=>ev, "Parameters (Eff)"=>p_eff,
                     "SE(Score Est)"=>ev_se, "SE(Parameters)"=>p_eff_se
                    )

        return PsisLoo(vals, pointwise, psis_object)

    end

    # Pointwise entries
    struct LooPoint{F<:AbstractFloat}
        estimate::F
        mcse::F
        p_eff::F
        pareto_k::F
    end

    abstract type AbstractLoo end

    struct PsisLoo{
        F<:AbstractFloat,
        AF<:AbstractArray{F},
        VF<:AbstractVector{F},
        I<:Integer,
        VI<:AbstractVector{I},
    } <: AbstractLoo
        estimates::Dict{String,F}
        pointwise::StructArray{LooPoint{F}}
        psis_object::Psis{F,AF,VF,I,VI}
    end

end

max_stars_repo_name: ParadaCarleton/Bombe

max_issues_repo_name: ParadaCarleton/Bombe

max_forks_repo_name: ParadaCarleton/Bombe


text:
using Test

#NOTE: @testset will not work for testing this.
using StructuralInheritance
## TEST BASIC STRUCTURAL INHERITENCE ##

SI = StructuralInheritance

  StructuralInheritance.@protostruct struct A
     f_a::Int
     f_b
  end

  StructuralInheritance.@protostruct struct B <: A
      f_c::Float32
      f_d
  end

  @test fieldnames(B) == (:f_a,:f_b,:f_c,:f_d) #Test names of fields and order
  @test fieldtype.(B,[1,2,3,4]) == [Int,Any,Float32,Any]
  @test B <: ProtoB && B <: ProtoA

# exception thrown when inherited class uses same names
  @test_throws Any StructuralInheritance.@protostruct struct C <: A
      f_a::Int
  end

  @test_throws Any StructuralInheritance.@protostruct struct D <: A
      f_b
  end

  # exception thrown trying to inherit from a concrete class not defined
  # by @protostruct
  @test_throws Any StructuralInheritance.@protostruct struct E <: Int
      f_b
  end

  @test_throws Any StructuralInheritance.@protostruct struct E <: Int
      f_b
  end

  #inheritence from any abstract class
  StructuralInheritance.@protostruct struct F <: Real
      f_b
  end

  @test F(3).f_b == 3

  StructuralInheritance.@protostruct struct G <: B
      f_e::Float32
      f_f
  end

  @test fieldnames(G) == (:f_a,:f_b,:f_c,:f_d,:f_e,:f_f) #Test names of fields and order
  @test fieldtype.(G,[1,2,3,4,5,6]) == [Int,Any,Float32,Any,Float32,Any]
  @test G <: ProtoG && G <: ProtoB && G <: ProtoA

#TEST MODULE SANITIZATION FACILITY
module MA
  using Main.StructuralInheritance
  @protostruct struct A
    f_a_MA::Int
  end
  @protostruct struct B <: A
    f_a::A
  end
end

@protostruct struct H <: MA.B
  f_c::A
end

@protostruct struct H_2 <: MA.B
  f_c::MA.A
end

@test fieldnames(H_2) == (:f_a_MA,:f_a,:f_c)
@test fieldtype.(H_2,[1,2,3]) == [Int,MA.A,MA.A]

@test fieldnames(H) == (:f_a_MA,:f_a,:f_c)
@test fieldtype.(H,[1,2,3]) == [Int,MA.A,A]

@test StructuralInheritance.@protostruct(struct I <: A
    f_c::Float32
    f_d
end,"ProtoType") == ProtoTypeI

@test StructuralInheritance.@protostruct(struct J
    f_c::Float32
    f_d
end,:("ProtoType")) == ProtoTypeJ

@test_throws Any StructuralInheritance.@protostruct(struct K
    f_c::Float32
    f_d
end,"") == ProtoTypeK

#TEST PARAMETRIC INHERITENCE
module M_paramfields
    using Main.StructuralInheritance
    using Test
    @test @protostruct(struct A
        f_a::Core.Main.Base.Array{Float16,3}
    end) == ProtoA

    @test @protostruct(struct B <: A
        f_b::Complex{Float64}
    end) == ProtoB

    @test @protostruct(struct C <: B
        field_c::Array{Float32,4}
    end) == ProtoC

    @test fieldnames(C) == (:f_a,:f_b,:field_c)
    @test fieldtype.(C,[1,2,3]) == [Array{Float16,3},Complex{Float64},Array{Float32,4}]

    @protostruct struct D{X,B}
        f1::Array{X,B}
    end

    @protostruct struct E{T} <: D{Complex{T},2} end

    @protostruct struct F <: E{Int} end

    @test fieldnames(F) == (:f1,)
    @test fieldtype(F,1) == Array{Complex{Int64},2}
end

@test StructuralInheritance.@protostruct(struct K{T}
    f_a::T
end) == ProtoK

@test fieldtype(K{Int},1) == Int
@test fieldnames(K) == (:f_a,)

@test @protostruct(struct L{T} <: K{T}
    f_b::T
end) == ProtoL

@test fieldnames(L) == (:f_a,:f_b)
@test fieldtype.(L{Real},[1,2]) == [Real,Real]

@test @protostruct(struct N{R} <: K{R}
    f_b::R
end) == ProtoN

@test fieldnames(N) == (:f_a,:f_b)
@test fieldtype.(N{Complex},[1,2]) == [Complex,Complex]

@protostruct mutable struct BB
        f_a::Int
end

@protostruct mutable struct CC{A,BA} <: BB
    f_b::A
    f_c::BA
end

@test @protostruct( mutable struct DD{C,D} <: CC{D,Complex}
    f_d::C
end) == ProtoDD

@test @protostruct( mutable struct DD_2{C,D} <: CC{D,Complex{C}}
    f_d::C
end) == ProtoDD_2

@test @protostruct( mutable struct O <: DD{Int,Real}
    f_e::Complex
end) == ProtoO

@test @protostruct( mutable struct O_2 <: DD_2{Int,Real}
    f_e::Complex
end) == ProtoO_2

@test fieldnames(O) == (:f_a,:f_b,:f_c,:f_d,:f_e)
@test fieldtype.(O,[1,2,3,4,5]) == [Int,Real,Complex,Int,Complex]

@test fieldnames(O_2) == (:f_a,:f_b,:f_c,:f_d,:f_e)
@test fieldtype.(O_2,[1,2,3,4,5]) == [Int,Real,Complex{Int},Int,Complex]

module M_literal_func
    using Main.StructuralInheritance
    f(x) = Int
    f(x::Array) = length(x) == 2 ? Complex : Array
    @protostruct struct P
        f1::f(2)
        f2::f([])
        f3::f([1,2])
    end
end

@protostruct struct Q <: M_literal_func.P
    f4::Real
end

@test fieldnames(Q) == (:f1,:f2,:f3,:f4)
@test fieldtype.(Q,[1,2,3,4]) == [Int,Array,Complex,Real]

@test SI.totuple(4 + 5im) == (4,5)
@test SI.totuple(4) == (4,)

@protostruct struct R
    ff::Int
    sf
    R(x) = new(x,x^2)
    R(x,y) = new(x,y)
end

@protostruct struct S <: R
    tf::Int
end

S(x) = S(SI.totuple(R(x))...,x^3)

@test S(2) == S(2,4,8)

module mutabilityTestings
    using Main.StructuralInheritance
    using Test
    @protostruct struct A_im; end;
    @protostruct mutable struct A_m; end;
    @test_throws Any @protostruct mutable struct B_m <: A_im; end;
    @test_throws Any @protostruct struct B_im <: A_m; end;
    @test_throws Any @protostruct mutable struct B_m <: A_im
    end "Proto" "hat"

    @test @protostruct(mutable struct B_m <: A_im end,
                       "Sudo",
                       true) == SudoB_m
    @test @protostruct(struct B_im <: A_m; end,
                       "A",
                       true) == AB_im;
end

module macroLiteral
    using Main.StructuralInheritance
    using Test
    module Inner
        using Main.StructuralInheritance

        macro m(x)
             x
         end
         @protostruct struct A
             f::@m(Int)
         end
    end

    @protostruct struct B <: Inner.A
        f2::Inner.@m Float64
    end

    @test fieldnames(B) == (:f,:f2)
    @test fieldtype.(B,[1,2]) == [Int,Float64]

end

module ParametricTypeConstraints
    using Test
    using Main.StructuralInheritance

    @protostruct struct A{T<:Number}
        ffa::T
    end

    @test_throws Any A("hmm")

    @test A(3).ffa == 3

    @protostruct struct B{T<:Real,A<:Array} <: A{T}
        ffb::A
    end

    @protostruct struct C{T,A,L} <: B{T,A}
        ffc::L
    end

    @test_throws Any C(4+3im,[],"huh")
    @test_throws Any C(2,4,"huh")

    @test StructuralInheritance.totuple(C(1,[4],"huh")) == (1,[4],"huh")

end

module MacroFields
    using Test
    using Main.StructuralInheritance

    macro e1()
        esc(quote f1::Int; f2::X; end)
    end

    macro e2()
        esc(quote f4::K; end)
    end

    macro e3()
        esc(quote; f5::Int; f6::T; end)
    end

    @protostruct struct A{X}
        @e1
    end

    @protostruct struct B{X,T,K} <: A{Array{X,T}}
        @e2()
    end

    @protostruct struct C{T} <: B{Complex{T},2,T}
        @e3
    end

    @protostruct struct D <: C{Float32} end

    fieldnames(D) == (:f1,:f2,:f4,:f5,:f6)
    fieldtype.(D,(1,2,3,4,5)) == (Int,Array{Complex{Float32},2},Float32,Int,Float32)
end

#whitebox tests
@test SI.qualifyname(Int64) == :(Core.Int64)
@test SI.inherits(:x) == false
@test SI.inherits(:(x <: y)) == true
@test SI.inherits(:(x{z} <: y)) == true
@test SI.inherits(:(x{y})) == false

@test SI.qualifyname(:Xib,[:c,:v,:l,:k]) == :(c.v.l.k.Xib)
@test SI.qualifyname(:(M.P.Xib),[:c,:v,:l,:k]) == :(c.v.l.k.M.P.Xib)

@test SI.detypevar(:x) == :x

@test SI.iscontainerlike(:(x.v)) == false
@test SI.iscontainerlike(4) == false

@test SI.ispath(:x) == false
@test SI.ispath(:(x.c)) == true

@test SI.flattenfields("str") == SI.FieldType[]
@test SI.flattenfields(:x) == SI.FieldType[:x]


text:
julia> "# Factorial"
       fac = f -> (n -> n < 2 ? 1 : n * f(n - 1))

julia> "# Fibonacci"
       fib = f -> (n -> n == 0 ? 0 : (n == 1 ? 1 : f(n - 1) + f(n - 2)))

julia> [Y(fac)(i) for i = 1:10]
10-element Array{Any,1}:
       1
       2
       6
      24
     120
     720
    5040
   40320
  362880
 3628800

julia> [Y(fib)(i) for i = 1:10]
10-element Array{Any,1}:
  1
  1
  2
  3
  5
  8
 13
 21
 34
 55

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
    ILMM(f, H)

An Instantaneous Linear Mixing Model (ILMM) -- a distribution over vector-
valued functions. Let `p` be the number of observed outputs, and `m` the number of latent
processes, then `H`, also known as the mixing matrix, is a `p x m` matrix whose column
space spans the output space. Given a sample from `f`, the value of the ILMM at any input `x` for all outputs is
```julia
H * [f((x, 1)), f((x, 2)), ..., f((x, m))]
```

# Arguments:
- f: a multi-output AbstractGP with `m` outputs.
- H: a `p x m` matrix representing a fixed basis of our p-dim target: h_1,...,h_m
"""
struct ILMM{Tf<:AbstractGP,TH<:AbstractMatrix{<:Real}} <: AbstractGP
    f::Tf
    H::TH
end

"""
    get_latent_gp(f::ILMM)

Returns the underlying latent space AbstractGP belonging to `f`.

```jldoctest
julia> f = ILMM(independent_mogp([GP(SEKernel())]), rand(2,2));

julia> latent_f = get_latent_gp(f);

julia> latent_f isa IndependentMOGP
true

julia> latent_f.fs == [GP(SEKernel())]
true
```
"""
get_latent_gp(f::ILMM) = f.f

noise_var(Σ::Diagonal{<:Real,<:Fill}) = FillArrays.getindex_value(Σ.diag)

reshape_y(y::AbstractVector{<:Real}, N::Int) = reshape(y, N, :)'

function unpack(fx::FiniteGP{<:ILMM,<:MOInputIsotopicByOutputs,<:Diagonal{<:Real,<:Fill}})
    f = fx.f.f
    H = fx.f.H
    σ² = noise_var(fx.Σy)
    x = fx.x.x

    # Check that the number of outputs requested agrees with the model.
    fx.x.out_dim == size(H, 1) || throw(error("out dim of x != out dim of f."))
    return f, H, σ², x
end

"""
    project(H, σ²)

Computes the projection `T` and `ΣT` given the mixing matrix and noise.
"""
function project(H::AbstractMatrix{Z}, σ²::Z) where {Z<:Real}
    Σ = σ² * I
    ΣT_inv = H' * inv(Σ) * H + 1e-9I
    T = cholesky(Symmetric(ΣT_inv)) \ H' * inv(Σ) # Σ easily invertible, no problem there
    ΣT = T * Σ * T'

    return T, ΣT
end

"""
    rand(rng::AbstractRNG, fx::FiniteGP{<:ILMM})

Sample from the ILMM, including the observation noise.
Follows generative structure of model 2 from [1].
Follows the AbstractGPs.jl API.
[1] - Bruinsma et al 2020.
"""
function AbstractGPs.rand(rng::AbstractRNG, fx::FiniteGP{<:ILMM})
    f, H, σ², x = unpack(fx)
    p, m = size(H)

    x_mo_input = MOInputIsotopicByOutputs(x, m)

    latent_rand = rand(rng, f(x_mo_input, 1e-12))
    ε = randn(rng, length(fx))
    return vec(reshape(latent_rand, length(x), m) * H') + sqrt(σ²) .* ε
end

# See AbstractGPs.jl API docs.
function AbstractGPs.rand(rng::AbstractRNG, fx::FiniteGP{<:ILMM}, N::Int)
    return reduce(hcat, [rand(rng, fx) for _ in 1:N])
end

# See AbstractGPs.jl API docs.
function Distributions._rand!(
    rng::AbstractRNG,
    fx::FiniteGP{<:ILMM},
    y::AbstractVecOrMat{<:Real}
)
    N = size(y, 2)
    if N == 1
        y .= AbstractGPs.rand(rng, fx)
    else
        y .= AbstractGPs.rand(rng, fx, N)
    end
end

function _intermediate_mean_and_var_and_cov_quantities(fx::FiniteGP{<:ILMM})
    f, H, σ², x = unpack(fx)
    p, m = size(H)
    n = length(x)

    x_mo_input = MOInputIsotopicByOutputs(x, m)

    latent_mean, latent_cov = mean_and_cov(f(x_mo_input))

    H_full = kron(H, Matrix(I, n, n))
    return H_full, latent_mean, latent_cov, σ²
end

# See AbstractGPs.jl API docs.
function AbstractGPs.mean_and_var(fx::FiniteGP{<:ILMM})
    H_full, latent_mean, latent_cov, σ² = _intermediate_mean_and_var_and_cov_quantities(fx)

    M = H_full * latent_mean
    V = AbstractGPs.diag_Xt_A_X(cholesky(latent_cov), H_full') .+ σ²

    return collect(vec(M)), V
end

# See AbstractGPs.jl API docs.
function AbstractGPs.mean_and_cov(fx::FiniteGP{<:ILMM})
    H_full, latent_mean, latent_cov, σ² = _intermediate_mean_and_var_and_cov_quantities(fx)

    M = H_full * latent_mean
    C = AbstractGPs.Xt_A_X(cholesky(latent_cov), H_full') + σ² * I

    return collect(vec(M)), C
end

# See AbstractGPs.jl API docs.
AbstractGPs.mean(fx::FiniteGP{<:ILMM}) = mean_and_var(fx)[1]

# See AbstractGPs.jl API docs.
AbstractGPs.var(fx::FiniteGP{<:ILMM}) = mean_and_var(fx)[2]

AbstractGPs.cov(fx::FiniteGP{<:ILMM}) = mean_and_cov(fx)[2]

# See AbstractGPs.jl API docs.
function AbstractGPs.logpdf(fx::FiniteGP{<:ILMM}, y::AbstractVector{<:Real})
    f, H, σ², x = unpack(fx)
    p, m = size(H)
    n = length(x)

    # # Projection step.
    Y = reshape_y(y, length(x))
    T, ΣT = project(H, σ²)
    Ty = RowVecs((T*Y)')
    Xproj, Yproj = prepare_isotopic_multi_output_data(x, Ty)
    ΣT = kron(ΣT, Matrix(I,n,n))

    return AbstractGPs.logpdf(f(Xproj, ΣT), Yproj) + regulariser(fx, Y)
end

"""
    regulariser(fx, y)

Computes the regularisation term of the logpdf.
See e.g. appendix A.4 of [1] - Bruinsma et al 2020.
"""
function regulariser(fx, Y::AbstractMatrix{<:Real})
    fs, H, σ², x = unpack(fx)
    p, m = size(H)
    n = size(Y, 2)

    # Projection step.
    T, ΣT = project(H, σ²)

    return -(n * ((p - m) * log(2π) + (p * log(σ²) - logdet(ΣT))) +
        sum((1/σ²) * abs2.(Y .- H*T*Y))) / 2
end

# See AbstractGPs.jl API docs.
function AbstractGPs.posterior(fx::FiniteGP{<:ILMM}, y::AbstractVector{<:Real})
    f, H, σ², x = unpack(fx)
    p, m = size(H)
    n = length(x)

    # # Projection step.
    Y = reshape_y(y, length(x))
    T, ΣT = project(H, σ²)
    Ty = RowVecs((T*Y)')
    Xproj, Yproj = prepare_isotopic_multi_output_data(x, Ty)
    ΣT = kron(ΣT, Matrix(I,n,n))

    f_posterior = AbstractGPs.posterior(f(Xproj, ΣT), Yproj)
    return ILMM(f_posterior, H)
end


text:
module ArDCA

using Random: randperm
using Distributed: @distributed
using Printf: @printf
using LinearAlgebra: rmul!, norm
using ExtractMacro: @extract
using NLopt: Opt,ftol_abs!,xtol_rel!,xtol_abs!,ftol_rel!,maxeval!,min_objective!,optimize
using Distributions: wsample
using LoopVectorization: @avx
using DCAUtils: read_fasta_alignment,remove_duplicate_sequences,compute_weights

export ardca,ArVar,ArAlg,ArNet,sample,sample_with_weights,epistatic_score,dms_single_site

include("types.jl")
include("ar.jl")
include("utils.jl")
include("dca.jl")

end # end module

max_issues_repo_name: pagnani/ArDCA

max_forks_repo_name: pagnani/ArDCA


text:
using HTTP

"""
Checks if the Neo4j database is accesible
"""
function check_accessibility(URI::String, port::String)
    r = HTTP.request("GET", "http://$(URI):$(port)")
    @assert r.status == 200 "Couldn't connect to the database"
end

"""
Returns the response to the query as a list of jolts
"""
function get_data(usr::String, passw::String, URI::String, port::String, query::String)
    transaction = HTTP.request(
        "POST",
        "http://$(usr):$(passw)@$(URI):$(port)/db/neo4j/tx",
        ["Content-Type" => "application/json", "Accept" => "application/vnd.neo4j.jolt"],
        "{\"statements\" : [ {
                \"statement\" : \"$(query)\"
            }
        ]}"
    )
    return split(String(transaction.body), "\n")
end


text: using YAAD, AutoGrad


text:
"""
    getDBRoot(fileName)
Get root node of .xml database
"""
function getDBRoot(fileName::String)::EzXML.Node
    xdoc = EzXML.readxml(fileName)
    dbRoot = EzXML.root(xdoc)
    return dbRoot
end

"""
    extractColumn(dbRoot, propName, n=0)
Extract a column of data providing the root node, and the property name as described in the .xml file. To access deeper levels use '/' e.g. "LiquidDensity/eqno"
This function corrects the type from String to Int64 or Float64.
"""
function extractColumn(dbRoot::EzXML.Node, propName::String, n::Int64 = 0)

    temp = EzXML.attributes.(findall(
        "//compounds/compound/" * propName,
        dbRoot,
    ))::Array{Array{EzXML.Node,1},1}

    if (n == 0)
        data = EzXML.nodecontent.(hcat(temp...)[end, :])
    else
        data = EzXML.nodecontent.(hcat(temp...)[n, :])
    end

    # type correction if possible, all values need to be of Int64/Float64 for conversion to take place
    if (sum(isnothing.(tryparse.(Int64, data))) == 0)
        data = parse.(Int64, data)
    elseif (sum(isnothing.(tryparse.(Float64, data))) == 0)
        data = parse.(Float64, data)
    end

    return data
end

"""
    getPropertyDataFrame(fileName)
Return a DataFrame containing component properties for a given filename
"""
function getPropertyDataFrame(fileName::String)::DataFrame # this function is not complete, still in progress
    dbRoot = getDBRoot(fileName)

    # some columns have a deeper structure, need to update column selected based on node depth
    columnNames = EzXML.nodename.(findall("//compounds/compound[1]/*", dbRoot))

    ID = extractColumn(dbRoot, "CompoundID")
    lengthID = length(ID)

    df = DataFrame(ID = ID)

    for name in columnNames
        data = extractColumn(dbRoot, name)
        (length(data) == lengthID) && (df[!, Symbol(name)] = data)
    end

    return df
end

@inline getScalar(props, propnum) = ScalarProperty(props[propnum]["units"], parse(Float64, props[propnum]["value"]))

function getVector(props, propnum)
    units = props[propnum]["units"]
    eqno =  parse(Int, EzXML.elements(props[propnum])[1]["value"])

    A_ = EzXML.findfirst("A", props[propnum])
    A = isnothing(A_) ? 0.0 : parse(Float64, A_["value"])
    B_ = EzXML.findfirst("B", props[propnum])
    B = isnothing(B_) ? 0.0 : parse(Float64, B_["value"])
    C_ = EzXML.findfirst("C", props[propnum])
    C = isnothing(C_) ? 0.0 : parse(Float64, C_["value"])
    D_ = EzXML.findfirst("D", props[propnum])
    D = isnothing(D_) ? 0.0 : parse(Float64, D_["value"])
    E_ = EzXML.findfirst("E", props[propnum])
    E = isnothing(E_) ? 0.0 : parse(Float64, E_["value"])

    Tmin_ = EzXML.findfirst("Tmin", props[propnum])
    Tmin = RangeBoundary(Tmin_["units"], parse(Float64, Tmin_["value"]))
    Tmax_ = EzXML.findfirst("Tmax", props[propnum])
    Tmax = RangeBoundary(Tmax_["units"], parse(Float64, Tmax_["value"]))

    return VectorProperty(units, eqno, A, B, C, D, E, Tmin, Tmax)
end

function extractComp(dbRoot, idx)
    tmp = EzXML.findfirst("//compounds/compound[" * string(idx) *"]", dbRoot)
    props = nodes(tmp)

    # LibraryIndex::Int
    libidx = parse(Int, props[2]["value"])
    # CompoundID::String
    cmpid = props[4]["value"]
    # StructureFormula::String
    struture = props[6]["value"]
    # Family::Int
    family = parse(Int, props[8]["value"])
    # CriticalTemperature::ScalarProperty
    Tc = getScalar(props, 10)
    # CriticalPressure::ScalarProperty
    Pc = getScalar(props, 12)
    # CriticalVolume::ScalarProperty
    Vc = getScalar(props, 14)
    # CriticalComperssibility::ScalarProperty
    Zc = getScalar(props, 16)
    # NormalBoilingPointTemperature::ScalarProperty
    Tb = getScalar(props, 18)
    # NormalMeltingPointTemperature::ScalarProperty
    Tf = getScalar(props, 20)
    # TripplePointTemperature::ScalarProperty
    T3 = getScalar(props, 22)
    # TripplePointPressure::ScalarProperty
    P3 = getScalar(props, 24)
    # MolecularWeight::ScalarProperty
    MW = getScalar(props, 26)
    # LiquidVolumeAtNormalBoilingPoint::ScalarProperty
    Vb = getScalar(props, 28)
    # AcentricityFactor::ScalarProperty
    omega = getScalar(props, 30)
    # SolubilityParameter::ScalarProperty
    lambda = getScalar(props, 32)
    # DipoleMoment::ScalarProperty
    mu = getScalar(props, 34)
    # HeatOfFormation::ScalarProperty
    hf = getScalar(props, 36)
    # GibbsEnergyOfFormation::ScalarProperty
    gf = getScalar(props, 38)
    # AbsEntropy::ScalarProperty
    sabs = getScalar(props, 40)
    # HeatOfFusionAtMeltingPoint::ScalarProperty
    hff = getScalar(props, 42)
    # HeatOfCombustion::ScalarProperty
    hc = getScalar(props, 44)
    # LiquidDensity::VectorProperty
    rhol = getVector(props, 46)
    # VaporPressure::VectorProperty
    Psat = getVector(props, 48)
    # HeatOfVaporization::VectorProperty
    hvap = getVector(props, 50)
    # LiquidHeatCapacityCp::VectorProperty
    cpliq = getVector(props, 52)
    # IdealHeatCapacityCp::VectorProperty
    cpig = getVector(props, 54)
    # SecondVirialCoefficient::VectorProperty
    virB = getVector(props, 56)
    # LiquidViscosity::VectorProperty
    muL = getVector(props, 58)
    # VaporViscosity::VectorProperty
    muV = getVector(props, 60)
    # LiquidThermalConductivity::VectorProperty
    kL = getVector(props, 62)
    # VaporThermalConductivity::VectorProperty
    kV = getVector(props, 64)
    # RPPHeatCapacityCp::VectorProperty
    cpRPP = getVector(props, 66)
    # RelativeStaticPermittivity::VectorProperty
    epsRel = getVector(props, 68)
    # AntoineVaporPressure::VectorProperty
    PAnt = getVector(props, 70)
    # LiquidViscosityRPS::VectorProperty
    muRPS = getVector(props, 72)
    # COSTALDVolume::ScalarProperty
    VCost = getScalar(props, 74)
    # DiameterLJ::ScalarProperty
    sigLJ = getScalar(props, 76)
    # EnergyLJ::ScalarProperty
    eLJ = getScalar(props, 78)
    # RacketParameter::ScalarProperty
    ZRa = getScalar(props, 80)
    # FullerVolume::ScalarProperty
    VFuller = getScalar(props, 82)
    # Parachor::ScalarProperty
    par = getScalar(props, 84)
    # SpecificGravity::ScalarProperty
    SG = getScalar(props, 86)
    # Charge::ScalarProperty
    charge = getScalar(props, 88)
    # CostaldAcentricFactor::ScalarProperty
    omegaCost = getScalar(props, 90)
    # WilsonVolume::ScalarProperty
    VWils = getScalar(props, 92)
    # ChaoSeaderAcentricFactor::ScalarProperty
    OmegaCS = getScalar(props, 94)
    # ChaoSeaderSolubilityParameter::ScalarProperty
    lambdaCS = getScalar(props, 96)
    # ChaoSeaderLiquidVolume::ScalarProperty
    VCS = getScalar(props, 98)
    # CAS::StringVal
    CAS = props[100]["value"]

    return PureComp(
        libidx,
        cmpid,
        struture,
        family,
        Tc,
        Pc,
        Vc,
        Zc,
        Tb,
        Tf,
        T3,
        P3,
        MW,
        Vb,
        omega,
        lambda,
        mu,
        hf,
        gf,
        sabs,
        hff,
        hc,
        rhol,
        Psat,
        hvap,
        cpliq,
        cpig,
        virB,
        muL,
        muV,
        kL,
        kV,
        cpRPP,
        epsRel,
        PAnt,
        muRPS,
        VCost,
        sigLJ,
        eLJ,
        ZRa,
        VFuller,
        par,
        SG,
        charge,
        omegaCost,
        VWils,
        OmegaCS,
        lambdaCS,
        VCS,
        CAS
    )
end

"""
_parse_property(compound_node)

parse a property and returns the proper struct from the ChemSep XML
"""
function parse_property(property)
    identifier = Symbol(property.name) #Symbol of the property name
    elems = length(EzXML.elements(property))
    if iszero(elems) #is value, float, int or String
        data = EzXML.attributes(property)
        name = data[1].content
        raw_value = data[end].content
        if length(data) == 3 #property is a float number
            units =  data[2].content
            value = parse(Float64,raw_value)
            return identifier,Property(name,units,value)
        elseif length(data) == 2 #value is String or Int
            if name in ("Index","Family") #only ints
                value = parse(Int,raw_value)
                return identifier,Property(name,"_",value)
            else
                return identifier,Property(name,"_",raw_value)
            end
        end
    elseif EzXML.firstelement(property).name == "eqno" #Variable Property
        data = EzXML.attributes(property)
        valuedata = EzXML.elements(property)
        len = length(valuedata)
        name = data[1].content
        units = data[2].content
        eqno = parse(Int,attributes(valuedata[1])[1].content)
        parameters = [parse(Float64,attributes(valuedata[i])[1].content) for i in 2:(len-2)]
        valid_range = [parse(Float64,attributes(valuedata[i])[2].content) for i in (len-1):len]
        value = VariableProperty(eqno,parameters,valid_range)
        return identifier,Property(name,units,value)
    elseif EzXML.firstelement(property).name == "group" #GroupContribution
        data = EzXML.attributes(property)
        name = data[1].content
        groups = EzXML.elements(property)
        len = length(groups)
        id_vector = Int[]
        value_vector = Int[]
        for group in groups
            raw_id,raw_value = EzXML.attributes(group)
            push!(id_vector,parse(Int,raw_id))
            push!(value_vector,parse(Int,raw_value))
        end
        _id = NTuple{len,Int}(id_vector)
        _value = NTuple{len,Int}(value_vector)
        value = GroupContribution{len}(_id,_value)
        return identifier,Property(name,"_",value)
    end
end

    """
    parse_compound(compound_node)

    parses the XML of a compound node. returns a Dict{Symbol,Property}

    """
function parse_compound(compound)
    symbols = Symbol[]
    properties = Property[]
    for property in EzXML.eachelement(compound)
        name, prop = parse_property(property)
        push!(symbols,name)
        push!(properties,prop)
    end
    return comp = Dict{Symbol,Property}((i,j) for (i,j) in zip(symbols,properties))
end
=======


text:
# cell value math functions
for f in [:sin, :cos, :tan, :cot, :abs, :exp, :log, :log10]
  @eval function $f(a::CellValue)
      return CellValue(a.domain, $f.(a.value))
  end
end

# face value math functions
for f in [:sin, :cos, :tan, :cot, :abs, :exp, :log, :log10]
  @eval function $f(a::FaceValue)
      return FaceValue(a.domain,
                  $f.(a.xvalue),
                  $f.(a.yvalue),
                  $f.(a.zvalue))
  end
end

# cell value operators
for f in [:+, :-, :*, :/, :^, :(==), :>, :(>=), :<, :(<=)]
  @eval function $f(a::CellValue, b::CellValue)
      return CellValue(a.domain, $f.(a.value, b.value))
  end
end

for f in [:+, :-, :*, :/, :^, :(==), :>, :(>=), :<, :(<=)]
  @eval function $f(a::CellValue, b::Real)
      return CellValue(a.domain, $f.(a.value, b))
  end
end

for f in [:+, :-, :*, :/, :^, :(==), :>, :(>=), :<, :(<=)]
  @eval function $f(a::Real, b::CellValue)
      return CellValue(b.domain, $f.(a, b.value))
  end
end

# face value operators
for f in [:+, :-, :*, :/, :^, :(==), :>, :(>=), :<, :(<=)]
  @eval function $f(a::FaceValue, b::FaceValue)
      return FaceValue(a.domain,
        $f.(a.xvalue, b.xvalue),
        $f.(a.yvalue, b.yvalue),
        $f.(a.zvalue, b.zvalue))
  end
end

for f in [:+, :-, :*, :/, :^, :(==), :>, :(>=), :<, :(<=)]
  @eval function $f(a::FaceValue, b::Real)
      return FaceValue(a.domain,
        $f.(a.xvalue, b),
        $f.(a.yvalue, b),
        $f.(a.zvalue, b))
  end
end

for f in [:+, :-, :*, :/, :^, :(==), :>, :(>=), :<, :(<=)]
  @eval function $f(a::Real, b::FaceValue)
      return FaceValue(b.domain,
        $f.(a, b.xvalue),
        $f.(a, b.yvalue),
        $f.(a, b.zvalue))
  end
end

function -(a::FaceValue)
  FaceValue(a.domain,
    -a.xvalue,
    -a.yvalue,
    -a.zvalue)
end

function -(a::CellValue)
  CellValue(a.domain,
    -a.value)
end

# # cell value operators
# for f in [:+, :-, :*, :/, :^, :(==), :>, :(>=), :<, :(<=)]
#   @eval function $f(a::CellValue, b::CellValue)
#       return CellValue(a.domain, $f.(a.value, b.value))
#   end
# end

# function broadcast(::typeof(*), a::FaceValue, b::FaceValue)
#   FaceValue(a.domain,
#     a.xvalue.*b.xvalue,
#     a.yvalue.*b.yvalue,
#     a.zvalue.*b.zvalue)
# end

# function broadcast(::typeof(/), a::FaceValue, b::FaceValue)
#   FaceValue(a.domain,
#     a.xvalue./b.xvalue,
#     a.yvalue./b.yvalue,
#     a.zvalue./b.zvalue)
# end

# function broadcast(::typeof(*), a::CellValue, b::CellValue)
# CellValue(a.domain,
#   a.value.*b.value)
# end

# function broadcast(::typeof(/), a::CellValue, b::CellValue)
# CellValue(a.domain,
#   a.value./b.value)
# end

# function broadcast(::typeof(*), a::CellVector, b::CellVector)
# CellVector(a.domain,
#   a.xvalue.*b.xvalue,
#   a.yvalue.*b.yvalue,
#   a.zvalue.*b.zvalue)
# end

# function broadcast(::typeof(/), a::CellVector, b::CellVector)
# CellVector(a.domain,
#   a.xvalue./b.xvalue,
#   a.yvalue./b.yvalue,
#   a.zvalue./b.zvalue)
# end

## To be thrown away

# function <(a::FaceValue, b::FaceValue)
#   FaceValue(a.domain,
#     a.xvalue.<b.xvalue,
#     a.yvalue.<b.yvalue,
#     a.zvalue.<b.zvalue)
# end

# function <=(a::FaceValue, b::FaceValue)
#   FaceValue(a.domain,
#     a.xvalue.<=b.xvalue,
#     a.yvalue.<=b.yvalue,
#     a.zvalue.<=b.zvalue)
# end

# function ==(a::FaceValue, b::FaceValue)
#   FaceValue(a.domain,
#     a.xvalue.==b.xvalue,
#     a.yvalue.==b.yvalue,
#     a.zvalue.==b.zvalue)
# end

# function >(a::FaceValue, b::FaceValue)
#   FaceValue(a.domain,
#     a.xvalue.>b.xvalue,
#     a.yvalue.>b.yvalue,
#     a.zvalue.>b.zvalue)
# end

# function >=(a::FaceValue, b::FaceValue)
#   FaceValue(a.domain,
#     a.xvalue.>=b.xvalue,
#     a.yvalue.>=b.yvalue,
#     a.zvalue.>=b.zvalue)
# end

# function +(a::FaceValue, b::FaceValue)
# FaceValue(a.domain,
#   a.xvalue+b.xvalue,
#   a.yvalue+b.yvalue,
#   a.zvalue+b.zvalue)
# end

# function *(a::FaceValue, b::FaceValue)
# FaceValue(a.domain,
#   a.xvalue.*b.xvalue,
#   a.yvalue.*b.yvalue,
#   a.zvalue.*b.zvalue)
# end

# function /(a::FaceValue, b::FaceValue)
# FaceValue(a.domain,
#   a.xvalue./b.xvalue,
#   a.yvalue./b.yvalue,
#   a.zvalue./b.zvalue)
# end

# function -(a::FaceValue, b::FaceValue)
# FaceValue(a.domain,
#   a.xvalue-b.xvalue,
#   a.yvalue-b.yvalue,
#   a.zvalue-b.zvalue)
# end

#
# function +(a::Real, b::FaceValue)
# FaceValue(b.domain,
#   a.+b.xvalue,
#   a.+b.yvalue,
#   a.+b.zvalue)
# end

# function +(a::FaceValue, b::Real)
# FaceValue(a.domain,
#   b.+a.xvalue,
#   b.+a.yvalue,
#   b.+a.zvalue)
# end

# function -(a::Real, b::FaceValue)
# FaceValue(b.domain,
#   a.-b.xvalue,
#   a.-b.yvalue,
#   a.-b.zvalue)
# end

# function -(a::FaceValue)
# FaceValue(a.domain,
#   -a.xvalue,
#   -a.yvalue,
#   -a.zvalue)
# end

# function -(a::FaceValue, b::Real)
# FaceValue(a.domain,
#   a.xvalue.-b,
#   a.yvalue.-b,
#   a.zvalue.-b)
# end

# function *(a::Real, b::FaceValue)
# FaceValue(b.domain,
#   a.*b.xvalue,
#   a.*b.yvalue,
#   a.*b.zvalue)
# end

# function *(a::FaceValue, b::Real)
# FaceValue(a.domain,
#   b.*a.xvalue,
#   b.*a.yvalue,
#   b.*a.zvalue)
# end

# function /(a::Real, b::FaceValue)
# FaceValue(b.domain,
#   a./b.xvalue,
#   a./b.yvalue,
#   a./b.zvalue)
# end

# function /(a::FaceValue, b::Real)
# FaceValue(a.domain,
#   a.xvalue./b,
#   a.yvalue./b,
#   a.zvalue./b)
# end

# # Cell Value operators
# function ==(a::CellValue, b::CellValue)
#   CellValue(a.domain,
#     a.value.==b.value)
# end

# function <(a::CellValue, b::CellValue)
#   CellValue(a.domain,
#     a.value.<b.value)
# end

# function <=(a::CellValue, b::CellValue)
#   CellValue(a.domain,
#     a.value.<=b.value)
# end

# function >(a::CellValue, b::CellValue)
#   CellValue(a.domain,
#     a.value.>b.value)
# end

# function >=(a::CellValue, b::CellValue)
#   CellValue(a.domain,
#     a.value.>=b.value)
# end

# function +(a::CellValue, b::CellValue)
# CellValue(a.domain,
#   a.value+b.value)
# end

# function *(a::CellValue, b::CellValue)
# CellValue(a.domain,
#   a.value.*b.value)
# end

# function /(a::CellValue, b::CellValue)
# CellValue(a.domain,
#   a.value./b.value)
# end

# function -(a::CellValue, b::CellValue)
# CellValue(a.domain,
#   a.value-b.value)
# end

# function +(a::Real, b::CellValue)
# CellValue(b.domain,
#   a.+b.value)
# end

# function +(a::CellValue, b::Real)
# CellValue(a.domain,
#   b.+a.value)
# end

# function -(a::Real, b::CellValue)
# CellValue(b.domain,
#   a-b.value)
# end

# function -(a::CellValue)
# CellValue(a.domain,
#   -a.value)
# end

# function -(a::CellValue, b::Real)
# CellValue(a.domain,
#   a.value-b)
# end

# function *(a::Real, b::CellValue)
# CellValue(b.domain,
#   a*b.value)
# end

# function *(a::CellValue, b::Real)
# CellValue(a.domain,
#   b*a.value)
# end

# function /(a::Real, b::CellValue)
# CellValue(b.domain,
#   a./b.value)
# end

# function /(a::CellValue, b::Real)
# CellValue(a.domain,
#   a.value/b)
# end

# # Cell Vector operators
# function +(a::CellVector, b::CellVector)
# CellVector(a.domain,
#   a.xvalue+b.xvalue,
#   a.yvalue+b.yvalue,
#   a.zvalue+b.zvalue)
# end

# function -(a::CellVector, b::CellVector)
# CellVector(a.domain,
#   a.xvalue-b.xvalue,
#   a.yvalue-b.yvalue,
#   a.zvalue-b.zvalue)
# end

# function +(a::Real, b::CellVector)
# CellVector(b.domain,
#   a.+b.xvalue,
#   a.+b.yvalue,
#   a.+b.zvalue)
# end

# function +(a::CellVector, b::Real)
# CellVector(a.domain,
#   b.+a.xvalue,
#   b.+a.yvalue,
#   b.+a.zvalue)
# end

# function -(a::Real, b::CellVector)
# CellVector(b.domain,
#   a.-b.xvalue,
#   a.-b.yvalue,
#   a.-b.zvalue)
# end

# function -(a::CellVector)
# CellVector(a.domain,
#   -a.xvalue,
#   -a.yvalue,
#   -a.zvalue)
# end

# function -(a::CellVector, b::Real)
# CellVector(a.domain,
#   a.xvalue-b,
#   a.yvalue-b,
#   a.zvalue-b)
# end

# function *(a::Real, b::CellVector)
# CellVector(b.domain,
#   a*b.xvalue,
#   a*b.yvalue,
#   a*b.zvalue)
# end

# function *(a::CellVector, b::Real)
# CellVector(a.domain,
#   b*a.xvalue,
#   b*a.yvalue,
#   b*a.zvalue)
# end

# function /(a::Real, b::CellVector)
# CellVector(b.domain,
#   a./b.xvalue,
#   a./b.yvalue,
#   a./b.zvalue)
# end

# function /(a::CellVector, b::Real)
# CellVector(a.domain,
#   a.xvalue/b,
#   a.yvalue/b,
#   a.zvalue/b)
# end


text:
methods_implemented = Dict{String,Vector{String}}()
methods_implemented["GaussianLikelihood"] = []
methods_implemented["StudentTLikelihood"] = ["AnalyticVI","AnalyticSVI","QuadratureVI","QuadratureSVI"] # ["QuadratureVI","QuadratureSVI"]
methods_implemented["LaplaceLikelihood"] = ["AnalyticVI","AnalyticSVI","QuadratureVI","QuadratureSVI"]
methods_implemented["HeteroscedasticLikelihood"] = []
methods_implemented["LogisticLikelihood"] = ["AnalyticVI","AnalyticSVI","QuadratureVI","QuadratureSVI"]
methods_implemented["BayesianSVM"] = ["AnalyticVI","AnalyticSVI"]
methods_implemented["LogisticSoftMaxLikelihood"] = ["AnalyticVI","AnalyticSVI"]# "NumericalVI","NumericalSVI"]
methods_implemented["SoftMaxLikelihood"] = ["MCIntegrationVI","MCIntegrationSVI"]
methods_implemented["PoissonLikelihood"] = ["AnalyticVI","AnalyticSVI"]

methods_implemented_VGP = deepcopy(methods_implemented)
push!(methods_implemented_VGP["StudentTLikelihood"],"GibbsSampling")
push!(methods_implemented_VGP["LaplaceLikelihood"],"GibbsSampling")
push!(methods_implemented_VGP["LogisticLikelihood"],"GibbsSampling")
push!(methods_implemented_VGP["LogisticSoftMaxLikelihood"],"GibbsSampling")
push!(methods_implemented_VGP["HeteroscedasticLikelihood"],"AnalyticVI")
methods_implemented_SVGP = deepcopy(methods_implemented)
methods_implemented_SVGP["GaussianLikelihood"] = ["AnalyticVI","AnalyticSVI"]

stoch(s::Bool,inference::String) = inference != "GibbsSampling" ? (s ? inference[1:end-2]*"SVI" : inference) : inference
addiargument(s::Bool,inference::String) = inference == "GibbsSampling" ? "nBurnin=0" : (s ? "b" : "")
addlargument(likelihood::String) = begin
    if (likelihood == "StudentTLikelihood")
        return "ν"
    else
        return ""
    end
end

function testconv(model::AbstractGP,problem_type::String,X::AbstractArray,y::AbstractArray)
    μ,Σ = predict_f(model,X,covf=true)
    y_pred = predict_y(model,X)
    py_pred = proba_y(model,X)
    if problem_type == "Regression"
        err=mean(abs.(y_pred-y))
        @info "Regression Error" err
        return err < 1.5
    elseif problem_type == "Classification"
        err=mean(y_pred.!=y)
        @info "Classification Error" err
        return err < 0.5
    elseif problem_type == "MultiClass"
        err = mean(y_pred.!=y)
        @info "Multiclass Error" err
        return err < 0.5
    elseif problem_type == "Event"
        err = mean(abs.(y_pred-y))
        @info "Event Error" err
        return err < 1.0
    end
end


text:
using ProtoBuf, Test

@testset "ProtoBuf Tests" begin
    include("testprotoc.jl")
    include("services/testsvc.jl")
    include("testmetalock.jl")
    include("testtypevers.jl")
    include("testutilapi.jl")
    include("testwellknown.jl")
    include("testcodec.jl")
end