text:
#using MathProgBase
import MathOptInterface
#using SparseArrays

########################################################
## this code is based on ModelReader in NLPModels
## and KNITRO.jl
########################################################

export OnePhaseSolver #, MathOptNLPModel, OnePhaseMathProgModel, NonlinearModel

const MOI  = MathOptInterface
const MOIU = MathOptInterface.Utilities

#const MPB = MathProgBase

# TODO
const SF = Union{MOI.ScalarAffineFunction{Float64},
                 MOI.ScalarQuadraticFunction{Float64}}
const VAF = MOI.VectorAffineFunction{Float64}
const VOV = MOI.VectorOfVariables

# ScalarAffineFunctions and VectorAffineFunctions
const SAF = MOI.ScalarAffineFunction{Float64}
const AF = Union{SAF, VAF}

const SS = Union{MOI.EqualTo{Float64},
                 MOI.GreaterThan{Float64},
                 MOI.LessThan{Float64},
                 MOI.Interval{Float64}}
# LinSets
const LS = Union{MOI.EqualTo{Float64},
                 MOI.GreaterThan{Float64},
                 MOI.LessThan{Float64}}
# VecLinSets
const VLS = Union{MOI.Nonnegatives,
                  MOI.Nonpositives,
                  MOI.Zeros}

##################################################
mutable struct VariableInfo
    lower_bound::Float64  # May be -Inf even if has_lower_bound == true
    has_lower_bound::Bool # Implies lower_bound == Inf
    lower_bound_dual_start::Union{Nothing, Float64}
    upper_bound::Float64  # May be Inf even if has_upper_bound == true
    has_upper_bound::Bool # Implies upper_bound == Inf
    upper_bound_dual_start::Union{Nothing, Float64}
    is_fixed::Bool        # Implies lower_bound == upper_bound and !has_lower_bound and !has_upper_bound.
    start::Union{Nothing, Float64}
end
VariableInfo() = VariableInfo(-Inf, false, nothing, Inf, false, nothing, false, nothing)

mutable struct OnePhaseProblem
    status::Symbol  # Final status
    #status::MOI.TerminationStatusCode
    # For MathProgBase
    x::Vector{Float64}  # Starting and final solution
    lambda::Vector{Float64}
    g::Vector{Float64}  # Final constraint values
    obj_val::Float64  # (length 1) Final objective
    solve_time::Float64

    # Custom attributes of the OnePhaseSolver
    iter::Class_iterate
    hist::Array{alg_history2,1}
    pars::Class_parameters

    function OnePhaseProblem()
        return new()
    end
end

##################################################
# EmptyNLPEvaluator for non-NLP problems.
struct EmptyNLPEvaluator <: MOI.AbstractNLPEvaluator end
MOI.features_available(::EmptyNLPEvaluator) = [:Grad, :Jac, :Hess]
MOI.initialize(::EmptyNLPEvaluator, features) = nothing
MOI.eval_objective(::EmptyNLPEvaluator, x) = NaN
function MOI.eval_constraint(::EmptyNLPEvaluator, g, x)
    @assert length(g) == 0
    return
end
function MOI.eval_objective_gradient(::EmptyNLPEvaluator, g, x)
    fill!(g, 0.0)
    return
end
MOI.jacobian_structure(::EmptyNLPEvaluator) = Tuple{Int64,Int64}[]
MOI.hessian_lagrangian_structure(::EmptyNLPEvaluator) = Tuple{Int64,Int64}[]
function MOI.eval_constraint_jacobian(::EmptyNLPEvaluator, J, x)
    @assert length(J) == 0
    return
end
function MOI.eval_hessian_lagrangian(::EmptyNLPEvaluator, H, x, s, mu)
    @assert length(H) == 0
    return
end
function MOI.eval_hessian_lagrangian(::Nothing, H, x, s, mu)
    @assert length(H) == 0
    return
end

empty_nlp_data() = MOI.NLPBlockData([], EmptyNLPEvaluator(), false)

mutable struct ConstraintInfo{F, S}
    func::F
    set::S
    dual_start::Union{Nothing, Float64}
end

ConstraintInfo(func, set) = ConstraintInfo(func, set, nothing)

mutable struct OnePhaseSolver <: MOI.AbstractOptimizer
    #inner::OnePhaseProblem
	inner::Union{OnePhaseProblem, Nothing}
	#inner::Union{Model, Nothing}

    # Problem data.
    #eval :: Union{MOI.AbstractNLPEvaluator, Nothing}
	variable_info::Vector{VariableInfo}
	nlp_data::MOI.NLPBlockData
	sense :: MOI.OptimizationSense
    objective::Union{MOI.SingleVariable,MOI.ScalarAffineFunction{Float64},MOI.ScalarQuadraticFunction{Float64},Nothing}
    linear_le_constraints::Vector{ConstraintInfo{MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}}}
    linear_ge_constraints::Vector{ConstraintInfo{MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}}}
    linear_eq_constraints::Vector{ConstraintInfo{MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}}}
	linear_int_constraints::Vector{ConstraintInfo{MOI.ScalarAffineFunction{Float64}, MOI.Interval{Float64}}}
    quadratic_le_constraints::Vector{ConstraintInfo{MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}}}
    quadratic_ge_constraints::Vector{ConstraintInfo{MOI.ScalarQuadraticFunction{Float64}, MOI.GreaterThan{Float64}}}
    quadratic_eq_constraints::Vector{ConstraintInfo{MOI.ScalarQuadraticFunction{Float64}, MOI.EqualTo{Float64}}}
	quadratic_int_constraints::Vector{ConstraintInfo{MOI.ScalarQuadraticFunction{Float64}, MOI.Interval{Float64}}}
    nlp_dual_start::Union{Nothing, Vector{Float64}}

	# Parameters.
	silent::Bool
	options::Dict{String, Any}

	# Solution attributes.
    solve_time::Float64
end

function OnePhaseSolver(; options...)
	options_dict = Dict{String, Any}()

	for (name, value) in options
        options_dict[string(name)] = value
    end

	onePhaseSolverModel = OnePhaseSolver(
        OnePhaseProblem(),
        [],
        empty_nlp_data(),
        MOI.FEASIBILITY_SENSE,
        nothing,
        [],
        [],
        [],
        [],
        [],
        [],
		[],
		[],
        nothing,
        false,
        Dict{String, Any}(),
        NaN,
    )
    set_options(onePhaseSolverModel, options)

	return onePhaseSolverModel
end

function set_options(model::OnePhaseSolver, options)
    for (name, value) in options
        sname = string(name)
        MOI.set(model, MOI.RawParameter(sname), value)
    end
    return
end

MOI.get(::OnePhaseSolver, ::MOI.SolverName) = "OnePhaseSolver"

"""
    MOI.is_empty(model::OnePhaseSolver )
"""

function MOI.is_empty(model::OnePhaseSolver)
    return isempty(model.variable_info) &&
           model.nlp_data.evaluator isa EmptyNLPEvaluator &&
           model.sense == MOI.FEASIBILITY_SENSE &&
           isempty(model.linear_le_constraints) &&
           isempty(model.linear_ge_constraints) &&
           isempty(model.linear_eq_constraints) &&
		   isempty(model.linear_int_constraints) &&
           isempty(model.quadratic_le_constraints) &&
           isempty(model.quadratic_ge_constraints) &&
           isempty(model.quadratic_eq_constraints) &&
		   isempty(model.quadratic_int_constraints)
end

MOI.get(model::OnePhaseSolver, ::MOI.SolveTime) = model.solve_time

function MOI.empty!(model::OnePhaseSolver)
    model.inner = nothing
    empty!(model.variable_info)
    model.nlp_data = empty_nlp_data()
    model.sense = MOI.FEASIBILITY_SENSE
    model.objective = nothing
    empty!(model.linear_le_constraints)
    empty!(model.linear_ge_constraints)
    empty!(model.linear_eq_constraints)
	empty!(model.linear_int_constraints)
    empty!(model.quadratic_le_constraints)
    empty!(model.quadratic_ge_constraints)
    empty!(model.quadratic_eq_constraints)
	empty!(model.quadratic_int_constraints)
    model.nlp_dual_start = nothing
end

"""
    column(x::MOI.VariableIndex)
Return the column associated with a variable.
"""
column(x::MOI.VariableIndex) = x.value

function MOI.add_variable(model::OnePhaseSolver)
    push!(model.variable_info, VariableInfo())
    return MOI.VariableIndex(length(model.variable_info))
end

function MOI.add_variables(model::OnePhaseSolver, n::Int)
    return [MOI.add_variable(model) for i in 1:n]
end

function MOI.is_valid(model::OnePhaseSolver, vi::MOI.VariableIndex)
    return column(vi) in eachindex(model.variable_info)
end

function MOI.Utilities.supports_default_copy_to(::OnePhaseSolver, copy_names::Bool)
    return !copy_names
end

function MOI.copy_to(model::OnePhaseSolver, src::MOI.ModelLike; copy_names = false)
    return MOI.Utilities.default_copy_to(model, src, copy_names)
end
# MathOptInterface constraints

##################################################
## Support constraints
MOI.supports_constraint(::OnePhaseSolver, ::Type{MOI.SingleVariable}, ::Type{MOI.LessThan{Float64}}) = true
MOI.supports_constraint(::OnePhaseSolver, ::Type{MOI.SingleVariable}, ::Type{MOI.GreaterThan{Float64}}) = true
MOI.supports_constraint(::OnePhaseSolver, ::Type{MOI.SingleVariable}, ::Type{MOI.EqualTo{Float64}}) = true
MOI.supports_constraint(::OnePhaseSolver, ::Type{MOI.ScalarAffineFunction{Float64}}, ::Type{MOI.LessThan{Float64}}) = true
MOI.supports_constraint(::OnePhaseSolver, ::Type{MOI.ScalarAffineFunction{Float64}}, ::Type{MOI.GreaterThan{Float64}}) = true
MOI.supports_constraint(::OnePhaseSolver, ::Type{MOI.ScalarAffineFunction{Float64}}, ::Type{MOI.EqualTo{Float64}}) = true
MOI.supports_constraint(::OnePhaseSolver, ::Type{MOI.ScalarAffineFunction{Float64}}, ::Type{MOI.Interval{Float64}}) = true
MOI.supports_constraint(::OnePhaseSolver, ::Type{MOI.ScalarQuadraticFunction{Float64}}, ::Type{MOI.LessThan{Float64}}) = true
MOI.supports_constraint(::OnePhaseSolver, ::Type{MOI.ScalarQuadraticFunction{Float64}}, ::Type{MOI.GreaterThan{Float64}}) = true
MOI.supports_constraint(::OnePhaseSolver, ::Type{MOI.ScalarQuadraticFunction{Float64}}, ::Type{MOI.EqualTo{Float64}}) = true
MOI.supports_constraint(::OnePhaseSolver, ::Type{MOI.ScalarQuadraticFunction{Float64}}, ::Type{MOI.Interval{Float64}}) = true

function has_upper_bound(model::OnePhaseSolver, vi::MOI.VariableIndex)
    return model.variable_info[column(vi)].has_upper_bound
end

function has_lower_bound(model::OnePhaseSolver, vi::MOI.VariableIndex)
    return model.variable_info[column(vi)].has_lower_bound
end

function is_fixed(model::OnePhaseSolver, vi::MOI.VariableIndex)
    return model.variable_info[column(vi)].is_fixed
end

function MOI.add_constraint(
    model::OnePhaseSolver, v::MOI.SingleVariable, lt::MOI.LessThan{Float64},
)
    vi = v.variable
    MOI.throw_if_not_valid(model, vi)
    if isnan(lt.upper)
        error("Invalid upper bound value $(lt.upper).")
    end
    if has_upper_bound(model, vi)
        throw(MOI.UpperBoundAlreadySet{typeof(lt), typeof(lt)}(vi))
    end
    if is_fixed(model, vi)
        throw(MOI.UpperBoundAlreadySet{MOI.EqualTo{Float64}, typeof(lt)}(vi))
    end
    col = column(vi)
    model.variable_info[col].upper_bound = lt.upper
    model.variable_info[col].has_upper_bound = true
    return MOI.ConstraintIndex{MOI.SingleVariable, MOI.LessThan{Float64}}(col)
end

function MOI.set(
    model::OnePhaseSolver,
    ::MOI.ConstraintSet,
    ci::MOI.ConstraintIndex{MOI.SingleVariable, MOI.LessThan{Float64}},
    set::MOI.LessThan{Float64},
)
    MOI.throw_if_not_valid(model, ci)
    model.variable_info[ci.value].upper_bound = set.upper
    return
end

function MOI.delete(
    model::OnePhaseSolver,
    ci::MOI.ConstraintIndex{MOI.SingleVariable, MOI.LessThan{Float64}},
)
    MOI.throw_if_not_valid(model, ci)
    model.variable_info[ci.value].upper_bound = Inf
    model.variable_info[ci.value].has_upper_bound = false
    return
end

function MOI.add_constraint(
    model::OnePhaseSolver, v::MOI.SingleVariable, gt::MOI.GreaterThan{Float64},
)
    vi = v.variable
    MOI.throw_if_not_valid(model, vi)
    if isnan(gt.lower)
        error("Invalid lower bound value $(gt.lower).")
    end
    if has_lower_bound(model, vi)
        throw(MOI.LowerBoundAlreadySet{typeof(gt), typeof(gt)}(vi))
    end
    if is_fixed(model, vi)
        throw(MOI.LowerBoundAlreadySet{MOI.EqualTo{Float64}, typeof(gt)}(vi))
    end
    col = column(vi)
    model.variable_info[col].lower_bound = gt.lower
    model.variable_info[col].has_lower_bound = true
    return MOI.ConstraintIndex{MOI.SingleVariable, MOI.GreaterThan{Float64}}(col)
end

function MOI.set(
    model::OnePhaseSolver,
    ::MOI.ConstraintSet,
    ci::MOI.ConstraintIndex{MOI.SingleVariable, MOI.GreaterThan{Float64}},
    set::MOI.GreaterThan{Float64},
)
    MOI.throw_if_not_valid(model, ci)
    model.variable_info[ci.value].lower_bound = set.lower
    return
end

function MOI.delete(
    model::OnePhaseSolver,
    ci::MOI.ConstraintIndex{MOI.SingleVariable, MOI.GreaterThan{Float64}},
)
    MOI.throw_if_not_valid(model, ci)
    model.variable_info[ci.value].lower_bound = -Inf
    model.variable_info[ci.value].has_lower_bound = false
    return
end

function MOI.add_constraint(
    model::OnePhaseSolver, v::MOI.SingleVariable, eq::MOI.EqualTo{Float64},
)
    vi = v.variable
    MOI.throw_if_not_valid(model, vi)
    if isnan(eq.value)
        error("Invalid fixed value $(eq.value).")
    end
    if has_lower_bound(model, vi)
        throw(MOI.LowerBoundAlreadySet{MOI.GreaterThan{Float64}, typeof(eq)}(vi))
    end
    if has_upper_bound(model, vi)
        throw(MOI.UpperBoundAlreadySet{MOI.LessThan{Float64}, typeof(eq)}(vi))
    end
    if is_fixed(model, vi)
        throw(MOI.LowerBoundAlreadySet{typeof(eq), typeof(eq)}(vi))
    end
    col = column(vi)
    model.variable_info[col].lower_bound = eq.value
    model.variable_info[col].upper_bound = eq.value
    model.variable_info[col].is_fixed = true
    return MOI.ConstraintIndex{MOI.SingleVariable, MOI.EqualTo{Float64}}(col)
end

function MOI.set(
    model::OnePhaseSolver,
    ::MOI.ConstraintSet,
    ci::MOI.ConstraintIndex{MOI.SingleVariable, MOI.EqualTo{Float64}},
    set::MOI.EqualTo{Float64},
)
    MOI.throw_if_not_valid(model, ci)
    model.variable_info[ci.value].lower_bound = set.value
    model.variable_info[ci.value].upper_bound = set.value
    return
end

function MOI.delete(
    model::OnePhaseSolver,
    ci::MOI.ConstraintIndex{MOI.SingleVariable, MOI.EqualTo{Float64}},
)
    MOI.throw_if_not_valid(model, ci)
    model.variable_info[ci.value].lower_bound = -Inf
    model.variable_info[ci.value].upper_bound = Inf
    model.variable_info[ci.value].is_fixed = false
    return
end

macro define_add_constraint(function_type, set_type, prefix)
    array_name = Symbol(string(prefix) * "_constraints")
    return quote
        function MOI.add_constraint(
            model::OnePhaseSolver, func::$function_type, set::$set_type,
        )
            check_inbounds(model, func)
            push!(model.$(array_name), ConstraintInfo(func, set))
            return MOI.ConstraintIndex{$function_type, $set_type}(length(model.$(array_name)))
        end
    end
end

@define_add_constraint(
    MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}, linear_le,
)

@define_add_constraint(
    MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}, linear_ge,
)

@define_add_constraint(
    MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}, linear_eq,
)

@define_add_constraint(
    MOI.ScalarAffineFunction{Float64}, MOI.Interval{Float64}, linear_int,
)

@define_add_constraint(
    MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}, quadratic_le,
)

@define_add_constraint(
    MOI.ScalarQuadraticFunction{Float64}, MOI.GreaterThan{Float64}, quadratic_ge,
)

@define_add_constraint(
    MOI.ScalarQuadraticFunction{Float64}, MOI.EqualTo{Float64}, quadratic_eq,
)

@define_add_constraint(
    MOI.ScalarQuadraticFunction{Float64}, MOI.Interval{Float64}, quadratic_int,
)

function MOI.set(
    model::OnePhaseSolver,
    ::MOI.ObjectiveFunction,
    func::Union{
        MOI.SingleVariable,
        MOI.ScalarAffineFunction,
        MOI.ScalarQuadraticFunction,
    },
)
    check_inbounds(model, func)
    model.objective = func
    return
end

#Utilities Begin Here
"""
    parser_SAF(fun, set, linrows, lincols, linvals, nlin, lin_lcon, lin_ucon)
Parse a `ScalarAffineFunction` fun with its associated set.
`linrows`, `lincols`, `linvals`, `lin_lcon` and `lin_ucon` are updated.
"""
function parser_SAF(fun, set, linrows, lincols, linvals, nlin, lin_lcon, lin_ucon)

  # Parse a ScalarAffineTerm{Float64}(coefficient, variable_index)
  for term in fun.terms
    push!(linrows, nlin + 1)
    push!(lincols, term.variable_index.value)
    push!(linvals, term.coefficient)
  end

  if typeof(set) in (MOI.Interval{Float64}, MOI.GreaterThan{Float64})
    push!(lin_lcon, -fun.constant + set.lower)
  elseif typeof(set) == MOI.EqualTo{Float64}
    push!(lin_lcon, -fun.constant + set.value)
  else
    push!(lin_lcon, -Inf)
  end

  if typeof(set) in (MOI.Interval{Float64}, MOI.LessThan{Float64})
    push!(lin_ucon, -fun.constant + set.upper)
  elseif typeof(set) == MOI.EqualTo{Float64}
    push!(lin_ucon, -fun.constant + set.value)
  else
    push!(lin_ucon, Inf)
  end
end

"""
    parser_VAF(fun, set, linrows, lincols, linvals, nlin, lin_lcon, lin_ucon)
Parse a `VectorAffineFunction` fun with its associated set.
`linrows`, `lincols`, `linvals`, `lin_lcon` and `lin_ucon` are updated.
"""
function parser_VAF(fun, set, linrows, lincols, linvals, nlin, lin_lcon, lin_ucon)

  # Parse a VectorAffineTerm{Float64}(output_index, scalar_term)
  for term in fun.terms
    push!(linrows, nlin + term.output_index)
    push!(lincols, term.scalar_term.variable_index.value)
    push!(linvals, term.scalar_term.coefficient)
  end

  if typeof(set) in (MOI.Nonnegatives, MOI.Zeros)
    append!(lin_lcon, -fun.constants)
  else
    append!(lin_lcon, -Inf * ones(set.dimension))
  end

  if typeof(set) in (MOI.Nonpositives, MOI.Zeros)
    append!(lin_ucon, -fun.constants)
  else
    append!(lin_ucon, Inf * ones(set.dimension))
  end
end

"""
    parser_MOI(moimodel)
Parse linear constraints of a `MOI.ModelLike`.
"""
function parser_MOI(moimodel)
  # Variables associated to linear constraints
  nlin = 0
  linrows = Int[]
  lincols = Int[]
  linvals = Float64[]
  lin_lcon = Float64[]
  lin_ucon = Float64[]

  contypes = MOI.get(moimodel, MOI.ListOfConstraints())
  for (F, S) in contypes
    F == MOI.SingleVariable && continue
	F <: AF
	S <: LS

    conindices = MOI.get(moimodel, MOI.ListOfConstraintIndices{F, S}())
    for cidx in conindices
      fun = MOI.get(moimodel, MOI.ConstraintFunction(), cidx)
      set = MOI.get(moimodel, MOI.ConstraintSet(), cidx)
      if typeof(fun) <: SAF
        parser_SAF(fun, set, linrows, lincols, linvals, nlin, lin_lcon, lin_ucon)
        nlin += 1
      end
      if typeof(fun) <: VAF
        parser_VAF(fun, set, linrows, lincols, linvals, nlin, lin_lcon, lin_ucon)
        nlin += set.dimension
      end
    end
  end
  coo = COO(linrows, lincols, linvals)
  nnzj = length(linvals)
  lincon = LinearConstraints(coo, nnzj)

  return nlin, lincon, lin_lcon, lin_ucon
end

function NLPModelsJuMP.parser_MOI(moimodel)
  # Variables associated to linear constraints
  nlin = 0
  linrows = Int[]
  lincols = Int[]
  linvals = Float64[]
  lin_lcon = Float64[]
  lin_ucon = Float64[]

  contypes = MOI.get(moimodel, MOI.ListOfConstraints())
  for (F, S) in contypes
    F == MOI.SingleVariable && continue
	F <: AF
	S <: LS

    conindices = MOI.get(moimodel, MOI.ListOfConstraintIndices{F, S}())
    for cidx in conindices
      fun = MOI.get(moimodel, MOI.ConstraintFunction(), cidx)
      set = MOI.get(moimodel, MOI.ConstraintSet(), cidx)
      if typeof(fun) <: SAF
        parser_SAF(fun, set, linrows, lincols, linvals, nlin, lin_lcon, lin_ucon)
        nlin += 1
      end
      if typeof(fun) <: VAF
        parser_VAF(fun, set, linrows, lincols, linvals, nlin, lin_lcon, lin_ucon)
        nlin += set.dimension
      end
    end
  end
  coo = NLPModelsJuMP.COO(linrows, lincols, linvals)
  nnzj = length(linvals)
  lincon = NLPModelsJuMP.LinearConstraints(coo, nnzj)

  return nlin, lincon, lin_lcon, lin_ucon
end

"""
    parser_JuMP(jmodel)
Parse variables informations of a `JuMP.Model`.
"""
function parser_JuMP(jmodel)

  # Number of variables and bounds constraints
  nvar = Int(num_variables(jmodel))
  vars = all_variables(jmodel)
  lvar = map(
    var -> JuMP.is_fixed(var) ? fix_value(var) : (JuMP.has_lower_bound(var) ? lower_bound(var) : -Inf),
    vars,
  )
  uvar = map(
    var -> JuMP.is_fixed(var) ? fix_value(var) : (JuMP.has_upper_bound(var) ? upper_bound(var) : Inf),
    vars,
  )

  # Initial solution
  x0 = zeros(nvar)
  for (i, val) in enumerate(start_value.(vars))
    if val !== nothing
      x0[i] = val
    end
  end

  return nvar, lvar, uvar, x0
end

"""
    parser_objective_MOI(moimodel, nvar)
Parse linear and quadratic objective of a `MOI.ModelLike`.
"""
function parser_objective_MOI(moimodel, nvar)

  # Variables associated to linear and quadratic objective
  type = "UNKNOWN"
  constant = 0.0
  vect = spzeros(Float64, nvar)
  rows = Int[]
  cols = Int[]
  vals = Float64[]

  fobj = MOI.get(moimodel, MOI.ObjectiveFunction{OBJ}())

  # Single Variable
  if typeof(fobj) == SV
    type = "LINEAR"
    vect[fobj.variable.value] = 1.0
  end

  # Linear objective
  if typeof(fobj) == SAF
    type = "LINEAR"
    constant = fobj.constant
    for term in fobj.terms
      vect[term.variable_index.value] = term.coefficient
    end
  end

  # Quadratic objective
  if typeof(fobj) == SQF
    type = "QUADRATIC"
    constant = fobj.constant
    for term in fobj.affine_terms
      vect[term.variable_index.value] = term.coefficient
    end
    for term in fobj.quadratic_terms
      i = term.variable_index_1.value
      j = term.variable_index_2.value
      if i >= j
        push!(rows, i)
        push!(cols, j)
      else
        push!(cols, j)
        push!(rows, i)
      end
      push!(vals, term.coefficient)
    end
  end
  return Objective(type, constant, vect, COO(rows, cols, vals), length(vals))
end

"""
    parser_linear_expression(cmodel, nvar, F)
Parse linear expressions of type `GenericAffExpr{Float64,VariableRef}`.
"""
function parser_linear_expression(cmodel, nvar, F)

  # Variables associated to linear expressions
  rows = Int[]
  cols = Int[]
  vals = Float64[]
  constants = Float64[]

  # Linear least squares model
  nlinequ = 0
  F_is_array_of_containers = F isa Array{<:AbstractArray}
  if F_is_array_of_containers
    @objective(
      cmodel,
      Min,
      0.0 +
      0.5 *
      sum(sum(Fi^2 for Fi in FF if typeof(Fi) == GenericAffExpr{Float64, VariableRef}) for FF in F)
    )
    for FF in F, expr in FF
      if typeof(expr) == GenericAffExpr{Float64, VariableRef}
        nlinequ += 1
        for (i, key) in enumerate(expr.terms.keys)
          push!(rows, nlinequ)
          push!(cols, key.index.value)
          push!(vals, expr.terms.vals[i])
        end
        push!(constants, expr.constant)
      end
    end
  else
    @objective(
      cmodel,
      Min,
      0.0 + 0.5 * sum(Fi^2 for Fi in F if typeof(Fi) == GenericAffExpr{Float64, VariableRef})
    )
    for expr in F
      if typeof(expr) == GenericAffExpr{Float64, VariableRef}
        nlinequ += 1
        for (i, key) in enumerate(expr.terms.keys)
          push!(rows, nlinequ)
          push!(cols, key.index.value)
          push!(vals, expr.terms.vals[i])
        end
        push!(constants, expr.constant)
      end
    end
  end
  moimodel = backend(cmodel)
  lls = parser_objective_MOI(moimodel, nvar)
  return lls, LinearEquations(COO(rows, cols, vals), constants, length(vals)), nlinequ
end

"""
    parser_nonlinear_expression(cmodel, nvar, F)
Parse nonlinear expressions of type `NonlinearExpression`.
"""
function parser_nonlinear_expression(cmodel, nvar, F)

  # Nonlinear least squares model
  nnlnequ = 0
  F_is_array_of_containers = F isa Array{<:AbstractArray}
  if F_is_array_of_containers
    nnlnequ = sum(sum(typeof(Fi) == NonlinearExpression for Fi in FF) for FF in F)
    if nnlnequ > 0
      @NLobjective(
        cmodel,
        Min,
        0.5 * sum(sum(Fi^2 for Fi in FF if typeof(Fi) == NonlinearExpression) for FF in F)
      )
    end
  else
    nnlnequ = sum(typeof(Fi) == NonlinearExpression for Fi in F)
    if nnlnequ > 0
      @NLobjective(cmodel, Min, 0.5 * sum(Fi^2 for Fi in F if typeof(Fi) == NonlinearExpression))
    end
  end
  ceval = cmodel.nlp_data == nothing ? nothing : NLPEvaluator(cmodel)
  (ceval != nothing) && (nnlnequ == 0) && MOI.initialize(ceval, [:Grad, :Jac, :Hess, :HessVec])  # Add :JacVec when available
  (ceval != nothing) &&
    (nnlnequ > 0) &&
    MOI.initialize(ceval, [:Grad, :Jac, :Hess, :HessVec, :ExprGraph])  # Add :JacVec when available

  if nnlnequ == 0
    Feval = nothing
  else
    Fmodel = JuMP.Model()
    @variable(Fmodel, x[1:nvar])
    JuMP._init_NLP(Fmodel)
    @objective(Fmodel, Min, 0.0)
    Fmodel.nlp_data.user_operators = cmodel.nlp_data.user_operators
    if F_is_array_of_containers
      for FF in F, Fi in FF
        if typeof(Fi) == NonlinearExpression
          expr = ceval.subexpressions_as_julia_expressions[Fi.index]
          replace!(expr, x)
          expr = :($expr == 0)
          JuMP.add_NL_constraint(Fmodel, expr)
        end
      end
    else
      for Fi in F
        if typeof(Fi) == NonlinearExpression
          expr = ceval.subexpressions_as_julia_expressions[Fi.index]
          replace!(expr, x)
          expr = :($expr == 0)
          JuMP.add_NL_constraint(Fmodel, expr)
        end
      end
    end
    Feval = NLPEvaluator(Fmodel)
    MOI.initialize(Feval, [:Grad, :Jac, :Hess, :HessVec])  # Add :JacVec when available
  end
  return ceval, Feval, nnlnequ
end

#Utilities End Here

function NLPModels.obj(nlp::MathOptNLPModel, x::AbstractVector)
  NLPModels.increment!(nlp, :neval_obj)
  if nlp.obj.type == "LINEAR"
    res = dot(nlp.obj.gradient, x) + nlp.obj.constant
  end
  if nlp.obj.type == "QUADRATIC"
    res =
      0.5 * coo_sym_dot(nlp.obj.hessian.rows, nlp.obj.hessian.cols, nlp.obj.hessian.vals, x, x) +
      dot(nlp.obj.gradient, x) +
      nlp.obj.constant
  end
  if nlp.obj.type == "NONLINEAR"
    res = MOI.eval_objective(nlp.eval, x)
  end
  return res
end

function append_to_hessian_sparsity!(
    ::Any,
    ::Union{MOI.SingleVariable,MOI.ScalarAffineFunction},
)
    return nothing
end

function append_to_hessian_sparsity!(
    hessian_sparsity,
    quad::MOI.ScalarQuadraticFunction,
)
    for term in quad.quadratic_terms
        push!(
            hessian_sparsity,
            (term.variable_index_1.value, term.variable_index_2.value),
        )
    end
end

function hessian_lagrangian_structure(model::OnePhaseSolver, nlp :: MathOptNLPModel)
    hessian_sparsity = Tuple{Int64,Int64}[]
    if model.objective !== nothing
        append_to_hessian_sparsity!(hessian_sparsity, model.obj)
    end
    for info in model.quadratic_le_constraints
        append_to_hessian_sparsity!(hessian_sparsity, info.func)
    end
    for info in model.quadratic_ge_constraints
        append_to_hessian_sparsity!(hessian_sparsity, info.func)
    end
    for info in model.quadratic_eq_constraints
        append_to_hessian_sparsity!(hessian_sparsity, info.func)
    end
    nlp_hessian_sparsity =
        MOI.hessian_lagrangian_structure(nlp.eval)
    append!(hessian_sparsity, nlp_hessian_sparsity)
    return hessian_sparsity
end

function hess_coord(nlp :: MathOptNLPModel, x :: Array{Float64};
  obj_weight :: Float64=1.0, y :: Array{Float64}=zeros(nlp.meta.ncon))
  NLPModels.increment!(nlp, :neval_hess)
  MOI.eval_hessian_lagrangian(nlp.eval, nlp.obj.hessian.vals, x, obj_weight, y)

  return (NLPModels.hess_structure(nlp)[1], NLPModels.hess_structure(nlp)[2], NLPModels.hess_coord(nlp, x, y, obj_weight=obj_weight))
  #return (nlp.obj.hessian.rows, nlp.obj.hessian.cols, nlp.obj.hessian.vals)
end

function NLPModels.hess_structure!(
  nlp::MathOptNLPModel,
  rows::AbstractVector{<:Integer},
  cols::AbstractVector{<:Integer},
)
  if nlp.obj.type == "QUADRATIC"
    for index = 1:(nlp.obj.nnzh)
      rows[index] = nlp.obj.hessian.rows[index]
      cols[index] = nlp.obj.hessian.cols[index]
    end
  end
  if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > 0)
    hesslag_struct = MOI.hessian_lagrangian_structure(nlp.eval)
    for index = (nlp.obj.nnzh + 1):(nlp.meta.nnzh)
      shift_index = index - nlp.obj.nnzh
      rows[index] = hesslag_struct[shift_index][1]
      cols[index] = hesslag_struct[shift_index][2]
    end
  end
  return rows, cols
end

############################
## END ModelReader CODE
############################

function status_One_Phase_To_JuMP(status::Symbol)
    # since our status are not equal to JuMPs we need to do a conversion
    if status == :Optimal
        return :Optimal
    elseif status == :primal_infeasible
        return :Infeasible
    elseif status == :dual_infeasible
        return :Unbounded
    elseif status == :MAX_IT || status === :MAX_TIME
        return :UserLimit
    else
        return :Error
    end
end

function create_pars_JuMP(options )
    pars = Class_parameters()
    for (param,value) in options
      what = split(String(param),"!") # we represent a parameter such as init.mu_scale as init!mu_scale because we cannot pass init.mu_scale as a parameter
      node = pars # root
      for i in 1:length(what)
          field = what[i]
          if i < length(what)
              node = getfield(node,Symbol(field))
          else # at the leaf
              setfield!(node,Symbol(field),value)
          end
      end
    end

    return pars
end

MOI.get(model::OnePhaseSolver, ::MOI.RawSolver) = model

function JuMP.optimize!(
    model::Model,
    # TODO: Remove the optimizer_factory and bridge_constraints
    # arguments when the deprecation error below is removed.
    optimizer_factory = nothing;
    bridge_constraints::Bool = true,
    ignore_optimize_hook = (model.optimize_hook === nothing),
    kwargs...,
)
    # The nlp_data is not kept in sync, so re-set it here.
    # TODO: Consider how to handle incremental solves.
    if model.nlp_data !== nothing
        MOI.set(model, MOI.NLPBlock(), JuMP._create_nlp_block_data(model))
        empty!(model.nlp_data.nlconstr_duals)
    end

    if optimizer_factory !== nothing
        # This argument was deprecated in JuMP 0.21.
        error(
            "The optimizer factory argument is no longer accepted by " *
            "`optimize!`. Call `set_optimizer` before `optimize!`.",
        )
    end

    # If the user or an extension has provided an optimize hook, call
    # that instead of solving the model ourselves
    if !ignore_optimize_hook
        return model.optimize_hook(model; kwargs...)
    end

    isempty(kwargs) || error(
        "Unrecognized keyword arguments: $(join([k[1] for k in kwargs], ", "))",
    )

    if mode(model) != DIRECT && MOIU.state(backend(model)) == MOIU.NO_OPTIMIZER
        throw(NoOptimizer())
    end

    try
	    m = backend(model)
		if m.mode == MathOptInterface.Utilities.AUTOMATIC && m.state == MathOptInterface.Utilities.EMPTY_OPTIMIZER
			MOIU.attach_optimizer(m)
		end

		# TODO: better error message if no optimizer is set
		@assert m.state == MathOptInterface.Utilities.ATTACHED_OPTIMIZER
        solver = m.optimizer.model
        t = time()
        nlp = MathOptNLPModel(model)

        pars = create_pars_JuMP(solver.options)

        iter, status, hist, t, err, timer = one_phase_solve(nlp,pars)

		solver.inner = OnePhaseProblem()

		solver.inner.status = status_One_Phase_To_JuMP(status)
		solver.inner.x = get_original_x(iter)
		solver.inner.obj_val = iter.cache.fval
		solver.inner.lambda = get_y(iter)
		solver.inner.solve_time = time() - t

		# custom one phase features
		solver.inner.pars = pars
		solver.inner.iter = iter
		solver.inner.hist = hist
    catch err
        # TODO: This error also be thrown also in MOI.set() if the solver is
        # attached. Currently we catch only the more common case. More generally
        # JuMP is missing a translation layer from MOI errors to JuMP errors.
        if err isa MOI.UnsupportedAttribute{MOI.NLPBlock}
            error(
                "The solver does not support nonlinear problems " *
                "(i.e., NLobjective and NLconstraint).",
            )
        else
            rethrow(err)
        end
    end

    return
end

function MOI.optimize!(solver :: OnePhaseSolver, jumpModel:: Model)
    t = time()
    nlp = MathOptNLPModel(jumpModel)
    pars = create_pars_JuMP(solver.options)

    iter, status, hist, t, err, timer = one_phase_solve(nlp,pars)
    solver.inner = OnePhaseProblem()
    solver.inner.status = status_One_Phase_To_JuMP(status)
    solver.inner.x = get_original_x(iter)
    solver.inner.obj_val = iter.cache.fval
    solver.inner.lambda = get_y(iter)
    solver.inner.solve_time = time() - t

    # custom one phase features
    solver.inner.pars = pars
    solver.inner.iter = iter
    solver.inner.hist = hist
end

function MOI.get(
    model::MathOptInterface.Utilities.CachingOptimizer,
    attr::MOI.TerminationStatus,
)
    return MOI.get(model.optimizer, attr)
end

function check_inbounds(model::OnePhaseSolver, var::MOI.SingleVariable)
    return MOI.throw_if_not_valid(model, var.variable)
end

function check_inbounds(model::OnePhaseSolver, aff::MOI.ScalarAffineFunction)
    for term in aff.terms
        MOI.throw_if_not_valid(model, term.variable_index)
    end
end

function check_inbounds(model::OnePhaseSolver, quad::MOI.ScalarQuadraticFunction)
    for term in quad.affine_terms
        MOI.throw_if_not_valid(model, term.variable_index)
    end
    for term in quad.quadratic_terms
        MOI.throw_if_not_valid(model, term.variable_index_1)
        MOI.throw_if_not_valid(model, term.variable_index_2)
    end
end

MOI.supports(::OnePhaseSolver, ::MOI.NLPBlock) = true

function MOI.supports(
    ::OnePhaseSolver, ::MOI.ObjectiveFunction{MOI.SingleVariable}
)
    return true
end

function MOI.supports(
    ::OnePhaseSolver, ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}
)
    return true
end

function MOI.supports(
    ::OnePhaseSolver, ::MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}
)
    return true
end

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

MOI.supports(::OnePhaseSolver, ::MOI.Silent) = true

MOI.supports(::OnePhaseSolver, ::MOI.RawParameter) = true

function MOI.get(model::OnePhaseSolver, ::MOI.ObjectiveFunction)
    return model.objective
end

function MOI.set(model::OnePhaseSolver, ::MOI.NLPBlock, nlp_data::MOI.NLPBlockData)
    model.nlp_data = nlp_data
    return
end

function MOI.set(
    model::OnePhaseSolver, ::MOI.ObjectiveSense, sense::MOI.OptimizationSense
)
    model.sense = sense
    return
end

MOI.get(model::OnePhaseSolver, ::MOI.ObjectiveSense) = model.sense

function MOI.set(model::OnePhaseSolver, ::MOI.Silent, value)
    model.silent = value
    return
end

MOI.get(model::OnePhaseSolver, ::MOI.Silent) = model.silent

function MOI.supports(
    ::OnePhaseSolver, ::MOI.VariablePrimalStart, ::Type{MOI.VariableIndex}
)
    return true
end

function MOI.set(
    model::OnePhaseSolver,
    ::MOI.VariablePrimalStart,
    vi::MOI.VariableIndex,
    value::Union{Real, Nothing},
)
    MOI.throw_if_not_valid(model, vi)
    model.variable_info[column(vi)].start = value
    return
end

function MOI.get(
    model::OnePhaseSolver, attr::MOI.VariablePrimal, vi::MOI.VariableIndex,
)
    MOI.check_result_index_bounds(model, attr)
    MOI.throw_if_not_valid(model, vi)
    return model.inner.x[column(vi)]
end

function MOI.set(model::OnePhaseSolver, ::MOI.TimeLimitSec, value::Real)
    MOI.set(model, MOI.RawParameter(TIME_LIMIT), Float64(value))
end

function MOI.set(model::OnePhaseSolver, p::MOI.RawParameter, value)
    model.options[p.name] = value
    return
end

function MOI.get(model::OnePhaseSolver, p::MOI.RawParameter)
    if haskey(model.options, p.name)
        return model.options[p.name]
    end
    error("RawParameter with name $(p.name) is not set.")
end

function MOI.get(model::OnePhaseSolver, ::MOI.TerminationStatus)
    if model.inner === nothing
        return MOI.OPTIMIZE_NOT_CALLED
    end
	status = model.inner.status
	return status
end

function MOI.get(model::OnePhaseSolver, ::MOI.RawStatusString)
	return string(model.inner.status)
end

# Ipopt always has an iterate available.
function MOI.get(model::OnePhaseSolver, ::MOI.ResultCount)
    return (model.inner !== nothing) ? 1 : 0
end

function MOI.get(model::OnePhaseSolver, attr::MOI.PrimalStatus)
    if !(1 <= attr.N <= MOI.get(model, MOI.ResultCount()))
        return MOI.NO_SOLUTION
    end

	status = model.inner.status
	return status
end

function MOI.get(model::OnePhaseSolver, attr::MOI.DualStatus)
    if !(1 <= attr.N <= MOI.get(model, MOI.ResultCount()))
        return MOI.NO_SOLUTION
    end
	status = model.inner.status
	return status
end

function MOI.get(model::OnePhaseSolver, attr::MOI.ObjectiveValue)
    MOI.check_result_index_bounds(model, attr)
    return model.inner.obj_val
end


text:
using ImageBasedDLForpTX
using Test

@testset "ImageBasedDLForpTX.jl" begin
    # Write your own tests here.
end

ext: jl


text:
# Algorithm 5.1
# Euler's method

function euler(f, a, b, N, α)
	result = Tuple{Float64, Float64}[]
	h = (b - a)/N
	t = α
	ω = α
	push!(result, (t, ω))
	for i ∈ 1:N
		ω = ω + h*f(t, ω)
		t = a + i*h
		push!(result, (t, ω))
	end
	return result
end

f(t, y) = y -t^2 + 1
println(euler(f, 0, 2, 10, 0.5))


text:
using Dice
using Dice: ifelse

code = @dice begin
    b = DistInt([DistBool(dicecontext(), false)])
    ifelse(flip(0.5), b, b)
end

bdd = compile(code)
@assert length(bdd.bits) > 0
infer(code, :bdd)


text:
# MIT License
# Copyright (c) 2017: Xavier Gandibleux, Anthony Przybylski, Gauthier Soleilhac, and contributors.
function solve_lexico(m::JuMP.Model, optimizer_factory, verbose; args...)
    #Retrieve objectives and their senses from vOptData
    vd = getvOptData(m)
    empty!(vd.Y_N) ; empty!(vd.X_E)
    objs = vd.objs
    objSenses = vd.objSenses
    nbObj = length(objs)

    #Check that problem is feasible and that no objective is unbounded
    for i = 1:nbObj
        JuMP.set_objective(m, objSenses[i], objs[i])
        JuMP.optimize!(m, i == 1 ? optimizer_factory : nothing, ignore_optimize_hook=true)
        status = JuMP.termination_status(m)
        if status != MOI.OPTIMAL
            return status
        end
    end

    #Create the constraints on the objectives
    cstr_rhs = [JuMP.@variable(m) for _ = 1:nbObj]
    cstr_obj = [objSenses[i] == MOI.MAX_SENSE ? JuMP.@constraint(m, objs[i] >= cstr_rhs[i]) : JuMP.@constraint(m, objs[i] <= cstr_rhs[i]) for i=1:nbObj]

    for p in Combinatorics.permutations(1:nbObj, nbObj)
        verbose && println("solving for objectives $p")
        solve_permutation(m, p, cstr_obj, cstr_rhs ; args...)
    end

    JuMP.delete.(m, cstr_obj)
    JuMP.delete.(m, cstr_rhs)

    return MOI.OPTIMAL
end

function solve_permutation(m::JuMP.Model, p, cstr_obj, cstr_rhs ; args...)

    vd = getvOptData(m)
    objs = vd.objs
    objSenses = vd.objSenses

    #Set the first objective of the permutation as an objective in the JuMP JuMP.Model
    JuMP.set_objective(m, objSenses[first(p)], objs[first(p)])

    #Solve with that objective
    JuMP.optimize!(m, ignore_optimize_hook=true)

    status = JuMP.termination_status(m)
    if status != MOI.OPTIMAL
        return status
    end

    for i = 2:length(p)
        fVal = JuMP.value(objs[p[i-1]]) #get the value for the last objective solved
        slack = objSenses[p[i-1]] == MOI.MAX_SENSE ? -1e-8 : 1e-8
        # JuMP.fix(cstr_rhs[p[i-1]], fVal - objs[p[i-1]].constant + slack)
        JuMP.fix(cstr_rhs[p[i-1]], fVal + slack)
        JuMP.set_objective(m, objSenses[p[i]], objs[p[i]]) #set the i-th objective of the permutation in the JuMP JuMP.Model
        JuMP.optimize!(m, ignore_optimize_hook=true) #and solve
    end

    varArray = JuMP.all_variables(m)
    #Store results in vOptData
    push!(vd.Y_N, map(JuMP.value, objs))
    push!(vd.X_E, JuMP.value.(varArray))

    for var in cstr_rhs
        JuMP.is_fixed(var) && JuMP.unfix.(var) #errors if calling unfix on an unfixed variable
    end

end

function solve_eps(m::JuMP.Model, optimizer_factory, ϵ::Float64, round_results, verbose, time_in_secs; args...)
    #Retrieve objectives and their senses from vOptData
    vd = getvOptData(m)
    empty!(vd.Y_N) ; empty!(vd.X_E)
    f1,f2 = vd.objs
    f1Sense,f2Sense = vd.objSenses
    varArray = JuMP.all_variables(m)

    #Set the first objective as an objective in the JuMP Model
    JuMP.set_objective(m, f1Sense, f1)

    R1 = f1Sense==MOI.MIN_SENSE ? (<=) : (>=)
    R2 = f2Sense==MOI.MIN_SENSE ? (<=) : (>=)
    weak_dom(a, b) = R1(a[1], b[1]) && R2(a[2], b[2])

    #Declare the epsilon-constraint (RHS will be set later)
    RHS = JuMP.@variable(m)
    eps = (f2Sense == MOI.MIN_SENSE) ? JuMP.@constraint(m, f2 <= RHS) : JuMP.@constraint(m, f2 >= RHS)

    #Solve with that objective
    JuMP.optimize!(m, optimizer_factory, ignore_optimize_hook=true)
    status = JuMP.termination_status(m)

    start = time()

    if status == MOI.OPTIMAL
        #While a solution exists
        while status == MOI.OPTIMAL

            # check if past time limit, then return
            if (time() - start > time_in_secs)
                #Sort X_E and Y_N
                s = sortperm(vd.Y_N, by = first)
                vd.X_E, vd.Y_N = vd.X_E[s], vd.Y_N[s]

                #Remove the epsilon constraint
                JuMP.delete(m, eps) ; JuMP.delete(m, RHS)
                return MOI.OPTIMAL
            end

            #Get the score on the objectives
            f1Val = JuMP.value(f1)
            f2Val = JuMP.value(f2)

            #If last solution found is dominated by this one
            if length(vd.Y_N) > 0
                if weak_dom((f1Val, f2Val), vd.Y_N[end])
                    verbose && println(vd.Y_N[end], "dominated by ($f1Val, $f2Val)")
                    pop!(vd.Y_N) ; pop!(vd.X_E) #Remove last solution from Y_N and X_E
                end
            end

            #Store results in vOptData
            push!(vd.X_E, JuMP.value.(varArray))
            push!(vd.Y_N, [f1Val, f2Val])

            verbose && print("z1 = ", f1Val, ", z2 = ", f2Val)

            #Set the RHS of the epsilon-constraint
            if f2Sense == MOI.MIN_SENSE
                JuMP.fix(RHS, f2Val - ϵ)
                verbose && println(". Solving with f2 <= ", f2Val - ϵ)
            else
                JuMP.fix(RHS, f2Val + ϵ)
                verbose && println(". Solving with f2 >= ", f2Val + ϵ)
            end

            #And solve again
            JuMP.optimize!(m, ignore_optimize_hook=true)
            status = JuMP.termination_status(m)
        end

        #Sort X_E and Y_N
        s = sortperm(vd.Y_N, by = first)
        vd.X_E, vd.Y_N = vd.X_E[s], vd.Y_N[s]

        #Remove the epsilon constraint
        JuMP.delete(m, eps) ; JuMP.delete(m, RHS)
    else
        return status
    end
    return MOI.OPTIMAL
end

function solve_dicho(m::JuMP.Model, optimizer_factory, round_results, verbose, time_in_secs; args...)
    vd = getvOptData(m)
    empty!(vd.Y_N) ; empty!(vd.X_E)
    f1, f2 = vd.objs
    f1Sense, f2Sense = vd.objSenses
    varArray = JuMP.all_variables(m)

    #Set the first objective as an objective in the JuMP JuMP.Model
    JuMP.set_objective(m, f1Sense, f1)
    verbose && println("solving for z1")

    #Solve with that objective
    JuMP.optimize!(m, optimizer_factory, ignore_optimize_hook=true)
    status = JuMP.termination_status(m)

    #If a solution exists
    if status == MOI.OPTIMAL

        yr_1 = JuMP.value(f1)
        yr_2 = JuMP.value(f2)

        #Store results in vOptData
        push!(vd.Y_N, round_results ? round.([yr_1, yr_2]) : [yr_1, yr_2])
        push!(vd.X_E, JuMP.value.(varArray))

        #Set the second objective as an objective in the JuMP JuMP.Model
        JuMP.set_objective(m, f2Sense, f2)
        verbose && println("solving for z2")

        #Solve with that objective
        JuMP.optimize!(m, ignore_optimize_hook=true)
        status = JuMP.termination_status(m)

        if status == MOI.OPTIMAL

            ys_1 = JuMP.value(f1)
            ys_2 = JuMP.value(f2)

            if !isapprox(yr_1, ys_1, atol=1e-3) || !isapprox(yr_2, ys_2, atol=1e-3)
                push!(vd.Y_N, round_results ? round.([ys_1, ys_2]) : [ys_1, ys_2])
                push!(vd.X_E, JuMP.value.(varArray))
                start = time()
                dichoRecursion(m, yr_1, yr_2, ys_1, ys_2, varArray, round_results, verbose, start, time_in_secs; args...)
            end

            #Sort X_E and Y_N
            s = sortperm(vd.Y_N, by = first)
            vd.Y_N = vd.Y_N[s]
            vd.X_E = vd.X_E[s]

            R1 = f1Sense==MOI.MIN_SENSE ? (<=) : (>=)
            R2 = f2Sense==MOI.MIN_SENSE ? (<=) : (>=)
            weak_dom(a, b) = R1(a[1], b[1]) && R2(a[2], b[2])

            #Filter X_E and Y_N :
            inds = Int[]
            for i = 1:length(vd.Y_N)-1
                if weak_dom(vd.Y_N[i], vd.Y_N[i+1])
                    push!(inds, i+1)
                elseif weak_dom(vd.Y_N[i+1], vd.Y_N[i])
                    push!(inds, i)
                end
            end
            # deleteat!(vd.Y_N, inds)
            # deleteat!(vd.X_E, inds)
        end
    end

    status
end

function dichoRecursion(m::JuMP.Model, yr_1, yr_2, ys_1, ys_2, varArray, round_results, verbose, start, time_in_secs; args...)

    vd = getvOptData(m)
    f1, f2 = vd.objs
    f1Sense, f2Sense = vd.objSenses

    λ1 = abs(yr_2 - ys_2)
    λ2 = abs(ys_1 - yr_1)

    if f1Sense==f2Sense
        lb = λ1*yr_1 + λ2*yr_2
        JuMP.set_objective(m, f1Sense, λ1*f1 + λ2*f2)
        verbose && println("solving for $λ1*f1 + $λ2*f2")
    else
        lb = λ1*yr_1 - λ2*yr_2
        JuMP.set_objective(m, f1Sense, λ1*f1 - λ2*f2)
        verbose && println("solving for $λ1*f1 - $λ2*f2")
    end

    JuMP.optimize!(m, ignore_optimize_hook=true)

    yt_1 = JuMP.value(f1)
    yt_2 = JuMP.value(f2)

    val = f1Sense == f2Sense ? λ1*yt_1 + λ2*yt_2 : λ1*yt_1 - λ2*yt_2

    # check if past time limit, then return
    if (time() - start > time_in_secs)
        push!(vd.Y_N, round_results ? round.([yt_1, yt_2]) : [yt_1, yt_2])
        push!(vd.X_E, JuMP.value.(varArray))
        return
    end

    if (f1Sense == MOI.MIN_SENSE && val < lb - 1e-4) || val > lb + 1e-4
        push!(vd.Y_N, round_results ? round.([yt_1, yt_2]) : [yt_1, yt_2])
        push!(vd.X_E, JuMP.value.(varArray))
        dichoRecursion(m, yr_1, yr_2, yt_1, yt_2, varArray, round_results, verbose, start, time_in_secs; args...)
        dichoRecursion(m, yt_1, yt_2, ys_1, ys_2, varArray, round_results, verbose, start, time_in_secs; args...)
    end

end

function solve_Chalmet(m::JuMP.Model, optimizer_factory, step, verbose ; args...)
    vd = getvOptData(m)
    empty!(vd.Y_N) ; empty!(vd.X_E)
    f1,f2 = vd.objs[1],vd.objs[2]
    f1Sense,f2Sense = vd.objSenses[1],vd.objSenses[2]
    varArray = JuMP.all_variables(m)

    #Set the first objective as an objective in the JuMP JuMP.Model
    JuMP.set_objective(m, f1Sense, f1)
    verbose && println("solving for z1")

    #Solve with that objective
    JuMP.optimize!(m, optimizer_factory, ignore_optimize_hook=true)
    status = JuMP.termination_status(m)

    #If a solution exists
    if status == MOI.OPTIMAL

        yr_1 = JuMP.value(f1)
        yr_2 = JuMP.value(f2)

        #Store results in vOptData
        push!(vd.Y_N, [yr_1, yr_2])
        push!(vd.X_E, JuMP.value.(varArray))

        #Set the second objective as an objective in the JuMP JuMP.Model
        JuMP.set_objective(m, f2Sense, f2)
        verbose && println("solving for z2")

        #Solve with that objective
        JuMP.optimize!(m, ignore_optimize_hook=true)
        status = JuMP.termination_status(m)

        if status == MOI.OPTIMAL

            ys_1 = JuMP.value(f1)
            ys_2 = JuMP.value(f2)

            if !isapprox(yr_1, ys_1, atol=1e-3) || !isapprox(yr_2, ys_2, atol=1e-3)
                push!(vd.Y_N, [ys_1, ys_2])
                push!(vd.X_E, JuMP.value.(varArray))
                #Declare the constraints on z1 and z2 (RHS will be set later)
                rhs_z1 = JuMP.@variable(m)
                rhs_z2 = JuMP.@variable(m)
                cstr_z1 = f1Sense == MOI.MIN_SENSE ? JuMP.@constraint(m, f1 <= rhs_z1) : JuMP.@constraint(m, f1 >= rhs_z1)
                cstr_z2 = f2Sense == MOI.MIN_SENSE ? JuMP.@constraint(m, f2 <= rhs_z2) : JuMP.@constraint(m, f2 >= rhs_z2)

                ChalmetRecursion(m, yr_1, yr_2, ys_1, ys_2, varArray, rhs_z1, rhs_z2, step, verbose; args...)

                #Sort X_E and Y_N
                s = sortperm(vd.Y_N, by = first)
                vd.Y_N = vd.Y_N[s]
                vd.X_E = vd.X_E[s]

                R1 = f1Sense==MOI.MIN_SENSE ? (<=) : (>=)
                R2 = f2Sense==MOI.MIN_SENSE ? (<=) : (>=)
                weak_dom(a, b) = R1(a[1], b[1]) && R2(a[2], b[2])

                #Filter X_E and Y_N :
                inds = Int[]
                for i = 1:length(vd.Y_N)-1
                    if weak_dom(vd.Y_N[i], vd.Y_N[i+1])
                        push!(inds, i+1)
                    elseif weak_dom(vd.Y_N[i+1], vd.Y_N[i])
                        push!(inds, i)
                    end
                end
                deleteat!(vd.Y_N, inds)
                deleteat!(vd.X_E, inds)

                JuMP.delete(m, cstr_z1)
                JuMP.delete(m, cstr_z2)
                JuMP.delete(m, rhs_z1)
                JuMP.delete(m, rhs_z2)
            end
        end
    end

    return MOI.OPTIMAL
end

function ChalmetRecursion(m::JuMP.Model, yr_1, yr_2, ys_1, ys_2, varArray, rhs_z1, rhs_z2, ϵ, verbose ; args...)

    vd = getvOptData(m)
    f1, f2 = vd.objs
    f1Sense, f2Sense = vd.objSenses
    JuMP.set_objective(m, f1Sense, f1Sense==f2Sense ? f1 + f2 : f1 - f2)
    lbz1 = f1Sense==MOI.MAX_SENSE ? min(yr_1, ys_1) : max(yr_1, ys_1)
    lbz2 = f2Sense==MOI.MAX_SENSE ? min(yr_2, ys_2) : max(yr_2, ys_2)

    if f1Sense == MOI.MIN_SENSE
        JuMP.fix(rhs_z1, lbz1 - ϵ)
        verbose && print("solving with z1 <= ", lbz1 - ϵ)
    else
        JuMP.fix(rhs_z1, lbz1 + ϵ)
        verbose && print("solving with z1 >= ", lbz1 + ϵ)
    end

    if f2Sense == MOI.MIN_SENSE
        JuMP.fix(rhs_z2, lbz2 - ϵ)
        verbose && println(" and z2 <= ", lbz2 - ϵ)
    else
        JuMP.fix(rhs_z2, lbz2 + ϵ)
        verbose && println(" and z2 >= ", lbz2 + ϵ)
    end

    JuMP.optimize!(m, ignore_optimize_hook=true)
    status = JuMP.termination_status(m)

    if status == MOI.OPTIMAL
        yt_1 = JuMP.value(f1)
        yt_2 = JuMP.value(f2)
        push!(vd.Y_N, [yt_1, yt_2])
        push!(vd.X_E, JuMP.value.(varArray))
        ChalmetRecursion(m, yr_1, yr_2, yt_1, yt_2, varArray, rhs_z1, rhs_z2, ϵ, verbose; args...)
        ChalmetRecursion(m, yt_1, yt_2, ys_1, ys_2, varArray, rhs_z1, rhs_z2, ϵ, verbose ; args...)
    end
end


text:
using DelimitedFiles

# First Half

exprpt = readdlm("inputs/day01.txt", Int)
[i * j for (i, j) in Iterators.product(exprpt, exprpt) if i + j == 2020][1]

# Second Half

exprpt = readdlm("inputs/day01.txt", Int)
[i * j * k for (i, j, k) in Iterators.product(exprpt, exprpt, exprpt) if i + j + k == 2020][1]


text:
module LibSpirvCross

import SPIRV_Cross_jll: libspirv_cross

using CEnum

const Ctm = Base.Libc.TmStruct
const Ctime_t = UInt
const Cclock_t = UInt

include("spirv_cross_common.jl")
include("spirv_cross_api.jl")

export Ctm, Ctime_t, Cclock_t

# export everything
foreach(names(@__MODULE__, all=true)) do s
   if startswith(string(s), r"spv|SPV|Spv")
      @eval export $s
   end
end

end # module


text:
# This file was generated, do not modify it. # hide
# hideall
c = """
3.26;0.829;1.676;0;1;1.453;3.770
2.189;0.58;0.863;0;0;1.348;3.115
2.125;0.638;0.831;0;0;1.348;3.531
3.027;0.331;1.472;1;0;1.807;3.510
2.094;0.827;0.86;0;0;1.886;5.390
3.222;0.331;2.177;0;0;0.706;1.819
3.179;0;1.063;0;0;2.942;3.947
3;0;0.938;1;0;2.851;3.513
2.62;0.499;0.99;0;0;2.942;4.402
2.834;0.134;0.95;0;0;1.591;3.021
2.405;0.134;0.843;0;0;1.769;3.210
2.728;0.223;0.953;0;0;1.591;2.371
2.512;0.223;0.929;1;0;1.769;3.919
2.834;0.134;1.237;0;0;1.859;3.030
2.819;0.331;1.271;0;1;0.981;2.736
2.126;0.251;1.114;0;0;0.143;2.157
2.834;0.134;1.322;0;0;1.199;2.413
3.014;0.56;1.781;0;0;-0.115;0.898
3.024;0.452;2.698;0;0;1.107;0.450
3.036;0.405;1.205;1;0;1.807;3.733
2.707;0.972;1.889;0;3;-1.169;2.976
2.978;1.246;1.103;0;1;3.988;6.535
3.111;0.732;0.923;0;0;4.068;5.643
"""
fpath, = mktemp()
write(fpath, c);


text:
using StatsPlots
using Test
using StableRNGs
using NaNMath
using Clustering

@testset "Grouped histogram" begin
    rng = StableRNG(1337)
    gpl = groupedhist( rand(rng, 1000), yscale=:log10, ylims=(1e-2, 1e4), bar_position = :stack)
    @test NaNMath.minimum(gpl[1][1][:y]) <= 1e-2
    @test NaNMath.minimum(gpl[1][1][:y]) > 0
    rng = StableRNG(1337)
    gpl = groupedhist( rand(rng, 1000), yscale=:log10, ylims=(1e-2, 1e4), bar_position = :dodge)
    @test NaNMath.minimum(gpl[1][1][:y]) <= 1e-2
    @test NaNMath.minimum(gpl[1][1][:y]) > 0

    data = [1,1,1,1,2,1]
    mask = (collect(1:6) .< 5)
    gpl1 = groupedhist(data[mask], group = mask[mask], color = 1)
    gpl2 = groupedhist(data[.!mask], group = mask[.!mask], color = 2)
    gpl12 = groupedhist(data, group = mask, nbins = 5, bar_position = :stack)
    @test NaNMath.maximum(gpl12[1][2][:y]) == NaNMath.maximum(gpl1[1][1][:y])
    data=[10 12; 1 1; 0.25 0.25]
    gplr = groupedbar(data)
    @test NaNMath.maximum(gplr[1][1][:y]) == 10
    @test NaNMath.maximum(gplr[1][2][:y]) == 12
    gplr = groupedbar(data, bar_position = :stack)
    @test NaNMath.maximum(gplr[1][1][:y]) == 22
    @test NaNMath.maximum(gplr[1][2][:y]) == 12
end # testset

@testset "dendrogram" begin
    # Example from https://en.wikipedia.org/wiki/Complete-linkage_clustering
    wiki_example = [
		0	17	21	31	23
		17	0	30	34	21
		21	30	0	28	39
		31	34	28	0	43
		23	21	39	43	0
	]
    clustering = hclust(wiki_example, linkage=:complete)

    xs, ys = StatsPlots.treepositions(clustering, true, :vertical)

    @test xs == [
        2.0  1.0  4.0  1.75
        2.0  1.0  4.0  1.75
        3.0  2.5  5.0  4.5
        3.0  2.5  5.0  4.5
    ]

    @test ys == [
         0.0   0.0   0.0  23.0
        17.0  23.0  28.0  43.0
        17.0  23.0  28.0  43.0
         0.0  17.0   0.0  28.0
    ]
end


text:
module Dice

include("core.jl")

include("backend/cudd.jl")
include("backend/ir.jl")

include("dsl.jl")
include("ocaml.jl")

# include("lib/int.jl")
# include("lib/tuple.jl")

end # module


text:
using st16SseqJuliaTools
using ArgParse2
using Statistics
using DataStructures
using BioSequences
using BioAlignments

function julia_main()::Cint
    parser = ArgumentParser(prog = "",
                            description = "Filter out contigs based on length median",
                            epilog = "")

    add_argument!(parser, "-i", "--infasta", type = String, help = "Input fasta")
    add_argument!(parser, "-o", "--outfasta", type = String, help = "Output fasta")
    add_argument!(parser, "-e", "--length-cutoff", type = Float64, help = "Fraction from maximumlength - shorter contigs will be discarded", default=0.6)
    add_argument!(parser, "-d", "--identity-cutoff", type = Float64, help = "Aligned fraction identity cut off for containment", default=0.97)
    add_argument!(parser, "-v", "--coverage-cutoff", type = Float64, help = "Aligned part length fraction in the shorter sequence", default=0.6)
    add_argument!(parser, "-l", "--length-3endcut", type = Int64, help = "How many 3end bases should be removed before merging", default=20)
    add_argument!(parser,
    "--filter-length",
    action = "store_true",
    default = false,
    help = "Turn on length filtering based on maximumlength")
    add_argument!(parser,
    "--remove-contained",
    action = "store_true",
    default = false,
    help = "Turn on removal of contained sequences")

    args = parse_args(parser)
    seqs = st16SseqJuliaTools.fasta_to_dict(args.infasta)
    lens = map(x -> length(x),collect(values(seqs)))
    llimit = args.length_cutoff*maximum(lens)
    seqs = sort(seqs, by = x -> length(seqs[x]),rev = true)
    #remove some of 3end
    tmp = OrderedDict{String,BioSequences.LongSequence{BioSequences.DNAAlphabet{4}}}()
    for k in keys(seqs)
        tmp[k] = seqs[k][1:end-args.length_3endcut]
    end
    seqs = tmp

    if args.filter_length
        inict = length(seqs)
        seqs = filter(x -> length(x[2]) >= llimit, seqs)
        gotct = length(seqs)
        dct = inict - gotct
        println(stderr,"Removed $dct sequences using length cutoff")
    end
    #len_filtered = filter(x -> len_sorted[x] >= llimit, len_sorted)
    #st16SseqJuliaTools.dict_to_fasta(len_sorted,args.outfasta)
    ks = collect(keys(seqs))
    len_ks = length(ks)
    out = OrderedDict{String,BioSequences.LongSequence{BioSequences.DNAAlphabet{4}}}()
    out[ks[1]] = seqs[ks[1]] #longest sequence is always kept
    contained = fill(false,len_ks)
    if args.remove_contained
        for i in 1:len_ks
            for j in i:len_ks
                if i != j
                    identityf,coveragef =  check_pair(seqs[ks[i]], seqs[ks[j]])
                    if round(identityf,digits=2) >= args.identity_cutoff && round(coveragef,digits=2) >= args.coverage_cutoff
                        contained[j] = true
                        println("aaaa")
                    end
                end
            end
        end
        contained_nr = count(contained)
        println(stderr,"Removed $contained_nr due to containements")
        for i in 1:len_ks
            k = ks[i]
            cnted = contained[i]
            if !cnted
                out[k]=seqs[k]
            end
        end
        seqs = out
    end
    st16SseqJuliaTools.dict_to_fasta(seqs,args.outfasta)
   return 0

end

julia_main()


text:
function is_equally_spaced(ta::TimeArray)
  timediffs = ta.timestamp[2:end] - ta.timestamp[1:end-1]
  all(timediffs .== timediffs[1])
end # is_equally_spaced

require_equally_spaced(ta::TimeArray) = is_equally_spaced(ta) ? nothing : error("
  Time series observations are not eqally spaced")

is_complete(ta::TimeArray) = !any(isnan(ta.values))

require_complete(ta::TimeArray) = is_complete(ta) ? nothing : error("
  Time series has missing values")

is_clean(ta::TimeArray) = is_equally_spaced(ta) && is_complete(ta)

require_clean(ta::TimeArray) = is_clean(ta) ? nothing : error("
  Time series has missing or unequally spaced values")


text:
using Graphs
using CausalInference
using StatsBase
using Random
using Test
Random.seed!(1)

@testset "backdoor" begin
for j in 1:100
    local d, v, Z
    d = 9
    d2 = 3
    adj = [ i < j ? rand()<0.1 : false for i in 1:d, j in 1:d]

    for i in 1:100
        gd = DiGraph(copy(adj))

        u, v = minmax(sample(1:d, 2, replace=false)...)
        Z = sample(1:d, d2, replace=false)
        (u in Z || v in Z) && continue
        #has_edge(gd, u, v) && continue
        bd = backdoor_criterion(gd, u, v, Z)
        hasap = has_a_path(gd, [u], Z, [])
        @test !(hasap && bd)
        @test hasap == any(has_path(gd, u, z) for z in Z)
        gd0 = DiGraph(copy(adj))
        for w in collect(outneighbors(gd0, u))
            rem_edge!(gd0, Edge(u, w)) || error()
        end
        ds = dsep(gd0, u, v, Z)
        if (ds && !hasap) !=  bd
            backdoor_criterion(gd, u, v, Z, verbose=true)

            @show (ds && !hasap),  bd
            display(adj)
            println(collect(edges(gd)))
            println(collect(edges(gd0)))
            println("$u $v $Z")
            error("")
        end
        @test (ds && !hasap) ==  bd

    end
end

# Pearl Causality Figure 3.4
dg = digraph([
    1 => 3
    3 => 7
    #
    2 => 5
    5 => 8
    #
    7 => 6
    6 => 8
    #
    1 => 4
    2 => 4
    4 => 7
    4 => 8])

    @test !backdoor_criterion(dg, 7, 8)
    @test backdoor_criterion(dg, 7, 8, [3,4])
    @test backdoor_criterion(dg, 7, 8, [4,5])

    # more test examples
    # http://causality.cs.ucla.edu/blog/index.php/category/back-door-criterion/

    dg = digraph([2=>3, 3=>4, 1=>3])
    @test backdoor_criterion(dg, 2, 4)
    @test backdoor_criterion(dg, 2, 4, [1])

    dg = digraph([1=>3, 1=>5, 2=>3, 2=>4, 4=>5])
    @test backdoor_criterion(dg, 2, 3)
    @test !backdoor_criterion(dg, 2, 3, [4])

    dg = digraph([1=>3, 1=>4, 2=>3, 2=>4])
    @test backdoor_criterion(dg, 2, 3)
    @test !backdoor_criterion(dg, 2, 3, [4])

    dg = digraph([1=>2, 2=>3])
    @test backdoor_criterion(dg, 1, 2)
    @test !backdoor_criterion(dg, 1, 2, [3])
    @test !backdoor_criterion(dg, 1, 3, [2])
    @test backdoor_criterion(dg, 1, 3)

    @test backdoor_criterion(dg, 2, 3)
    @test backdoor_criterion(dg, 2, 3, [1])

    dg = digraph([1=>3, 2=>3])
    @test backdoor_criterion(dg, 1, 3)
    @test backdoor_criterion(dg, 1, 3, [2])

    dg = digraph([2=>3, 1=>2, 1=>3])
    @test !backdoor_criterion(dg, 2, 3)
    @test backdoor_criterion(dg, 2, 3, [1])

    dg = digraph([2=>3, 1=>4, 4=>2, 1=>3])
    @test !backdoor_criterion(dg, 2, 3)
    @test backdoor_criterion(dg, 2, 3, [4])

    dg = digraph([2=>3, 1=>2, 1=>4, 4=>3])
    @test !backdoor_criterion(dg, 2, 3)
    @test backdoor_criterion(dg, 2, 3, [4])

end


text:
# voronoi.jl
# algorithm for creating a voronoi tessellation

using GeometryTypes
using StaticArrays

#include("predicates.jl")

imod(a,b) = mod(a-1,b) + 1

function flip23(teta::Tetrahedron, tetb::Tetrahedron)::NTuple{3, Tetrahedron}
    i::Int = 1
    # Faster to use [teta...] than the raw type??
    a::Int = setdiff([teta.vertices...], [tetb.vertices...])[1]
    while teta.vertices[1] != a
        teta = permute(teta,i)
        i += 1
    end
    e::Int = setdiff([tetb.vertices...], [teta.vertices...])[1]
    while tetb.vertices[4] != e
        tetb = permute(tetb,i)
        i += 1
    end
    (Tetrahedron([teta.vertices[2],teta.vertices[4],e,a],
                  teta.neighbors[2],teta.neighbors[4],tetb.neighbors[4],teta.neighbors[1]),
     Tetrahedron([teta.vertices[4],teta.vertices[3],e,a],
                  teta.neighbors[4],teta.neighbors[3],tetb.neighbors[4],teta.neighbors[1]),
     Tetrahedron([teta.vertices[3],teta.vertices[2],e,a],
                  teta.neighbors[3],teta.neighbors[2],tetb.neighbors[4],teta.neighbors[1]))
 end

function flip32(teta::Tetrahedron, tetb::Tetrahedron, tetc::Tetrahedron)::NTuple{2, Tetrahedron}
    # 3 times faster this way!!
    a::Int, e::Int = intersect([teta.vertices...], [tetb.vertices...], [tetc.vertices...])
    i::Int = 1
    while teta.vertices[4] != a
        teta = permute(teta,i)
        i += 1
    end
    while teta.vertices[3] != e
        teta = rotate(teta,4)
    end
    b::Int = teta.vertices[1]
    d::Int = teta.vertices[2]

    while tetb.vertices[4] != a
        tetb = permute(tetb,i)
        i += 1
    end
    while tetb.vertices[3] != e
        tetb = rotate(tetb,4)
    end

    while tetc.vertices[4] != a
        tetc = permute(tetc,i)
        i += 1
    end
    while tetc.vertices[3] != e
        tetc = rotate(tetc,4)
    end
    c::Int = (tetb.vertices[2] == b) ? tetb.vertices[1] : tetc.vertices[1]
    cn::Int = (tetb.vertices[2] == b) ? tetb.neighbors[1] : tetc.neighbors[1]

    (Tetrahedron([a,b,c,d],
                 [teta.neighbors[4],teta.neighbors[1],cn,teta.neighbors[2]]),
     Tetrahedron([b,c,d,e],
                 [teta.neighbors[1],cn,teta.neighbors[2],teta.neighbors[3]]))
 end

function flip14(tet::Tetrahedron, e::Int, offset::Int)::NTuple{4,Tetrahedron}
    (Tetrahedron([tet.vertices[1],tet.vertices[2],tet.vertices[3],e],
                 [offset+3,       offset+2,       offset+4,tet.neighbors[4]]),
     Tetrahedron([tet.vertices[1],tet.vertices[3],tet.vertices[4],e],
                 [offset+3,       offset+4,       offset+1,tet.neighbors[2]]),
     Tetrahedron([tet.vertices[3],tet.vertices[2],tet.vertices[4],e],
                 [offset+4,       offset+2,       offset+1,tet.neighbors[1]]),
     Tetrahedron([tet.vertices[1],tet.vertices[4],tet.vertices[2],e],
                 [offset+3,       offset+1,       offset+2,tet.neighbors[3]]))
end

# TODO: If needed, define a flip41 - not sure that I need the functionality to delete a point
#function flip41(teta::Tetrahedron, tetb::Tetrahedron, tetc::Tetrahedron, tetd::Tetrahedron)::Tetrahedron
#     e::Int = intersect([teta.vertices...],[tetb.vertices...],[tetc.vertices...],[tetd.vertices...])[1]
#     i:Int = 1
#     while teta.vertices[4] != e
#         teta = permute(teta,i)
#         i += 1
#     end
#     while tetb.vertices[4] != e
#         tetb = permute(tetb,i)
#         i += 1
#     end
#     while tetc.vertices[4] != e
#         tetc = permute(tetc,i)
#         i += 1
#     end
#     while tetd.vertices[4] != e
#         tetd = permute(tetd,i)
#         i += 1
#     end
#     d::Int = setdiff([tetb.vertices...],[teta.vertices...])[1]
#     Tetrahedron([teta.vertices[1],teta.vertices[2],teta.vertices[3],d],
#                 [])
#end

function flip44(teta::Tetrahedron, tetb::Tetrahedron, tetc::Tetrahedron, tetd::Tetrahedron)::NTuple{4,Tetrahedron}
    # Tets need to be in the correct order prior to this function!!
    # abcd, aced, fecd, fcbd, where bced is the degenerate plane
    a::Int, b::Int, c::Int, d::Int = teta.vertices
    f::Int = tetc.vertices[1]
    e::Int = tetc.vertices[2]
    (Tetrahedron([a,b,e,d],
                 [teta.neighbors[1],teta.neighbors[2],tetc.neighbors[2],teta.neighbors[4]]),
     Tetrahedron([a,b,c,e],
                 [teta.neighbors[1],teta.neighbors[2],teta.neighbors[3],tetc.neighbors[2]]),
     Tetrahedron([f,d,e,b],
                 [tetc.neighbors[1],teta.neighbors[4],tetc.neighbors[2],teta.neighbors[2]]),
     Tetrahedron([f,e,c,b],
                 [tetc.neighbors[1],tetc.neighbors[2],teta.neighbors[3],teta.neighbors[2]]))
end


text:
export assign_chunk!,
       distribute_task!

"""
    distribute_task!(S::LazySet{N}, v::SharedVector{N}...) where {N<:Real}

Distribute the assignment of each chunk among the available processes.

### Input

- `S` -- convex set
- `v` -- variable number of shared vectors

### Output

Nothing.

### Notes

Use this function to distribute a given task acting on a set `S` and a pool `v`
of shared vectors.

The task for each processor is distributed through `remotecall_wait` using a
function `assign_chunk!` that should be defined elsewhere. The vectors `v`
contain one or more shared vectors in which the values of the task are written.

Typically, the function `assign_chunk!` is a wrapper around some problem-specific
`process_chunk!` function.
"""
function distribute_task!(S::LazySet{N}, v::SharedVector{N}...) where {N<:Real}
    @sync begin
        for p in procs(v[1])
            @async remotecall_wait(assign_chunk!, p, S, v...)
        end
    end
end

"""
    assign_chunk!(S::LazySet{N}, v::SharedVector{N}...) where {N<:Real}

Return the function that assigns the work for each process.

### Input

- `S` -- convex set
- `v` -- variable number of shared vectors

### Output

The function `process_chunk!` that equally distributes the load for each worker.

### Notes

Use this function to distribute a given task acting on a set `S` and a pool `v`
of shared vectors. The tasks are equally distributed among the number of processes.

See also [`distribute_task!`](@ref).
"""
function assign_chunk!(S::LazySet{N}, v::SharedVector{N}...) where {N<:Real}
    return process_chunk!(S, _prange(v[1]), v...)
end

"""
    _prange(v::SharedVector{N}) where {N<:Real}

Returns the indexes assigned to a process.

### Input

- `v` -- shared vector of length `n`

### Output

The indices range assigned to each process.

### Notes

The indices are assigned such that the vector is equally distributed among the
processes. If the worker is not assigned a piece, the unit range `1:0` is
returned.
"""
function _prange(v::SharedVector{N}) where {N<:Real}
    idx = indexpids(v)
    if idx == 0
        # this worker is not assigned a piece
        return 1:0
    end
    nchunks = length(procs(v))
    splits = [round(Int, s) for s in range(0, stop=length(v), length=nchunks+1)]
    splits[idx]+1:splits[idx+1]
end


text:
function dicttest()
  com = LanguageDynamics.ZipfTravellerPopulation{LanguageDynamics.MomentumSelector}(7, 2.5)

  for i = 1:200
    LanguageDynamics.inject!(LanguageDynamics.MomentumSelector(0.01, 0.02, 1.0, 4, 0.1), com)
  end

  #for i = 1:100
  #  eject!(com)
  #end

  while com.nactions < 0.5*200*100#_000
    LanguageDynamics.travel!(com)
    LanguageDynamics.rendezvous!(com)
  end
end

function arrtest()
  com = PoolPopulation{EmptyInteractor}()

  for i = 1:500
    inject!(EmptyInteractor(), com)
  end

  #for i = 1:100
  #  eject!(com)
  #end

  for i = 1:1000
    rendezvous!(com)
  end
end


text:
# using DataStructures

# # Arthur Choi, Guy Van den Broeck, and Adnan Darwiche. Tractable learning for structured probability
# # spaces: A case study in learning preference distributions. In Proceedings of IJCAI, 2015.
# "Calculate the probability of the logic formula given by sdd for the psdd"
# function pr_constraint(psdd_node::ProbΔNode, sdd_node::Union{ProbΔNode, StructLogicalΔNode})
#     cache = Dict{Tuple{ProbΔNode, Union{ProbΔNode, StructLogicalΔNode}}, Float64}()

#     return pr_constraint(psdd_node, sdd_node, cache)
# end
# function pr_constraint(psdd_node::ProbΔNode, sdd_node::Union{ProbΔNode, StructLogicalΔNode},
#                        cache::Dict{Tuple{ProbΔNode, Union{ProbΔNode, StructLogicalΔNode}}, Float64})::Float64
#     if (psdd_node, sdd_node) in keys(cache) # Cache hit
#         return cache[psdd_node, sdd_node]
#     elseif psdd_node isa ProbLiteral # Boundary cases
#         if sdd_node isa Union{ProbLiteral, StructLiteralNode} # Both are literals, just check whether they agrees with each other
#             if literal(psdd_node) == literal(sdd_node)
#                 return get!(cache, (psdd_node, sdd_node), 1.0)
#             else
#                 return get!(cache, (psdd_node, sdd_node), 0.0)
#             end
#         else
#             pr_constraint(psdd_node, sdd_node.children[1], cache)
#             if length(sdd_node.children) > 1
#                 pr_constraint(psdd_node, sdd_node.children[2], cache)
#                 return get!(cache, (psdd_node, sdd_node), 1.0)
#             else
#                 return get!(cache, (psdd_node, sdd_node),
#                     literal(sdd_node.children[1]) == literal(psdd_node) ? 1.0 : 0.0
#                 )
#             end
#         end
#     elseif psdd_node.children[1] isa ProbLiteral # The psdd is true
#         theta = exp(psdd_node.log_thetas[1])
#         return get!(cache, (psdd_node, sdd_node),
#             theta * pr_constraint(psdd_node.children[1], sdd_node, cache) +
#             (1.0 - theta) * pr_constraint(psdd_node.children[2], sdd_node, cache)
#         )
#     else # Both psdds are not trivial
#         prob = 0.0
#         for (prob⋀_node, log_theta) in zip(psdd_node.children, psdd_node.log_thetas)
#             p = prob⋀_node.children[1]
#             s = prob⋀_node.children[2]

#             theta = exp(log_theta)
#             for sdd⋀_node in sdd_node.children
#                 r = sdd⋀_node.children[1]
#                 t = sdd⋀_node.children[2]
#                 prob += theta * pr_constraint(p, r, cache) * pr_constraint(s, t, cache)
#             end
#         end
#         return get!(cache, (psdd_node, sdd_node), prob)
#     end
# end

# "Entropy of the distribution of the input psdd."
# function psdd_entropy(psdd_node::ProbΔNode)::Float64
#     psdd_entropy_cache = Dict{ProbΔNode, Float64}()

#     return psdd_entropy(psdd_node, psdd_entropy_cache)
# end
# function psdd_entropy(psdd_node::Prob⋁, psdd_entropy_cache::Dict{ProbΔNode, Float64})::Float64
#     if psdd_node in keys(psdd_entropy_cache)
#         return psdd_entropy_cache[psdd_node]
#     elseif psdd_node.children[1] isa ProbLiteral
#         return get!(psdd_entropy_cache, psdd_node,
#             - exp(psdd_node.log_thetas[1]) * psdd_node.log_thetas[1] -
#             exp(psdd_node.log_thetas[2]) * psdd_node.log_thetas[2]
#         )
#     else
#         local_entropy = 0.0
#         for (prob⋀_node, log_prob) in zip(psdd_node.children, psdd_node.log_thetas)
#             p = prob⋀_node.children[1]
#             s = prob⋀_node.children[2]

#             local_entropy += exp(log_prob) * (psdd_entropy(p, psdd_entropy_cache) +
#                 psdd_entropy(s, psdd_entropy_cache) - log_prob)
#         end
#         return get!(psdd_entropy_cache, psdd_node, local_entropy)
#     end
# end
# function psdd_entropy(psdd_node::Prob⋀, psdd_entropy_cache::Dict{ProbΔNode, Float64})::Float64
#     return get!(psdd_entropy_cache, psdd_node.children[1], psdd_entropy(psdd_node.children[1], psdd_entropy_cache)) +
#         get!(psdd_entropy_cache, psdd_node.children[2], psdd_entropy(psdd_node.children[2], psdd_entropy_cache))
# end
# function psdd_entropy(psdd_node::ProbLiteral, psdd_entropy_cache::Dict{ProbΔNode, Float64})::Float64
#     return get!(psdd_entropy_cache, psdd_node, 0.0)
# end

# "KL divergence calculation for psdds that are not necessarily identical"
# function psdd_kl_divergence(psdd_node1::ProbΔNode, psdd_node2::ProbΔNode)::Float64
#     kl_divergence_cache = Dict{Tuple{ProbΔNode, ProbΔNode}, Float64}()
#     pr_constraint_cache = Dict{Tuple{ProbΔNode, Union{ProbΔNode, StructLogicalΔNode}}, Float64}()

#     return psdd_kl_divergence(psdd_node1, psdd_node2, kl_divergence_cache, pr_constraint_cache)
# end
# function psdd_kl_divergence(psdd_node1::ProbΔNode, psdd_node2::ProbΔNode,
#         kl_divergence_cache::Dict{Tuple{ProbΔNode, ProbΔNode}, Float64})::Float64
#     pr_constraint_cache = Dict{Tuple{ProbΔNode, Union{ProbΔNode, StructLogicalΔNode}}, Float64}()

#     return psdd_kl_divergence(psdd_node1, psdd_node2, kl_divergence_cache, pr_constraint_cache)
# end
# function psdd_kl_divergence(psdd_node1::ProbΔNode, psdd_node2::ProbΔNode,
#         kl_divergence_cache::Dict{Tuple{ProbΔNode, ProbΔNode}, Float64},
#         pr_constraint_cache::Dict{Tuple{ProbΔNode, Union{ProbΔNode, StructLogicalΔNode}}, Float64})
#     @assert !(psdd_node1 isa Prob⋀ || psdd_node2 isa Prob⋀) "Prob⋀ not a valid PSDD node for KL-Divergence"

#     # Check if both nodes are normalized for same vtree node
#     @assert variables(psdd_node1.origin.vtree) == variables(psdd_node2.origin.vtree) "Both nodes not normalized for same vtree node"

#     if (psdd_node1, psdd_node2) in keys(kl_divergence_cache) # Cache hit
#         return kl_divergence_cache[(psdd_node1, psdd_node2)]
#     elseif psdd_node1.children[1] isa ProbLiteral
#         if psdd_node2 isa ProbLiteral
#             psdd_kl_divergence(psdd_node1.children[1], psdd_node2, kl_divergence_cache, pr_constraint_cache)
#             psdd_kl_divergence(psdd_node1.children[2], psdd_node2, kl_divergence_cache, pr_constraint_cache)
#             if literal(psdd_node1.children[1]) == literal(psdd_node2)
#                 return get!(kl_divergence_cache, (psdd_node1, psdd_node2),
#                     psdd_node1.log_thetas[1] * exp(psdd_node1.log_thetas[1])
#                 )
#             else
#                 return get!(kl_divergence_cache, (psdd_node1, psdd_node2),
#                     psdd_node1.log_thetas[2] * exp(psdd_node1.log_thetas[2])
#                 )
#             end
#         else
#             # The below four lines actually assign zero, but still we need to
#             # call it.
#             psdd_kl_divergence(psdd_node1.children[1], psdd_node2.children[1], kl_divergence_cache, pr_constraint_cache)
#             psdd_kl_divergence(psdd_node1.children[1], psdd_node2.children[2], kl_divergence_cache, pr_constraint_cache)
#             psdd_kl_divergence(psdd_node1.children[2], psdd_node2.children[1], kl_divergence_cache, pr_constraint_cache)
#             psdd_kl_divergence(psdd_node1.children[2], psdd_node2.children[2], kl_divergence_cache, pr_constraint_cache)
#             # There are two possible matches
#             if literal(psdd_node1.children[1]) == literal(psdd_node2.children[1])
#                 return get!(kl_divergence_cache, (psdd_node1, psdd_node2),
#                     exp(psdd_node1.log_thetas[1]) * (psdd_node1.log_thetas[1] - psdd_node2.log_thetas[1]) +
#                     exp(psdd_node1.log_thetas[2]) * (psdd_node1.log_thetas[2] - psdd_node2.log_thetas[2])
#                 )
#             else
#                 return get!(kl_divergence_cache, (psdd_node1, psdd_node2),
#                     exp(psdd_node1.log_thetas[1]) * (psdd_node1.log_thetas[1] - psdd_node2.log_thetas[2]) +
#                     exp(psdd_node1.log_thetas[2]) * (psdd_node1.log_thetas[2] - psdd_node2.log_thetas[1])
#                 )
#             end
#         end
#     else # the normal case
#         kl_divergence = 0.0

#         # loop through every combination of prim and sub
#         for (prob⋀_node1, log_theta1) in zip(psdd_node1.children, psdd_node1.log_thetas)
#             for (prob⋀_node2, log_theta2) in zip(psdd_node2.children, psdd_node2.log_thetas)
#                 p = prob⋀_node1.children[1]
#                 s = prob⋀_node1.children[2]

#                 r = prob⋀_node2.children[1]
#                 t = prob⋀_node2.children[2]

#                 theta1 = exp(log_theta1)

#                 p11 = pr_constraint(s, t, pr_constraint_cache)
#                 p12 = pr_constraint(p, r, pr_constraint_cache)

#                 p13 = theta1 * (log_theta1 - log_theta2)

#                 p21 = psdd_kl_divergence(p, r, kl_divergence_cache, pr_constraint_cache)
#                 p31 = psdd_kl_divergence(s, t, kl_divergence_cache, pr_constraint_cache)

#                 kl_divergence += p11 * p12 * p13 + theta1 * (p11 * p21 + p12 * p31)
#             end
#         end
#         return get!(kl_divergence_cache, (psdd_node1, psdd_node2), kl_divergence)
#     end
# end
# function psdd_kl_divergence(psdd_node1::ProbLiteral, psdd_node2::ProbLiteral,
#         kl_divergence_cache::Dict{Tuple{ProbΔNode, ProbΔNode}, Float64},
#         pr_constraint_cache::Dict{Tuple{ProbΔNode, Union{ProbΔNode, StructLogicalΔNode}}, Float64})
#     # Check if literals are over same variables in vtree
#    @assert variables(psdd_node1.origin.vtree) == variables(psdd_node2.origin.vtree) "Both nodes not normalized for same vtree node"

#     if (psdd_node1, psdd_node2) in keys(kl_divergence_cache) # Cache hit
#         return kl_divergence_cache[psdd_node1, psdd_node2]
#     else
#         # In this case probability is 1, kl divergence is 0
#         return get!(kl_divergence_cache, (psdd_node1, psdd_node2), 0.0)
#     end
# end
# function psdd_kl_divergence(psdd_node1::Prob⋁, psdd_node2::ProbLiteral,
#         kl_divergence_cache::Dict{Tuple{ProbΔNode, ProbΔNode}, Float64},
#         pr_constraint_cache::Dict{Tuple{ProbΔNode, Union{ProbΔNode, StructLogicalΔNode}}, Float64})
#     @assert variables(psdd_node1.origin.vtree) == variables(psdd_node2.origin.vtree) "Both nodes not normalized for same vtree node"

#     if (psdd_node1, psdd_node2) in keys(kl_divergence_cache) # Cache hit
#         return kl_divergence_cache[psdd_node1, psdd_node2]
#     else
#         psdd_kl_divergence(psdd_node1.children[1], psdd_node2, kl_divergence_cache, pr_constraint_cache)
#         psdd_kl_divergence(psdd_node1.children[2], psdd_node2, kl_divergence_cache, pr_constraint_cache)
#         if literal(psdd_node1.children[1]) == literal(psdd_node2)
#             return get!(kl_divergence_cache, (psdd_node1, psdd_node2),
#                 psdd_node1.log_thetas[1] * exp(psdd_node1.log_thetas[1])
#             )
#         else
#             return get!(kl_divergence_cache, (psdd_node1, psdd_node2),
#                 psdd_node1.log_thetas[2] * exp(psdd_node1.log_thetas[2])
#             )
#         end
#     end
# end
# function psdd_kl_divergence(psdd_node1::ProbLiteral, psdd_node2::Prob⋁,
#         kl_divergence_cache::Dict{Tuple{ProbΔNode, ProbΔNode}, Float64},
#         pr_constraint_cache::Dict{Tuple{ProbΔNode, Union{ProbΔNode, StructLogicalΔNode}}, Float64})
#     @assert variables(psdd_node1.origin.vtree) == variables(psdd_node2.origin.vtree) "Both nodes not normalized for same vtree node"

#     if (psdd_node1, psdd_node2) in keys(kl_divergence_cache) # Cache hit
#         return kl_divergence_cache[psdd_node1, psdd_node2]
#     else
#         psdd_kl_divergence(psdd_node1, psdd_node2.children[1], kl_divergence_cache, pr_constraint_cache)
#         psdd_kl_divergence(psdd_node1, psdd_node2.children[2], kl_divergence_cache, pr_constraint_cache)
#         if literal(psdd_node1) == literal(psdd_node2.children[1])
#             return get!(kl_divergence_cache, (psdd_node1, psdd_node2),
#                 -psdd_node2.log_thetas[1]
#             )
#         else
#             return get!(kl_divergence_cache, (psdd_node1, psdd_node2),
#                 -psdd_node2.log_thetas[2]
#             )
#         end
#     end
# end


text:
"""
Helper infrastructure to compile and sample models using `cmdstan`.

[`StanModel`](@ref) wraps a model definition (source code), while [`stan_sample`](@ref) can
be used to sample from it.

[`stan_compile`](@ref) can be used to pre-compile a model without sampling, or for
debugging. A [`StanModelError`](@ref) is thrown if this fails, which contains the error
messages from `stanc`.

# Example

```julia
using StanRun
model = StanModel("/path/to/model.stan") # directory should be writable, for compilation
data = (N = 100, x = randn(N, 1000))     # in a format supported by stan_dump
chains = stan_sample(model, data, 5)     # 5 chain paths and log files
```
"""
module StanRun

export StanModel, StanModelError, stan_sample, stan_compile

using ArgCheck: @argcheck
using Distributed: pmap
using DocStringExtensions: FIELDS, SIGNATURES, TYPEDEF
using UnPack: @unpack
using StanDump: stan_dump

const CMDSTAN_HOME_VAR = "JULIA_CMDSTAN_HOME"

function get_cmdstan_home()
    get(ENV, CMDSTAN_HOME_VAR) do
        throw(ErrorException("The environment variable $CMDSTAN_HOME_VAR needs to be set."))
    end
end

struct StanModel{S <: AbstractString}
    source_path::S
    cmdstan_home::S
end

function Base.show(io::IO, model::StanModel)
    @unpack source_path, cmdstan_home = model
    print(io, "Stan model at $(source_path)",
          "\n    (CmdStan home: $(model.cmdstan_home))")
end

"""
$(SIGNATURES)

Replace the extension of `path` (including the `'.'`) with `new_ext`, which can be any
string (not necessarily an extension, with the dot).

When `verified_ext` is given, it the original extension is checked to be equivalent.
Defaults to `".stan"`.

Internal, not exported.
"""
function replace_ext(path::AbstractString, new_ext, verified_ext = ".stan")
    basename, ext = splitext(path)
    verified_ext ≢ nothing && @argcheck ext == verified_ext
    basename * new_ext
end

"""
$(SIGNATURES)

Executable path corresponding to a source file, or a model.

Internal, not exported.
"""
function executable_path(source_path::AbstractString)
    replace_ext(source_path, Sys.iswindows() ? ".exe" : "")
end

executable_path(model::StanModel) = executable_path(model.source_path)

"""
$(SIGNATURES)

Define a model by its Stan source code location, which needs to end in `".stan"`.

Its directory needs to be *writable*, as it will contain the compiled executable (generated
on demand if it does not exist, or if the source code is more recent).

`cmdstan_home` should specify the directory of the `cmdstan` installation. The default is
obtained from the environment variable `$(CMDSTAN_HOME_VAR)`."
"""
function StanModel(source_path; cmdstan_home = get_cmdstan_home())
    StanModel(source_path, cmdstan_home)
end

"""
$(TYPEDEF)

Error thrown when a Stan model fails to compile. Accessing fields directly is part of the
API.

$(FIELDS)
"""
struct StanModelError <: Exception
    model::StanModel
    message::String
end

function Base.showerror(io::IO, e::StanModelError)
    print(io, "error when compiling ", e.model, ":\n",
          e.message)
end

"""
$(SIGNATURES)

Ensure that a compiled model executable exists, and return its path.

If compilation fails, a `StanModelError` is returned instead.

See [`stan_compile`](@ref) for the documentation of keyword arguments.

Internal, not exported.
"""
function ensure_executable(model::StanModel; debug::Bool = false, dry_run::Bool = false)
    @unpack cmdstan_home = model
    exec_path = executable_path(model)
    error_output = IOBuffer()
    cmd = `make -f $(cmdstan_home)/makefile -C $(cmdstan_home) $(exec_path)`
    if debug
        @info "Stan compilation information" cmdstan_home cmd exec_path
    end
    if dry_run
        nothing
    else
        is_ok = cd(cmdstan_home) do
            success(pipeline(cmd; stderr = error_output))
        end
        if is_ok
            exec_path
        else
            throw(StanModelError(model, String(take!(error_output))))
        end
    end
end

"""
$(SIGNATURES)

Default `output_base`, in the same directory as the model. Internal, not exported.
"""
default_output_base(model::StanModel) = replace_ext(model.source_path, "", ".stan")

sample_file_path(output_base::AbstractString, id::Int) = output_base * "_chain_$(id).csv"

log_file_path(output_base::AbstractString, id::Int) = output_base * "_chain_$(id).log"

"""
$(SIGNATURES)

Make a Stan command. Internal, not exported.
"""
function stan_cmd_and_paths(exec_path::AbstractString, data_file::AbstractString,
                            output_base::AbstractString, id::Integer,
                            sample_options, output_options)
    sample_file = sample_file_path(output_base, id)
    log_file = log_file_path(output_base, id)
    pipeline(`$(exec_path) sample id=$(id) $(sample_options) data file=$(data_file) output file=$(sample_file) $(output_options)`;
             stdout = log_file), (sample_file, log_file)
end

"""
$(SIGNATURES)

Compile a model, throwing an error if it failed.

When `debug = true`, write the command that would be executed into the log.

When `dry_run = true`, don't compile.
"""
function stan_compile(model; debug::Bool = false, dry_run::Bool = false)
    ensure_executable(model; debug = debug, dry_run = dry_run)
    nothing
end

function stan_sample(model, data::NamedTuple, n_chains::Integer;
                     output_base = default_output_base(model),
                     data_file = output_base * ".data.R",
                     rm_samples = true, sample_options = (), output_options = ())
    stan_dump(data_file, data; force = true)
    stan_sample(model, data_file, n_chains; output_base = output_base,
                rm_samples = rm_samples, sample_options = sample_options,
                output_options = output_options)
end

"""
$(SIGNATURES)

Sample `n_chains` from `model` using `data_file`. Return the full paths of the sample files
and logs as pairs. In case of an error with a chain, the first value is `nothing`.

`output_base` is used to write the data file (using `StanDump.stan_dump`) and to determine
the resulting names for the sampler output. It defaults to the source file name without the
extension.

When `data` is provided as a `NamedTuple`, it is written using `StanDump.stan_dump` first.

When `rm_samples` (default: `true`), remove potential pre-existing sample files after
compiling the model.

`sample_options` and `output_options` are either strings, or iterables (empty by default),
and are pasted in after `sample` or `output`, respectively, in the command line.
"""
function stan_sample(model::StanModel, data_file::AbstractString, n_chains::Integer;
                     output_base = default_output_base(model),
                     rm_samples = true, sample_options = (), output_options = ())
    exec_path = ensure_executable(model)
    rm_samples && rm.(find_samples(model))
    cmds_and_paths = [stan_cmd_and_paths(exec_path, data_file, output_base, id,
                                         sample_options, output_options)
                      for id in 1:n_chains]
    pmap(cmds_and_paths) do cmd_and_path
        cmd, (sample_path, log_path) = cmd_and_path
        success(cmd) ? sample_path : nothing, log_path
    end
end

"""
$(SIGNATURES)

Return filenames of CSV files (with MCMC samples, this is not checked) matching
`output_base` from the model.

Part of the API, but not exported.
"""
function find_samples(output_base::AbstractString)
    dir, basename = splitdir(output_base)
    rx = Regex(basename * raw"_chain_\d+.csv")
    joinpath.(Ref(dir), filter(file -> occursin(rx, file), readdir(dir)))
end

find_samples(model::StanModel) = find_samples(default_output_base(model))

end # module


text:
function TransSphere(ksi,eta,z,F,Topo,Param)
X=zeros(1,3);
X[1:3]=0.25*((1-ksi)*(1-eta)*F.P[1:3,1]+
  (1+ksi)*(1-eta)*F.P[1:3,2]+
  (1+ksi)*(1+eta)*F.P[1:3,3]+
  (1-ksi)*(1+eta)*F.P[1:3,4]);
r=norm(X);
X=X/r;
X=X*(Param.RadPrint+z);
return X
end


text:
# Load Julia packages (libraries) needed for clip

using StatisticalRethinking, StructuralCausalModels
using StatsPlots

ProjDir = @__DIR__

println()
N = 100
df = DataFrame();
df[!, :A] = rand(Normal(), N)
df[!, :M] = [rand(Normal(-a), 1)[1] for a in df[:, :A]]
df[!, :D] = [rand(Normal(m), 1)[1] for m in df[:, :M]]

StatsPlots.cornerplot(Array(df), label=names(df))
savefig("$(ProjDir)/AMD_df_1.png")

println()

for i in 1:3
  include(scm_path("..", "examples", "SR", "chain", "AMD_m$i.jl"))
end

if success(rc)

  r1 = plotcoef([m5_1s, m5_2s, m5_3s], [:bA, :bM], "$(ProjDir)/AMD_chain_1.png",
    "Particles (Normal) estimates")
  display(r1)

end

d = OrderedDict(
  :M => [:A],
  :D => [:M]
);

dag = DAG("chain", d; df=df);
show(dag)

display(dag.s); println()

shipley_t = shipley_test(dag)
if !isnothing(shipley_t)
  display(shipley_t)
end
println()

f = [:A]; s = [:D]; sel = vcat(f, s)
c = [:M]

e = d_separation(dag, f, s)
println("d_separation($(dag.name), $f, $s) = $e")

e = d_separation(dag, f, s; c=c)
println("d_separation($(dag.name), $f, $s; c=$c) = $e")


text:
### A Pluto.jl notebook ###
# v0.14.1

using Markdown
using InteractiveUtils

# ╔═╡ e1a58fce-68e7-11eb-28e0-1dfcb1dd2da9
begin
	push!(LOAD_PATH, "../../../../Julia/Julia_Projects/02_modules/")
	using Letters
end

# ╔═╡ eb6d37b4-68e7-11eb-2edb-75552408475d
begin
	push!(LOAD_PATH, "../../../../Julia/Julia_Projects/02_modules/")
	import Numbers
end

# ╔═╡ f64b4784-68e7-11eb-063c-b373fbf05182
include("../../../../Julia/Julia_Projects/02_modules/Letters.jl")

# ╔═╡ c559673c-68e7-11eb-1d60-f53b2158b21c
md"""
### Loading modules with using
"""

# ╔═╡ 0c4ce83a-68e8-11eb-2677-eb6c5dea141e
names(Letters)

# ╔═╡ 1104b66e-68e8-11eb-3283-0387f513e1b5
names(Letters, all=true)

# ╔═╡ 15a40d8c-68e8-11eb-2cc3-fb7be8f74647
Letters.rand_string()

# ╔═╡ 1ca2cd9e-68e8-11eb-100c-fdee2f7cd4b2
begin
	# using Letters: my_name, MY_NAME # from here we no longer need to prefix these "objects" with Letters
	# my_name
	# my_name()
end

# ╔═╡ e251a020-68e7-11eb-2f5f-71145240b647
md"""
### Loading modules with import
"""

# ╔═╡ 36a82734-68e8-11eb-0570-69c26a8f6aea
names(Numbers)

# ╔═╡ 3bc25f70-68e8-11eb-2476-8f0b891d8f09
names(Numbers,all =true)

# ╔═╡ 3ba009fa-68e8-11eb-212c-c139505444bb
Numbers.half_rand()

# ╔═╡ 3b61ec42-68e8-11eb-240d-c7211691226a
Numbers.MY_NAME

# ╔═╡ 3b09e948-68e8-11eb-274e-afc96a0f2c36
Numbers.my_name()

# ╔═╡ ec2b956a-68e7-11eb-0cae-ad84f259746a
md"""
### Loading modules with include
"""

# ╔═╡ 6bec9754-68e8-11eb-3424-454616577b81
# using Letters

# ╔═╡ 7e24e980-68e8-11eb-1f06-d1e1e3ac8f12
# using Main.Letters

# ╔═╡ 8704a19e-68e8-11eb-29ea-57bac063da49
# rand_string()

# ╔═╡ Cell order:
# ╟─c559673c-68e7-11eb-1d60-f53b2158b21c
# ╠═e1a58fce-68e7-11eb-28e0-1dfcb1dd2da9
# ╠═0c4ce83a-68e8-11eb-2677-eb6c5dea141e
# ╠═1104b66e-68e8-11eb-3283-0387f513e1b5
# ╠═15a40d8c-68e8-11eb-2cc3-fb7be8f74647
# ╠═1ca2cd9e-68e8-11eb-100c-fdee2f7cd4b2
# ╟─e251a020-68e7-11eb-2f5f-71145240b647
# ╠═eb6d37b4-68e7-11eb-2edb-75552408475d
# ╠═36a82734-68e8-11eb-0570-69c26a8f6aea
# ╠═3bc25f70-68e8-11eb-2476-8f0b891d8f09
# ╠═3ba009fa-68e8-11eb-212c-c139505444bb
# ╠═3b61ec42-68e8-11eb-240d-c7211691226a
# ╠═3b09e948-68e8-11eb-274e-afc96a0f2c36
# ╟─ec2b956a-68e7-11eb-0cae-ad84f259746a
# ╠═f64b4784-68e7-11eb-063c-b373fbf05182
# ╠═6bec9754-68e8-11eb-3424-454616577b81
# ╠═7e24e980-68e8-11eb-1f06-d1e1e3ac8f12
# ╠═8704a19e-68e8-11eb-29ea-57bac063da49

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks


text:
abstract type AbstractBC{T} <: AbstractDiffEqLinearOperator{T} end

abstract type AtomicBC{T} <: AbstractBC{T} end

"""
Robin, General, and in general Neumann, Dirichlet and Bridge BCs
are not necessarily linear operators.  Instead, they are affine
operators, with a constant term Q*x = Qa*x + Qb.
"""
abstract type AffineBC{T} <: AtomicBC{T} end

struct NeumannBC{N} end
struct Neumann0BC{N} end
struct DirichletBC{N} end
struct Dirichlet0BC{N} end
struct PeriodicBC{T} <: AtomicBC{T}
    PeriodicBC(T::Type) = new{T}()
end
struct RobinBC{T, V<:AbstractVector{T}} <: AffineBC{T}
    a_l::V
    b_l::T
    a_r::V
    b_r::T
    function RobinBC(l::NTuple{3,T}, r::NTuple{3,T}, dx::T, order = 1) where {T}
        αl, βl, γl = l
        αr, βr, γr = r

        s = calculate_weights(1, one(T), Array(one(T):convert(T,order+1))) #generate derivative coefficients about the boundary of required approximation order

        a_l = -s[2:end]./(αl*dx/βl + s[1])
        a_r = s[end:-1:2]./(αr*dx/βr - s[1]) # for other boundary stencil is flippedlr with *opposite sign*

        b_l = γl/(αl+βl*s[1]/dx)
        b_r = γr/(αr-βr*s[1]/dx)

        return new{T, typeof(a_l)}(a_l, b_l, a_r, b_r)
    end
    function RobinBC(l::AbstractVector{T}, r::AbstractVector{T}, dx::AbstractVector{T}, order = 1) where {T}
        αl, βl, γl = l
        αr, βr, γr = r

        s_index = Array(one(T):convert(T,order+1))
        dx_l, dx_r = dx[1:length(s_index)], dx[(end-length(s_index)+1):end]

        s = calculate_weights(1, one(T), s_index) #generate derivative coefficients about the boundary of required approximation order
        denom_l = αl+βl*s[1]/dx_l[1]
        denom_r = αr-βr*s[1]/dx_r[end]

        a_l = -βl.*s[2:end]./(denom_l*dx_l[2:end])
        a_r = βr.*s[end:-1:2]./(denom_r*dx_r[1:(end-1)]) # for other boundary stencil is flippedlr with *opposite sign*

        b_l = γl/denom_l
        b_r = γr/denom_r

        return new{T, typeof(a_l)}(a_l, b_l, a_r, b_r)
    end
end

"""
q = GeneralBC(α_leftboundary, α_rightboundary, dx::T, approximation_order)

-------------------------------------------------------------------------------------

Implements a generalization of the Robin boundary condition, where α is a vector of coefficients.
Represents a condition of the form α[1] + α[2]u[0] + α[3]u'[0] + α[4]u''[0]+... = 0
Implemented in a similar way to the RobinBC (see above).
This time there are multiple stencils for multiple derivative orders - these can be written as a matrix S.
All components that multiply u(0) are factored out, turns out to only involve the first column of S, s̄0. The rest of S is denoted S`. the coeff of u(0) is s̄0⋅ᾱ[3:end] + α[2].
the remaining components turn out to be ᾱ[3:end]⋅(S`ū`) or equivalently (transpose(ᾱ[3:end])*S`)⋅ū`. Rearranging, a stencil q_a to be dotted with ū` upon extension can readily be found, along with a constant component q_b
"""
struct GeneralBC{T, L<:AbstractVector{T}, R<:AbstractVector{T}} <:AffineBC{T}
    a_l::L
    b_l::T
    a_r::R
    b_r::T
    function GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dx::T, order = 1) where {T}
        nl = length(αl)
        nr = length(αr)
        S_l = zeros(T, (nl-2, order+nl-2))
        S_r = zeros(T, (nr-2, order+nr-2))

        for i in 1:(nl-2)
            S_l[i,:] = [transpose(calculate_weights(i, one(T), Array(one(T):convert(T, order+i)))) transpose(zeros(T, Int(nl-2-i)))]./(dx^i) #am unsure if the length of the dummy_x is correct here
        end

        for i in 1:(nr-2)
            S_r[i,:] = [transpose(calculate_weights(i, convert(T, order+i), Array(one(T):convert(T, order+i)))) transpose(zeros(T, Int(nr-2-i)))]./(dx^i)
        end
        s0_l = S_l[:,1] ; Sl = S_l[:,2:end]
        s0_r = S_r[:,end] ; Sr = S_r[:,(end-1):-1:1]

        denoml = αl[2] .+ αl[3:end] ⋅ s0_l
        denomr = αr[2] .+ αr[3:end] ⋅ s0_r

        a_l = -transpose(transpose(αl[3:end]) * Sl) ./denoml
        a_r = reverse(-transpose(transpose(αr[3:end]) * Sr) ./denomr)

        b_l = -αl[1]/denoml
        b_r = -αr[1]/denomr
        new{T, typeof(a_l), typeof(a_r)}(a_l,b_l,a_r,b_r)
    end

    function GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dx::AbstractVector{T}, order = 1) where {T}

        nl = length(αl)
        nr = length(αr)
        dx_l, dx_r = (dx[1:(order+nl-2)], reverse(dx[(end-order-nr+3):end]))
        S_l = zeros(T, (nl-2, order+nl-2))
        S_r = zeros(T, (nr-2, order+nr-2))

        for i in 1:(nl-2)
            S_l[i,:] = [transpose(calculate_weights(i, one(T), Array(one(T):convert(T, order+i)))) transpose(zeros(T, Int(nl-2-i)))]./(dx_l.^i)
        end

        for i in 1:(nr-2)
            S_r[i,:] = [transpose(calculate_weights(i, convert(T, order+i), Array(one(T):convert(T, order+i)))) transpose(zeros(T, Int(nr-2-i)))]./(dx_r.^i)
        end
        s0_l = S_l[:,1] ; Sl = S_l[:,2:end]
        s0_r = S_r[:,end] ; Sr = S_r[:,(end-1):-1:1]

        denoml = αl[2] .+ αl[3:end] ⋅ s0_l
        denomr = αr[2] .+ αr[3:end] ⋅ s0_r

        a_l = -transpose(transpose(αl[3:end]) * Sl) ./denoml
        a_r = reverse(-transpose(transpose(αr[3:end]) * Sr) ./denomr)

        b_l = -αl[1]/denoml
        b_r = -αr[1]/denomr
        new{T,  typeof(a_l), typeof(a_r)}(a_l,b_l,a_r,b_r)
    end
end

#implement Neumann and Dirichlet as special cases of RobinBC
NeumannBC(α::NTuple{2,T}, dx::Union{AbstractVector{T}, T}, order = 1) where T = RobinBC((zero(T), one(T), α[1]), (zero(T), one(T), α[2]), dx, order)
DirichletBC(αl::T, αr::T) where T = RobinBC((one(T), zero(T), αl), (one(T), zero(T), αr), 1.0, 2.0 )
#specialized constructors for Neumann0 and Dirichlet0
Dirichlet0BC(T::Type) = DirichletBC(zero(T), zero(T))
Neumann0BC(dx::Union{AbstractVector{T}, T}, order = 1) where T = NeumannBC((zero(T), zero(T)), dx, order)

# other acceptable argument signatures
#RobinBC(al::T, bl::T, cl::T, dx_l::T, ar::T, br::T, cr::T, dx_r::T, order = 1) where T = RobinBC([al,bl, cl], [ar, br, cr], dx_l, order)

Base.:*(Q::AffineBC, u::AbstractVector) = BoundaryPaddedVector(Q.a_l ⋅ u[1:length(Q.a_l)] + Q.b_l, Q.a_r ⋅ u[(end-length(Q.a_r)+1):end] + Q.b_r, u)
Base.:*(Q::PeriodicBC, u::AbstractVector) = BoundaryPaddedVector(u[end], u[1], u)

Base.size(Q::AtomicBC) = (Inf, Inf) #Is this nessecary?

gettype(Q::AbstractBC{T}) where T = T