text:
"""
Runs the Carleo & Troyer (2017) algorithm.
"""

function main(NetSettings)
    @unpack mfree, mc_trials, dim_rbm, n, repetitions, calc_stat, γ_init, γ_decay,
    save_figures, modelname, use_meter = NetSettings

    HamiltWeights = init_weights(WEIGHTS, modelname, n)
    NetParams = init_params(NETPARAMS, NetSettings)
    CGParams = init_cg_params(CGPARAMS, NetSettings)
    VarDervs = init_vardervs(VARDERIVS, NetSettings)

    if mfree
        OptParams = init_optim_mfree(OPTPARAMSMFREE, NetSettings)
    else
        OptParams = init_optim(OPTPARAMS, NetSettings)
    end

    elocs = SharedArray(zeros(Complex{Float64}, mc_trials, nworkers()))
    # Stores the energy samples of each worker in a column.
    forces_w = SharedArray(zeros(Complex{Float64}, dim_rbm, nworkers()))
    # Stores the forces of each worker in a column.
    warm_states = zeros(Int, n, nworkers())
    # Saves last visited state in a chain.
    mean_e = zeros(repetitions)
    γ = γ_init

    progress_gui = Progress(repetitions)  # Progress bar
    for i=1:repetitions
        @sync for (id, worker) in enumerate(workers())
            @async warm_states[:, id] = remotecall_fetch(mc_sampling!,
                worker,
                NetParams,
                VarDervs.o_tot,
                elocs,
                HamiltWeights,
                NetSettings,
                id,
                warm_state=warm_states[:, id]
                )
            end

        mean_e[i] = mean(real(elocs))
        cg_iters = optimize!(CGParams, VarDervs, OptParams, elocs, forces_w, γ, NetSettings)
        update_netparams!(NetParams, OptParams, NetSettings)
        γ *= γ_decay
        if use_meter
            ProgressMeter.next!(progress_gui; showvalues = [(:iter,i),
                (:E,mean_e[i]), (:CG_iters, cg_iters)])
        end
    end
    @everywhere GC.gc()
    eloc_and_stats = [0.0, init_stat(STATISTICS, n)]
    if calc_stat
        eloc_and_stats = stats_calculator(NetParams, HamiltWeights, warm_states, NetSettings)
    end
    return mean_e, NetParams, eloc_and_stats[1], eloc_and_stats[2]
end

"""
    run_NQS(NetSettings)

# Arguments
    - `NetSettings`: The settings struct containing all possible network settings.
"""
function run_NQS(NetSettings)
    @unpack modelname, n, test_iter, calc_stat, α, mag0, writetofile, save_figures = NetSettings
    init_print(NetSettings)

    energy = zeros(Float64, 0)
    for i=1:test_iter
        println("RUN $i")
        energy, NetParams, final_eloc, GroundStats = main(NetSettings)

        if writetofile
            if calc_stat
                namebase = "$modelname-n=$n-mag0=$mag0-alpha=$α-run=$i-"
                for k in fieldnames(STATISTICS)
                    writedlm(namebase * "statistics-$k.txt", getfield(GroundStats, k))
                end
                writedlm(namebase * "mc_and_final_eloc.txt", [energy; final_eloc])
                if save_figures
                    plot_statistics(GroundStats, namebase, n)
                end
            else
                writedlm(namebase * "mc_eloc.txt", energy)
            end
            for k in fieldnames(NETPARAMS)
                writedlm(namebase * "netparams-$k.txt", getfield(NetParams, k))
            end
        end
    end
    return energy
end

function init_print(NetSettings)
    @unpack n, modelname, pbc, γ_init, γ_decay, regulator, mag0, α, repetitions,
        mc_trials, therm_steps, stat_samples, mfree, iterative_inverse,
        writetofile =  NetSettings
    println("## Initializing $n spin system
## model = $modelname
## pbc = $pbc
## learning rate = $γ_init
## lrate-decay = $γ_decay
## S-matrix diag shift = $regulator
## workers = $(nworkers())
## ∑s_z=0-sector sampling = $mag0
## α = $α
## iterations = $repetitions
## samples per worker = $mc_trials
## thermalization steps = $therm_steps
## final samples per worker = $stat_samples
## matrixfree = $mfree
## conjugate gradient inversion = $iterative_inverse
## saving results = $writetofile")
end


text:
using Documenter, MultivariateBases

makedocs(
    sitename = "MultivariateBases",
    # See https://github.com/JuliaDocs/Documenter.jl/issues/868
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    # See https://github.com/JuliaOpt/JuMP.jl/issues/1576
    strict = true,
    pages = [
        "Introduction" => "index.md",
    ],
    # The following ensures that we only include the docstrings from
    # this module for functions define in Base that we overwrite.
    modules = [MultivariateBases]
)

deploydocs(
    repo   = "github.com/JuliaAlgebra/MultivariateBases.jl.git",
)

ext: jl

max_stars_repo_name: JuliaAlgebra/PolynomialBases

max_issues_repo_name: JuliaAlgebra/PolynomialBases

max_forks_repo_name: JuliaAlgebra/PolynomialBases


text:
module CodeCosts

import InteractiveUtils: gen_call_with_extracted_types_and_kwargs

export @code_costs

struct CostsSummary
    cost_threshold::Int
    zero::Int              # cost == 0
    cheap::Int             # cost == 1
    middle::Vector{Int}    # 2 <= cost <= 5
    expensive::Vector{Int} # 6 <= cost
end

function CostsSummary(cost_threshold::Int, costs::Vector{Int})
    zero = 0
    cheap = 0
    middle = []
    expensive = []
    for c in costs
        if c == 0
            zero += 1
        elseif c == 1
            cheap += 1
        elseif c <= 5
            push!(middle, c)
        else
            push!(expensive, c)
        end
    end
    CostsSummary(cost_threshold, zero, cheap,
                 sort!(middle, rev=true), sort!(expensive, rev=true))
end

struct Hint{ID} end

struct CodeCostsInfo
    ci::Core.CodeInfo
    costs::Vector{Int}
    summary::CostsSummary
    hints::Vector{Tuple{Int,Hint}}
end

macro code_costs(ex0...)
    thecall = gen_call_with_extracted_types_and_kwargs(__module__, :code_costs, ex0)
    quote
        results = $thecall
        results isa CodeCostsInfo ? results : results[1]
    end
end

if isdefined(Core.Compiler, :OptimizationParams)
    function code_costs(f, types; debuginfo::Symbol=:default,
                        params::Core.Compiler.OptimizationParams=Core.Compiler.OptimizationParams(),
                        hints::Symbol=:default)

        mi = Base.method_instances(f, types)[1]
        ci = code_typed(f, types; debuginfo=debuginfo)[1][1]

        interp = Core.Compiler.NativeInterpreter()
        opt = Core.Compiler.OptimizationState(mi, params, interp)

        cost(statement) = 0
        function cost(statement::Expr)
            Core.Compiler.statement_cost(statement, -1, ci, opt.sptypes, opt.slottypes, opt.params)
        end
        raw_costs = map(cost, ci.code)
        cost_threshold = opt.params.inline_cost_threshold
        summary = CostsSummary(cost_threshold, costs(raw_costs, ci))
        hint_messages = Vector{Tuple{Int,Hint}}()
        CodeCostsInfo(ci, raw_costs, summary, hint_messages)
    end
else
    function code_costs(f, types; debuginfo::Symbol=:default,
                        params::Core.Compiler.Params=Core.Compiler.Params(typemax(UInt)),
                        hints::Symbol=:default)

        mi = Base.method_instances(f, types)[1]
        if VERSION >= v"1.1"
            ci = code_typed(f, types; debuginfo=debuginfo)[1][1]
        else
            ci = code_typed(f, types)[1][1]
        end

        opt = Core.Compiler.OptimizationState(mi, params)
        opt.src.inlineable = true
        sptypes = :sptypes in fieldnames(typeof(opt)) ? opt.sptypes : opt.sp

        cost(statement) = 0
        function cost(statement::Expr)
            Core.Compiler.statement_cost(statement, -1, ci, sptypes, opt.slottypes, opt.params)
        end
        raw_costs = map(cost, ci.code)
        cost_threshold = opt.params.inline_cost_threshold
        summary = CostsSummary(cost_threshold, costs(raw_costs, ci))
        hint_messages = Vector{Tuple{Int,Hint}}()
        CodeCostsInfo(ci, raw_costs, summary, hint_messages)
    end
end

function _code_hash(s::String)
    h = UInt(0)
    for c in codeunits(replace(s, r"\e\[[^m]*m" => ""))
        h = (0x20 < c) & (c < 0x80) ? hash(c, h) : h
    end
    h
end

function Base.show(io::IO, costsinfo::CodeCostsInfo; debuginfo::Symbol=:source)
    buf_s = IOBuffer()
    buf_n = IOBuffer()
    if VERSION >= v"1.1"
        show(IOContext(buf_s, io), costsinfo.ci, debuginfo=debuginfo)
        debuginfo !== :none && show(buf_n, costsinfo.ci, debuginfo=:none)
    else
        show(IOContext(buf_s, io), costsinfo.ci)
    end
    println(io, "CodeCostsInfo(")

    seekstart(buf_s)
    seekstart(buf_n)
    idx = 0
    nhash = _code_hash(readline(buf_n))
    for line in readlines(buf_s)
        l = VERSION >= v"1.1" ? line : replace(line, "\e[1G" => "\e[6G")

        if VERSION >= v"1.1" && debuginfo !== :none
            shash = _code_hash(l)
            if shash !== nhash
                println(io, " "^5, l)
                continue
            end
            n = readline(buf_n)
            nhash = _code_hash(n)
        end

        if 0 < idx <= length(costsinfo.costs)
            # multi-line statement
            if !occursin(r"(?:^|\e\[6G)(?:\e\[[^m]*m|\s)*(?:\d+ |│|└)", l)
                println(io, " "^5, l)
                continue
            end
            # skip code_coverage_effect
            if Meta.isexpr(costsinfo.ci.code[idx], :code_coverage_effect)
                idx += 1
                continue
            end

            c = costsinfo.costs[idx]
            cs = lpad(string(c), 4)
            if c == 1
                println(io, cs, " ", l)
            else
                color = c == 0 ? Base.info_color() :
                        c <= 5 ? Base.warn_color() : Base.error_color()
                printstyled(io, cs, color=color)
                println(io, " ", l)
            end
        else
            println(io, " "^5, l)
        end
        idx += 1
    end
    print(io, ", ")
    print(io, costsinfo.summary)
    print(io, ")")
end

function Base.show(io::IO, s::CostsSummary)
    println(io, "CodeCostsSummary(")
    mids, exps = sum(s.middle), sum(s.expensive)
    total = s.cheap + mids + exps
    pad(n) = lpad(string(n), max(ndigits(total), ndigits(s.zero)))
    printstyled(io, "     zero: ", pad(s.zero),  "|", color=Base.info_color())
    println(io)
    print(      io, "    cheap: ", pad(s.cheap), "| ", costs_string(fill(1, s.cheap)))
    println(io)
    printstyled(io, "   middle: ", pad(mids),    "| ", costs_string(s.middle), color=Base.warn_color())
    println(io)
    printstyled(io, "expensive: ", pad(exps),    "| ", costs_string(s.expensive), color=Base.error_color())
    println(io)
    print(      io, "    total: ", pad(total),   "| ", s.cost_threshold, " (default threshold)")
    println(io)
    print(io, ")")
end

function costs_string(costs::Vector{Int})
    buf = IOBuffer()
    len = 0
    for c in costs
        cs = string(c)
        n = min(50, c - length(cs))
        print(buf, cs, "="^n)
        len += c
        if len >= 75
            print(buf, "...")
            break
        end
    end
    String(take!(buf))
end

function costs(raw_costs::Vector{Int}, ci::Core.CodeInfo)
    cs = Int[]
    for (i, c) in enumerate(raw_costs)
        Meta.isexpr(ci.code[i], :code_coverage_effect) || push!(cs, c)
    end
    cs
end
costs(costsinfo::CodeCosts.CodeCostsInfo) = costs(costsinfo.costs, costsinfo.ci)

end # module


text:
using Random
Random.seed!(7)

dt = 0.01
integration = :rk4

###################
## Parallel Park ##
###################

opts = SolverOptions()
opts.verbose = false
opts.cost_tolerance = 1e-5
opts.cost_tolerance_intermediate = 1e-5
opts.constraint_tolerance = 1e-5
opts.resolve_feasible = true
opts.outer_loop_update_type = :default
opts.R_infeasible = 10

# Set up model, objective, and solver
model, obj = TrajectoryOptimization.Dynamics.dubinscar_parallelpark
solver = Solver(model, obj, integration=integration, N=101, opts=opts)
U0 = ones(solver.model.m,solver.N-1)
X0 = line_trajectory(solver)

results, stats = TrajectoryOptimization.solve(solver,U0)
# results_inf, stats_inf = TrajectoryOptimization.solve(solver,X0,U0)

max_violation(results)
results_new = copy(results)

# newton_solve!(results_new,solver)

newton_results = NewtonResults(solver)
newton_active_set!(newton_results,results_new,solver)
# sum(newton_results.active_set)
# sum(vcat(results_new.active_set...))
# sum(newton_results.active_set_ineq)
# newton_results.s[findall(x->x != 0.0, newton_results.s)]
# findall(x->x < 0.0, vcat(results.C...)[newton_results.active_set_ineq])
update_newton_results!(newton_results,results_new,solver)
newton_step!(results_new,newton_results,solver,1.0,0.00000001)
max_violation(results_new)


text:
using Base.Test

@testset "AlphaGo" begin
  include("test_go.jl")
  include("test_mcts.jl")
  include("test_mcts_player.jl")
  include("test_features.jl")
end


text:
type StateSpace
  timeseries::Array{Float64, 1}
  embedding::Array{Float64, 2}
  embeddingdim::Int64
  embeddinglag::Int64
  triangulation_vertices::Array{Float64, 2}
  triangulation_indices::Array{Int64, 2}
  indices_simplices::Array{Int64} # Indices of simplices forming the triangulation of the state space
  indices_images::Array{Int64}   # Indices of the images of the simplices under the map
  P::Array{Float64, 2}   # Transition probability matrix
end

"""
 Extract a simplex from a triangulated state space instance at the index specified by 'index'
"""
function getsimplex(statespace, index)
  statespace.embedding[statespace.indices_simplices[index, :], :].'
end

max_stars_repo_name: kahaaga/InvariantDistribution

max_issues_repo_name: kahaaga/InvariantDistribution

max_forks_repo_name: kahaaga/InvariantDistribution


text:
using MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities

# indices
const VI = MOI.VariableIndex
const CI = MOI.ConstraintIndex
# supported functions
const SAF = MOI.ScalarAffineFunction{Float64}
const SQF = MOI.ScalarQuadraticFunction{Float64}
const VAF = MOI.VectorAffineFunction{Float64}
const VECTOR = MOI.VectorOfVariables
# supported sets
const BOUNDS = Union{MOI.EqualTo{Float64}, MOI.GreaterThan{Float64},
                     MOI.LessThan{Float64}, MOI.Interval{Float64}}
const VAR_TYPES = Union{MOI.ZeroOne, MOI.Integer}
const SOS1 = MOI.SOS1{Float64}
const SOS2 = MOI.SOS2{Float64}
# other MOI types
const AFF_TERM = MOI.ScalarAffineTerm{Float64}
const QUAD_TERM = MOI.ScalarQuadraticTerm{Float64}
const VEC_TERM = MOI.VectorAffineTerm{Float64}

const PtrMap = Dict{Ptr{Cvoid}, Union{VarRef, ConsRef}}
const ConsTypeMap = Dict{Tuple{DataType, DataType}, Set{ConsRef}}

mutable struct Optimizer <: MOI.AbstractOptimizer
    inner::SCIPData
    reference::PtrMap
    constypes::ConsTypeMap
    binbounds::Dict{VI,BOUNDS} # only for binary variables
    params::Dict{String,Any}
    start::Dict{VI,Float64} # can be partial
    moi_separator::Any # ::Union{CutCbSeparator, Nothing}
    objective_sense::MOI.OptimizationSense

    function Optimizer(; kwargs...)
        scip = Ref{Ptr{SCIP_}}(C_NULL)
        @SCIP_CALL SCIPcreate(scip)
        @assert scip[] != C_NULL
        @SCIP_CALL SCIPincludeDefaultPlugins(scip[])
        @SCIP_CALL SCIP.SCIPcreateProbBasic(scip[], "")

        scip_data = SCIPData(scip, Dict(), Dict(), 0, 0, Dict(), Dict(), Dict(), [])

        o = new(scip_data, PtrMap(), ConsTypeMap(), Dict(), Dict(), Dict(), nothing, MOI.MIN_SENSE)
        finalizer(free_scip, o)

        # Set all parameters given as keyword arguments, replacing the
        # delimiter, since "/" is used by all SCIP parameters, but is not
        # allowed in Julia identifiers.
        for (key, value) in kwargs
            name = replace(String(key),"_" => "/")
            MOI.set(o, MOI.RawOptimizerAttribute(name), value)
        end
        return o
    end
end

free_scip(o::Optimizer) = free_scip(o.inner)

Base.cconvert(::Type{Ptr{SCIP_}}, o::Optimizer) = o
# Protect Optimizer from GC for ccall with Ptr{SCIP_} argument.
Base.unsafe_convert(::Type{Ptr{SCIP_}}, o::Optimizer) = o.inner.scip[]

## convenience functions (not part of MOI)

"Return pointer to SCIP variable."
function var(o::Optimizer, v::VI)::Ptr{SCIP_VAR}
    return var(o.inner, VarRef(v.value))
end

"Return var/cons reference of SCIP variable/constraint."
ref(o::Optimizer, ptr::Ptr{Cvoid}) = o.reference[ptr]

"Return pointer to SCIP constraint."
function cons(o::Optimizer, c::CI{F,S})::Ptr{SCIP_CONS} where {F,S}
    return cons(o.inner, ConsRef(c.value))
end

"Extract bounds from sets."
bounds(set::MOI.EqualTo{Float64}) = (set.value, set.value)
bounds(set::MOI.GreaterThan{Float64}) = (set.lower, nothing)
bounds(set::MOI.LessThan{Float64}) = (nothing, set.upper)
bounds(set::MOI.Interval{Float64}) = (set.lower, set.upper)

"Make set from bounds."
from_bounds(::Type{MOI.EqualTo{Float64}}, lower, upper) = MOI.EqualTo{Float64}(lower)
from_bounds(::Type{MOI.GreaterThan{Float64}}, lower, upper) = MOI.GreaterThan{Float64}(lower)
from_bounds(::Type{MOI.LessThan{Float64}}, lower, upper) = MOI.LessThan{Float64}(upper)
from_bounds(::Type{MOI.Interval{Float64}}, lower, upper) = MOI.Interval{Float64}(lower, upper)

"Register pointer in mapping, return var/cons reference."
function register!(o::Optimizer, ptr::Ptr{Cvoid}, ref::R) where R <: Union{VarRef, ConsRef}
    @assert !haskey(o.reference, ptr)
    o.reference[ptr] = ref
    return ref
end

"Register constraint in mapping, return constraint reference."
function register!(o::Optimizer, c::CI{F,S}) where {F,S}
    cr = ConsRef(c.value)
    if haskey(o.constypes, (F, S))
        push!(o.constypes[F,S], cr)
    else
        o.constypes[F,S] = Set([cr])
    end
    return c
end

"Go back from solved stage to problem modification stage, invalidating results."
function allow_modification(o::Optimizer)
    if !(SCIPgetStage(o) in (SCIP_STAGE_PROBLEM, SCIP_STAGE_SOLVING))
        @SCIP_CALL SCIPfreeTransform(o)
    end
    return nothing
end

## general queries and support

MOI.get(::Optimizer, ::MOI.SolverName) = "SCIP"

MOI.supports_incremental_interface(::Optimizer) = true

function _throw_if_invalid(o::Optimizer, ci::CI{F, S}) where {F, S}
    if !haskey(o.constypes, (F, S)) || !in(ConsRef(ci.value), o.constypes[F, S])
        throw(MOI.InvalidIndex(ci))
    end
    return nothing
end

function MOI.get(o::Optimizer, param::MOI.RawOptimizerAttribute)
    return get_parameter(o.inner, param.name)
end

function MOI.set(o::Optimizer, param::MOI.RawOptimizerAttribute, value)
    set_parameter(o.inner, param.name, value)
    o.params[param.name] = value
    return nothing
end

MOI.supports(o::Optimizer, ::MOI.Silent) = true

function MOI.get(o::Optimizer, ::MOI.Silent)
    return MOI.get(o, MOI.RawOptimizerAttribute("display/verblevel")) == 0
end

function MOI.set(o::Optimizer, ::MOI.Silent, value)
    param = MOI.RawOptimizerAttribute("display/verblevel")
    if value
        MOI.set(o, param, 0) # no output at all
    else
        MOI.set(o, param, 4) # default level
    end
end

MOI.supports(o::Optimizer, ::MOI.TimeLimitSec) = true

function MOI.get(o::Optimizer, ::MOI.TimeLimitSec)
    raw_value = MOI.get(o, MOI.RawOptimizerAttribute("limits/time"))
    if raw_value == SCIPinfinity(o)
        return nothing
    else
        return raw_value
    end
end

function MOI.set(o::Optimizer, ::MOI.TimeLimitSec, value)
    if value === nothing
        MOI.set(o, MOI.RawOptimizerAttribute("limits/time"), SCIPinfinity(o))
    else
        MOI.set(o, MOI.RawOptimizerAttribute("limits/time"), value)
    end
end

MOI.supports(::Optimizer, ::MOI.SolverVersion) = true

MOI.get(::Optimizer, ::MOI.SolverVersion) = "v" * string(SCIP_versionnumber())

## model creation, query and modification

function MOI.is_empty(o::Optimizer)
    return length(o.inner.vars) == 0 && length(o.inner.conss) == 0
end

function MOI.empty!(o::Optimizer)
    # free the underlying problem
    free_scip(o.inner)
    # clear auxiliary mapping structures
    o.reference = PtrMap()
    o.constypes = ConsTypeMap()
    o.binbounds = Dict()
    o.start = Dict()
    # manually recreate empty o.inner (formerly done by creating a new mscip before ManagedSCIP was removed)
    scip = Ref{Ptr{SCIP_}}(C_NULL)
    @SCIP_CALL SCIPcreate(scip)
    @assert scip[] != C_NULL
    @SCIP_CALL SCIPincludeDefaultPlugins(scip[])
    @SCIP_CALL SCIP.SCIPcreateProbBasic(scip[], "")
    # create a new problem
    o.inner = SCIPData(scip, Dict(), Dict(), 0, 0, Dict(), Dict(), Dict(), [])
    # reapply parameters
    for pair in o.params
        set_parameter(o.inner, pair.first, pair.second)
    end
    return nothing
end

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

MOI.get(o::Optimizer, ::MOI.Name) = unsafe_string(SCIPgetProbName(o))
MOI.set(o::Optimizer, ::MOI.Name, name::String) = @SCIP_CALL SCIPsetProbName(o, name)

function MOI.get(o::Optimizer, ::MOI.NumberOfConstraints{F,S}) where {F,S}
    return haskey(o.constypes, (F, S)) ? length(o.constypes[F, S]) : 0
end

function MOI.get(o::Optimizer, ::MOI.ListOfConstraintTypesPresent)
    return collect(keys(o.constypes))
end

function MOI.get(o::Optimizer, ::MOI.ListOfConstraintIndices{F, S}) where {F, S}
    list_indices = Vector{CI{F,S}}()
    if !haskey(o.constypes, (F, S))
        return list_indices
    end
    for cref in o.constypes[F, S]
        push!(list_indices, CI{F,S}(cref.val))
    end
    return sort!(list_indices, by=v->v.value)
end

function set_start_values(o::Optimizer)
    if isempty(o.start)
        # no primal start values are given
        return
    end

    # create new partial solution object
    sol__ = Ref{Ptr{SCIP_SOL}}(C_NULL)
    @SCIP_CALL SCIPcreatePartialSol(o, sol__, C_NULL)
    @assert sol__[] != C_NULL

    # set all given values
    sol_ = sol__[]
    for (vi, value) in o.start
        @SCIP_CALL SCIPsetSolVal(o, sol_, var(o, vi), value)
    end

    # submit the candidate
    stored_ = Ref{SCIP_Bool}(FALSE)
    @SCIP_CALL SCIPaddSolFree(o, sol__, stored_)
    @assert sol__[] == C_NULL
end

function MOI.optimize!(o::Optimizer)
    set_start_values(o)
    if o.objective_sense == MOI.FEASIBILITY_SENSE
        MOI.set(o, MOI.ObjectiveFunction{SAF}(), SAF([], 0.0))
    end
    @SCIP_CALL SCIPsolve(o)
    return nothing
end

function MOI.delete(o::Optimizer, ci::CI{F, S}) where {F, S}
    _throw_if_invalid(o, ci)
    allow_modification(o)
    delete!(o.constypes[F, S], ConsRef(ci.value))
    if isempty(o.constypes[F, S])
        delete!(o.constypes, (F, S))
    end
    delete!(o.reference, cons(o, ci))
    delete(o.inner, ConsRef(ci.value))
    return nothing
end

function MOI.get(o::Optimizer, ::MOI.ListOfVariableAttributesSet)
    attributes = MOI.AbstractVariableAttribute[MOI.VariableName()]
    if !isempty(o.start)
        push!(attributes, MOI.VariablePrimalStart())
    end
    return attributes
end

MOI.get(::Optimizer, ::MOI.ListOfModelAttributesSet) = MOI.AbstractModelAttribute[MOI.Name(), MOI.ObjectiveSense()]

function MOI.get(
    ::Optimizer,
    ::MOI.ListOfConstraintAttributesSet{F,S},
) where {F,S}
    attributes = MOI.AbstractConstraintAttribute[]
    if F != MOI.VariableIndex
        return push!(attributes, MOI.ConstraintName())
    end
    return attributes
end

function MOI.get(::Optimizer, ::MOI.ListOfOptimizerAttributesSet)
    attributes = MOI.ListOfOptimizerAttributesSet[]
    timelim = MOI.get(o, MOI.TimeLimitSec())
    if timelim !== nothing
        push!(attributes, MOI.TimeLimitSec())
    end
    return attributes
end

include(joinpath("MOI_wrapper", "variable.jl"))
include(joinpath("MOI_wrapper", "constraints.jl"))
include(joinpath("MOI_wrapper", "linear_constraints.jl"))
include(joinpath("MOI_wrapper", "quadratic_constraints.jl"))
include(joinpath("MOI_wrapper", "sos_constraints.jl"))
include(joinpath("MOI_wrapper", "indicator_constraints.jl"))
include(joinpath("MOI_wrapper", "nonlinear_constraints.jl"))
include(joinpath("MOI_wrapper", "objective.jl"))
include(joinpath("MOI_wrapper", "results.jl"))
include(joinpath("MOI_wrapper", "conshdlr.jl"))
include(joinpath("MOI_wrapper", "sepa.jl"))


text:
using Gadfly

set_default_plot_size(6inch, 3inch)

plot(layer(y=[4.72067, 6.93034, 6.25414, 0.158186, 9.77226, 4.21936, 19.0383, 19.9981, 5.03324, 19.7333, 11.115, 8.74216, 8.49436, 15.4645, 5.6238, 4.18945, 5.02758, 0.407497, 5.75403, 17.1902],
           Geom.bar),
     layer(x=[1.53902, 12.8079, 17.4709, 5.57165, 15.0263, 12.8977, 1.55653, 16.9637, 1.7127, 11.0641, 9.267, 3.71643, 2.23962, 19.5262, 1.03229, 10.7606, 9.11384, 5.5879, 3.56492, 10.9797],
           Geom.line))


text:
export NEAT_populate!, NEAT_evaluate!, GA_NEAT_populate!, NEAT_tournament, NEAT_random_top

function NEAT_tournament(pop::Array{NEATIndiv}, t_size::Int)
    if length(pop) == 1
        return pop[1]
    end
    inds = shuffle!(collect(1:length(pop)))
    sort(pop[inds[1:t_size]])[end]
end

function NEAT_random_top(pop::Array{NEATIndiv}, threshold::Float64)
    sort!(pop, by = x -> x.fitness, rev = true)
    max_index = Integer(ceil(threshold * length(pop)))
    pop[rand(1:max_index)]
end

function NEAT_populate!(e::Evolution, selection::Function)
    # Create offsprings and add to the population
    new_pop::Array{NEATIndiv} = []
    if e.cfg["elites"]>0
        sort!(e.population, by = x -> x.fitness, rev = true)
        for i in 1:e.cfg["elites"]
            push!(new_pop, e.population[i])
        end
    end
    # println("\n-------\nTotal fitness: ", e.cfg["total_fitness"], " Pop: ", length(e.population))
    for s in values(e.cfg["Species"])
        # println("\nFitness: ", s.fitness_val, " Members: ", length(s.members))
        if e.cfg["total_fitness"] > 0
            nb_offspring = Integer(round(
                e.cfg["n_population"] * s.fitness_val / e.cfg["total_fitness"],
            ))
        else
            nb_offspring = length(s.members) #Same species size
        end
        # println("Children: ", nb_offspring)
        reproduction!(s, selection, nb_offspring, e.cfg)
        append!(new_pop, s.members)
    end

    e.population = new_pop
    if e.cfg["verbose"]
        println(
            "\n",
            e.gen,
            " - Pop: ",
            length(e.population),
            "  Species: ",
            length(e.cfg["Species"]),
            " Inno: ",
            e.cfg["innovation_max"],
        )
        println("Best fitness: ", maximum(getfield.(e.population, :fitness))[1])
    end
end

function NEAT_evaluate!(e::Evolution, fitness::Function)
    # Build
    build!.(e.population)
    # Assign species
    if length(e.cfg["Species"]) == 0
        s = Species(e.cfg)
    end
    for i in e.population
        find_species!(i, e.cfg)
    end

    # update distance threshold to keep species number stable
    update_threshold!(e.cfg)

    # Update species dict
    for i in keys(e.cfg["Species"])
        if length(e.cfg["Species"][i].members) == 0
            pop!(e.cfg["Species"], i)
        end
    end

    # Compute fitness
    e.cfg["total_fitness"] = 0
    if e.cfg["use_max_fitness"]
        for s in values(e.cfg["Species"])
            e.cfg["total_fitness"] += compute_fitness_max!(s, fitness)
        end
    else
        for s in values(e.cfg["Species"])
            e.cfg["total_fitness"] += compute_fitness_mean!(s, fitness)
        end
    end
end

function update_threshold!(cfg::Dict)
    if length(cfg["Species"]) == 0
        println("NO SPECIES")
        return nothing
    end
    species_pop = []
    for s in values(cfg["Species"])
        push!(species_pop, length(s.members))
    end
    if cfg["verbose"]
        if length(species_pop) > 1
            sort!(species_pop)
            println(
                "Indiv/species | Mean: ",
                sum(species_pop) / length(species_pop),
                " | Med: ",
                species_pop[Integer(floor(length(species_pop) / 2))],
            )
        else
            println("Only 1 species: ", sum(species_pop), " elements")
        end
    end

    nb_species = length(cfg["Species"])
    nb_excess = nb_species - cfg["target_species_number"]
    dist_mod = cfg["dist_mod"] * abs(nb_excess) / cfg["target_species_number"]
    if nb_excess < 0
        cfg["dist_threshold"] -= dist_mod
        if cfg["verbose"]
            println(
                "Dist Threshold - ",
                dist_mod,
                " -> ",
                cfg["dist_threshold"],
            )
        end
    elseif nb_excess > 0
        cfg["dist_threshold"] += dist_mod
        if cfg["verbose"]
            println(
                "Dist Threshold + ",
                dist_mod,
                " -> ",
                cfg["dist_threshold"],
            )
        end
    end

    if cfg["dist_threshold"] < cfg["dist_mod"]
        cfg["dist_threshold"] = cfg["dist_mod"]
    end
end

function GA_NEAT_populate!(e; args...)
    Cambrian.ga_populate(
        e,
        mutation = NeuroEvolution.mutate,
        crossover = NeuroEvolution.crossover,
        selection = x ->
            tournament_selection(x, e.cfg["tournament_size"]),
    )
end


text:
#derivative logic

#ForwardDiff compiles one method per function, so separating the
#differentiation logic from the property logic allows the differentials
#to be compiled only once

"""
    ∂f∂T(model,V,T,z=SA[1.0])

returns `f` and `∂f/∂T` at constant total volume and composition, where f is the total helmholtz energy, given by `eos(model,V,T,z)`

"""
function ∂f∂T(model,V,T,z=SA[1.0])
    f(∂T) = eos(model,V,∂T,z)
    return Solvers.derivative(f,T)
end

"""
    ∂f∂V(model,V,T,z=SA[1.0])

returns `f` and `∂f/∂V` at constant temperature and composition, where f is the total helmholtz energy, given by `eos(model,V,T,z)`, and V is the total volume

"""
function ∂f∂V(model,V,T,z)
    f(∂V) = eos(model,∂V,T,z)
    return Solvers.derivative(f,V)
end

#returns a tuple of the form ([∂f∂V,∂f∂T],f),using the least amount of computation
"""
    ∂f(model,V,T,z)

returns zeroth order (value) and first order derivative information of the total helmholtz energy (given by `eos(model,V,T,z)`).
the result is given in two values:

```julia
grad_f,fval = ∂2f(model,V,T,z)
```

where:

```julia
fval   = f(V,T) = eos(model,V,T,z)

grad_f = [ ∂f/∂V; ∂f/∂T]

```

Where `V` is the total volume, `T` is the temperature and `f` is the total helmholtz energy.
"""
function ∂f(model,V,T,z)
    f(∂V,∂T) = eos(model,∂V,∂T,z)
    _f,_df = Solvers.fgradf2(f,V,T)
    return _df,_f
end

#returns p and ∂p∂V at constant T
#it doesnt do a pass over temperature, so its
#faster that d2f when only requiring d2fdV2

"""
    p∂p∂V(model,V,T,z=SA[1.0])

returns `p` and `∂p/∂V` at constant temperature, where p is the pressure = `pressure(model,V,T,z)` and `V` is the total Volume.

"""
function p∂p∂V(model,V,T,z=SA[1.0])
    f(∂V) = pressure(model,∂V,T,z)
    p,∂p∂V = Solvers.f∂f(f,V)
    return p,∂p∂V
end

"""
    ∂2f(model,V,T,z)

returns zeroth order (value), first order and second order derivative information of the total helmholtz energy (given by `eos(model,V,T,z)`).
the result is given in three values:

```
hess_f,grad_f,fval = ∂2f(model,V,T,z)
```

where:
```
fval   = f(V,T) = eos(model,V,T,z)

grad_f = [ ∂f/∂V; ∂f/∂T]

hess_f = [ ∂²f/∂V²; ∂²f/∂V∂T
          ∂²f/∂V∂T; ∂²f/∂V²]
 ```

Where `V` is the total volume, `T` is the temperature and `f` is the total helmholtz energy.
"""
function ∂2f(model,V,T,z)
    f(_V,_T) = eos(model,_V,_T,z)
    _f,_∂f,_∂2f =  Solvers.∂2(f,V,T)
    return (_∂2f,_∂f,_f)
end

"""
    ∂2p(model,V,T,z)

returns zeroth order (value), first order and second order derivative information of the pressure.
the result is given in three values:

```
hess_p,grad_p,pval = ∂2p(model,V,T,z)
```

where:
```
pval   = p(V,T) = pressure(model,V,T,z)

grad_p = [ ∂p/∂V; ∂p/∂T]

hess_p = [ ∂²p/∂V²; ∂²p/∂V∂T
          ∂²p/∂V∂T; ∂²p/∂V²]
 ```

Where `V` is the total volume, `T` is the temperature and `p` is the pressure.
"""
function ∂2p(model,V,T,z)
    f(_V,_T) = pressure(model,_V,_T,z)
    _f,_∂f,_∂2f =  Solvers.∂2(f,V,T)
    return (_∂2f,_∂f,_f)
end

"""
    f_hess(model,V,T,z)

returns the second order volume (`V`) and temperature (`T`) derivatives of the total helmholtz energy (given by `eos(model,V,T,z)`). the result is given in a 2x2 `SMatrix`, in the form:

```
[ ∂²f/∂V²  ∂²f/∂V∂T
 ∂²f/∂V∂T  ∂²f/∂T²]
 ```

use this instead of the ∂2f if you only need second order information. ∂2f also gives zeroth and first order derivative information, but due to a bug in the used AD, it allocates more than necessary.
"""
function f_hess(model,V,T,z)
    f(w) = eos(model,first(w),last(w),z)
    V,T = promote(V,T)
    VT_vec = SVector(V,T)
    return ForwardDiff.hessian(f,VT_vec)
end

"""
    ∂²³f(model,V,T,z=SA[1.0])

returns `∂²A/∂V²` and `∂³A/∂V³`, in a single ForwardDiff pass. used mainly in `crit_pure` objective function

"""
function ∂²³f(model,V,T,z=SA[1.0])
    p(∂²A∂V²) =  pressure(model,∂²A∂V²,T,z)
    f(∂A∂V) = Solvers.derivative(p,∂A∂V)
    ∂²A∂V², ∂³A∂V³ = Solvers.f∂f(f,V)
    return ∂²A∂V², ∂³A∂V³
end


text:
"""
# Description

Given an expression, `e`, and `Set{ColumnName}`, `column_names`, this
function will recursively descend the expression adding any symbols to
`column_names` that should be considered column names based on their syntactic
position in the expression. The Boolean flag is `true` if the symbol is a
dynamic name rather than a static name.

See docs/expression_rules.md for a high-level description of the logic of how
potential column names are meant to be identified.

# Arguments

1. `e::Any`: An expression.
2. `column_names::Set{ColumnName}`: A set of symbols that will be
    mutated to contain all column names found in `e`.

# Return Values

1. `column_names::Set{ColumnName}`: The input `column_names`, which
    has been mutated to contain the column names.

# Examples

```
julia>  _find_column_names!(:(a + sin(b)), Set{ColumnName}())
```
"""
function _find_column_names!(
    @nospecialize(e::Any),
    column_names::Set{ColumnName},
)::Nothing
    if isa(e, Expr) && e.head == :(=)
        for i in 2:length(e.args)
            _find_column_names!(e.args[i], column_names)
        end
        nothing
    elseif isa(e, Expr) && e.head == :kw
        _find_column_names!(e.args[2], column_names)
        nothing
    elseif isa(e, Expr) && e.head == :call
        j = isa(e.args[1], Symbol) ? 2 : 1
        for i in j:length(e.args)
            _find_column_names!(e.args[i], column_names)
        end
        nothing
    elseif isa(e, Expr) && e.head == :$
        nothing
    elseif isa(e, Expr) && e.head == :comparison
        for i in 1:length(e.args)
            if i % 2 == 1
                _find_column_names!(e.args[i], column_names)
            end
        end
        nothing
    elseif isa(e, Expr) && e.head == :macrocall && isa(e.args[1], GlobalRef) && e.args[1].name == Symbol("@cmd")
        if startswith(e.args[3], '$')
            push!(column_names, ColumnName(Symbol(strip(e.args[3], '$')), true))
        else
            push!(column_names, ColumnName(Symbol(e.args[3]), false))
        end
        nothing
    elseif isa(e, Expr)
        for i in 1:length(e.args)
            _find_column_names!(e.args[i], column_names)
        end
        nothing
    elseif isa(e, Symbol)
        push!(column_names, ColumnName(e, false))
        nothing
    else
        nothing
    end
end

"""
# Description

Find all column names occuring in an expression, `e`, and return them as a
tuple of symbols.

# Arguments

1. `e::Any`: An expression.

# Return Values

1. `column_names::NTuple{N, ColumnName} where N`: A tuple containing
    all column names found in `e` and a flag indicating whether they are
    static or dynamic names.

# Examples
```
julia> find_column_names(:(a + b + sin(c) + `\$d`))
((:a, false), (:d, true), (:c, false), (:b, false))
```
"""
function find_column_names(@nospecialize(e::Any))::Tuple
    column_names = Set{ColumnName}()
    _find_column_names!(e, column_names)
    (column_names..., )
end


text:
#############################################################
#
#  add extra quaternion operations for Quaternions
#  some should probably be merged back into Quaternions
#
#############################################################

#=
"""
    mean(qvec, [method = 0])

Function to find the mean of a set of rotations.

Args:

    qvec    - Vector of `Rotation`s

    method  - if 0 (default), the returned rotation minimizes sum_i=1:n (||R(Qi) - R(Qbar)||)  with || the Frobenius norm, or
              equivalently minimizes sum_i=1:n (sin(phi_i / 2))^2, where phi_i is the angle of Qe_i, with  Qbar = Q_i x Qe_i

            - if 1, the returned rotation minimizes sum_i=1:n (phi_i ^2) where phi_i is the angle of Qe_i, with  Qbar = Q_i x Qe_i
              this should be about the same as method == 0 if phi_i are all small (note: still not implemented)
"""
=#

"""
    mean(r::AbstractVector{Rotation})  ->  r̄

Find the mean `r̄` of a set of rotations `r`. The returned rotation minimizes `sum_i=1:n (||r[i] - r̄||)`  with `||` the Frobenius norm, or
equivalently minimizes `sum_i=1:n (sin(ϕ[i] / 2))^2`, where `ϕ[i] = rotation_angle(r[i] / r̄)`.
"""
function Statistics.mean(qvec::AbstractVector{QuatRotation{T}}, method::Integer = 0) where T
    #if (method == 0)
        M = zeros(4, 4)
        for i in 1:length(qvec)
            q = qvec[i]
            Qi = @SVector [q.w, q.x, q.y, q.z]  # convert types to ensure we don't get quaternion multiplication
            M .+= Qi * (Qi')
        end
        evec = eigfact(Symmetric(M), 4:4)
        Qbar = QuatRotation(evec.vectors[1], evec.vectors[2], evec.vectors[3], evec.vectors[4]) # This will renormalize the quaternion...
    #else
    #    error("I haven't coded this")
    #end

    return Qbar
end

function Statistics.mean(vec::AbstractVector{R}) where R<:Rotation
    R(mean(convert(Vector{QuatRotation{eltype(R)}}, vec)))
end


text:
"""
    LdagLSparseOpProblem <: AbstractProblem

Problem or finding the steady state of a ℒdagℒ matrix
"""
struct LdagLSparseOpProblem{B, SM} <: HermitianMatrixProblem where {B<:Basis,
                                                 SM<:SparseMatrixCSC}
    HilbSpace::B            # 0
    H::SM                   # 1
    H_t::SM                 # 2
    HdH::SM                 # 3
    L_ops::Vector{SM}       # 4
    L_ops_t::Vector{SM}     # 5

    LdH_ops::Vector{SM}     # 6
    HdL_ops::Vector{SM}     # 7
    LdL_ops_t::Matrix{SM}   # 8
    ρss
end

basis(prob::LdagLSparseOpProblem) = prob.HilbSpace

"""
    LdagLSparseOpProblem([T=STD_REAL_PREC], lindbladian)

Creates a problem for minimizing the cost function 𝒞 = ∑|ρ(σ)|²|⟨⟨σ|ℒ'ℒ |ρ⟩⟩|².
Computes |⟨⟨σ|ℒ'ℒ |ρ⟩⟩| by computing on the fly commutators with the
Hamiltonian and with the collapse operators.

`lindbladian` can either be the lindbladian on a graph, a QuantumOptics superoperator
or the Hamiltonian and a vector of collapse operators.

`T=STD_REAL_PREC` by default is the numerical precision used. It should match that of
the network.
"""
LdagLSparseOpProblem(args...) = LdagLSparseOpProblem(STD_REAL_PREC, args...)
LdagLSparseOpProblem(T::Type{<:Number}, gl::GraphLindbladian) =
    LdagLSparseOpProblem(T, basis(gl), SparseOperator(hamiltonian(gl)), jump_operators(gl))
LdagLSparseOpProblem(T::Type{<:Number}, ham::DataOperator, cops) =
    LdagLSparseOpProblem(T, ham.basis_l, ham, cops)
function LdagLSparseOpProblem(T::Type{<:Number}, Hilb::Basis, Ham::DataOperator, c_ops_q)
    # Fix complex numbers
    if real(T) == T
        T = Complex{T}
    end
    # Generate H_eff
    H_eff   = deepcopy(T.(Ham.data))
    ST = typeof(H_eff)

    c_ops       = Vector{ST}(undef, length(c_ops_q))
    c_ops_trans = Vector{ST}(undef, length(c_ops_q))
    for i=1:length(c_ops)
        c_ops[i]       = c_ops_q[i].data
        c_ops_trans[i] = transpose(c_ops[i])
        H_eff         -= 0.5im * (c_ops[i]'*c_ops[i])
    end

    LdH_ops     = Vector{ST}(undef, length(c_ops))
    HdL_ops     = Vector{ST}(undef, length(c_ops))
    for i=1:length(c_ops)
        LdH_ops[i]     = c_ops[i]'*H_eff
        HdL_ops[i]     = LdH_ops[i]'
    end

    LdL_ops_t   = Matrix{ST}(undef, length(c_ops), length(c_ops))
    for i=1:length(c_ops)
        for j=1:length(c_ops)
            LdL_ops_t[i,j] = transpose(c_ops[i]'*c_ops[j])
        end
    end

    LdagLSparseOpProblem{typeof(Hilb), ST}(
                 Hilb,                  # 0
                 H_eff,                 # 1
                 transpose(H_eff),
                 H_eff'*H_eff,          # 3
                 c_ops,
                 c_ops_trans,           # 5
                 LdH_ops,
                 HdL_ops,               # 7
                 LdL_ops_t,             # 8
                 0.0)
end

function compute_Cloc(prob::LdagLSparseOpProblem, net::MatrixNet, 𝝝, lnψ=net(𝝝), 𝝝p=deepcopy(𝝝))
    # Quantities of the problem
    H = prob.H
    H_t = prob.H_t
    HdH = prob.HdH
    c_ops = prob.L_ops
    c_ops_trans = prob.L_ops_t
    LdH_ops = prob.LdH_ops
    HdL_ops = prob.HdL_ops
    LdL_ops_t = prob.LdL_ops_t

    σ  = row(𝝝)
    σt = col(𝝝)
    set_index!(𝝝p, index(𝝝))
    𝝝p_row = row(𝝝p)
    𝝝p_col = col(𝝝p)

    i_σt = index(σt)
    i_σ  = index(σ)

    C_loc = zero(Complex{real(out_type(net))})

    # ⟨σ|ρHᴴH|σt⟩
    set!(𝝝p_row, toint(σ))
    for row_id = HdH.colptr[i_σt]:(HdH.colptr[i_σt+1]-1)
      i_σ_p = HdH.rowval[row_id]
      set_index!(𝝝p_col, i_σ_p)

      log_ratio = net(𝝝p) - lnψ
      ΔE  =  HdH.nzval[row_id] * exp(log_ratio) #mat[i_σ, i_σ_p]
      C_loc  += ΔE #2.0*real(ΔE)
    end

    # ⟨σ|HᴴHρ|σt⟩ (using hermitianity of HdH)
    set!(𝝝p_col, toint(σt))
    for row_id = HdH.colptr[i_σ]:(HdH.colptr[i_σ+1]-1)
      i_σ_p = HdH.rowval[row_id]
      set_index!(𝝝p_row, i_σ_p)

      log_ratio = logψ(net, 𝝝p) - lnψ
      ΔE  =  conj(HdH.nzval[row_id]) * exp(log_ratio) #mat[i_σ_p, i_σt]
      C_loc  += ΔE #2.0*real(ΔE)
    end

    # ⟨σ|HᴴρHᴴ|σt⟩
    for L_row_id = H.colptr[i_σ]:(H.colptr[i_σ+1]-1)
      i_σ_p = H.rowval[L_row_id]
      L_σp_σ = H.nzval[L_row_id]
      set_index!(𝝝p_row, i_σ_p)

      for Ldag_row_id = H_t.colptr[i_σt]:(H_t.colptr[i_σt+1]-1)
        # Find nonzero elements s by doing <i_σp|L|i_σ> and backconvert
        i_σt_p = H_t.rowval[Ldag_row_id]
        set_index!(𝝝p_col, i_σt_p)

        log_ratio = logψ(net, 𝝝p)- lnψ
        #@assert (mat[i_σ_p, i_σt_p]/exp(lnψ) - log_ratio) == 0
        ΔE  =  - conj(L_σp_σ) * conj(H_t.nzval[Ldag_row_id]) * exp(log_ratio) # mat[i_σ_p, i_σt_p]
        C_loc  += ΔE #2.0*real(ΔE)
      end
    end

    # ⟨σ|HρH|σt⟩
    for L_row_id = H_t.colptr[i_σ]:(H_t.colptr[i_σ+1]-1)
      i_σ_p = H_t.rowval[L_row_id]
      L_σp_σ = H_t.nzval[L_row_id]
      set_index!(𝝝p_row, i_σ_p)

      for Ldag_row_id = H.colptr[i_σt]:(H.colptr[i_σt+1]-1)
        # Find nonzero elements s by doing <i_σp|L|i_σ> and backconvert
        i_σt_p = H.rowval[Ldag_row_id]
        set_index!(𝝝p_col, i_σt_p)

        log_ratio = logψ(net, 𝝝p) - lnψ

        ΔE  =  - L_σp_σ * H.nzval[Ldag_row_id] * exp(log_ratio) # mat[i_σ_p, i_σt_p]
        C_loc  += ΔE #2.0*real(ΔE)
      end
    end

    # L rho Ldag H #ok
    # -im ⟨σ|L ρ LᴴH|σt⟩
    for i=1:length(LdH_ops)
      LdH = LdH_ops[i]
      L   = c_ops_trans[i]

      for ext_row_id = L.colptr[i_σ]:(L.colptr[i_σ+1]-1)
        i_σ_p = L.rowval[ext_row_id]
        set_index!(𝝝p_row, i_σ_p)
        val_σ_p = L.nzval[ext_row_id]

        for int_row_id = LdH.colptr[i_σt]:(LdH.colptr[i_σt+1]-1)
          i_σt_p = LdH.rowval[int_row_id]
          set_index!(𝝝p_col, i_σt_p)

          log_ratio = logψ(net, 𝝝p) - lnψ
          #@assert (mat[i_σ_p, i_σt_p] - log_ratio) == 0
          ΔE  = -1.0im * val_σ_p * LdH.nzval[int_row_id] *  exp(log_ratio) # mat[i_σ_p, i_σt_p]
          C_loc  += ΔE # 2.0*real(ΔE)
        end
      end
    end

    # +im ⟨σ|HᴴL ρ Lᴴ|σt⟩
    for i=1:length(LdH_ops)
      LdH = LdH_ops[i]
      L   = c_ops_trans[i]

      for ext_row_id = L.colptr[i_σt]:(L.colptr[i_σt+1]-1)
        i_σ_p = L.rowval[ext_row_id]
        set_index!(𝝝p_col, i_σ_p)
        val_σ_p = L.nzval[ext_row_id]

        for int_row_id = LdH.colptr[i_σ]:(LdH.colptr[i_σ+1]-1)
          i_σt_p = LdH.rowval[int_row_id]
          set_index!(𝝝p_row, i_σt_p)

          log_ratio = logψ(net, 𝝝p) - lnψ
          #@assert (mat[i_σt_p, i_σ_p] - log_ratio) == 0

          ΔE  = 1.0im * conj(val_σ_p) * conj(LdH.nzval[int_row_id]) *  exp(log_ratio) # mat[i_σt_p, i_σ_p]
          C_loc  += ΔE # 2.0*real(ΔE)
        end
      end
    end

    # Ldag rho Hdag L #ok
    # ⟨σ|Lᴴ ρ HᴴL|σt⟩
    for i=1:length(LdH_ops)
      HdL = HdL_ops[i]
      L   = c_ops[i]

      for ext_row_id = L.colptr[i_σ]:(L.colptr[i_σ+1]-1)
        i_σ_p = L.rowval[ext_row_id]
        set_index!(𝝝p_row, i_σ_p)
        val_σ_p = L.nzval[ext_row_id]

        for int_row_id = HdL.colptr[i_σt]:(HdL.colptr[i_σt+1]-1)
          i_σt_p = HdL.rowval[int_row_id]
          set_index!(𝝝p_col, i_σt_p)

          log_ratio = logψ(net, 𝝝p) - lnψ
          #@assert (mat[i_σ_p, i_σt_p] - log_ratio) == 0

          ΔE  = 1.0im * conj(val_σ_p) * HdL.nzval[int_row_id] * exp(log_ratio) # mat[i_σ_p, i_σt_p]
          C_loc  += ΔE # 2.0*real(ΔE)
        end
      end
    end

    # ⟨σ|Lᴴ H ρ L|σt⟩
    for i=1:length(LdH_ops)
      HdL = HdL_ops[i]
      L   = c_ops[i]

      for ext_row_id = L.colptr[i_σt]:(L.colptr[i_σt+1]-1)
        i_σ_p = L.rowval[ext_row_id]
        set_index!(𝝝p_col, i_σ_p)
        val_σ_p = L.nzval[ext_row_id]

        for int_row_id = HdL.colptr[i_σ]:(HdL.colptr[i_σ+1]-1)
          i_σt_p = HdL.rowval[int_row_id]
          set_index!(𝝝p_row, i_σt_p)

          log_ratio = logψ(net, 𝝝p) - lnψ
          #@assert (mat[i_σt_p, i_σ_p] - log_ratio) == 0

          ΔE  = -1.0im * val_σ_p * conj(HdL.nzval[int_row_id]) * exp(log_ratio) # mat[i_σt_p, i_σ_p]
          C_loc  += ΔE # 2.0*real(ΔE)
        end
      end
    end

    for LdL=LdL_ops_t
      for ext_row_id = LdL.colptr[i_σ]:(LdL.colptr[i_σ+1]-1)
        i_σ_p = LdL.rowval[ext_row_id]
        set_index!(𝝝p_row, i_σ_p)
        val_σ_p = LdL.nzval[ext_row_id]
        for int_row_id = LdL.colptr[i_σt]:(LdL.colptr[i_σt+1]-1)
          i_σt_p = LdL.rowval[int_row_id]
          set_index!(𝝝p_col, i_σt_p)

          log_ratio = logψ(net, 𝝝p) - lnψ
          #@assert (mat[i_σ_p, i_σt_p] - log_ratio) == 0

          ΔE  = val_σ_p * conj(LdL.nzval[int_row_id]) * exp(log_ratio) #mat[i_σ_p, i_σt_p]
          C_loc  += ΔE
        end
      end
    end

    return C_loc
end

Base.show(io::IO, p::LdagLSparseOpProblem) = print(io,
    "LdagLSparseOpProblem on space : $(basis(p)), computing the energy of LdagL with H, jump_ops")


text:
"""
# KM3NeT Data Definitions v1.3.1-8-gae7dfb2
https://git.km3net.de/common/km3net-dataformat
"""

module W2list_genhen
  const W2LIST_GENHEN_GLOBAL_GEN_WEIGHT = 0
  const W2LIST_GENHEN_EG = 1
  const W2LIST_GENHEN_SIG = 2
  const W2LIST_GENHEN_COLUMN_DEPTH = 3
  const W2LIST_GENHEN_P_EARTH = 4
  const W2LIST_GENHEN_REFF = 5
  const W2LIST_GENHEN_BX = 7
  const W2LIST_GENHEN_BY = 8
  const W2LIST_GENHEN_ICHAN = 9
  const W2LIST_GENHEN_CC = 10
end


text:
using Statistics

# Read in file of helper functions
# This file incluse: cleantext, countwords
include("../word_count_helpers.jl")

# Load the file names
dataLoc = "../../../data/word_count/";
fnames = dataLoc.*readdir(dataLoc)

# Create function for processing each file
function getwordcounts(fname)
    # Read in file and clean the text
    f = open(fname)
    text=cleantext(readlines(f))

    # Count number of times each word appears
    return countwords(text)
end

# Find word counts for all books
allcounts = getwordcounts.(fnames)

# Get overall counts for each word
overallcounts = merge(+,allcounts...)

# Iterate through each set of word counts
for counts in allcounts

    # Calculate normalized score
    thresh = quantile(collect(values(counts)), .98)
    counts = filter!(word-> word.second > thresh, counts)
    normscore = merge(/,counts,filter(word->haskey(counts,word.first), overallcounts))

    # Sort and print the top 5 words with their normalized counts
    rankedwords = sort(collect(normscore), by=x->x[2], rev=true)
    println(rankedwords[1:min(5,length(counts)),:])
end

max_stars_repo_name: llsc-supercloud/teaching-examples

max_issues_repo_name: llsc-supercloud/teaching-examples

max_forks_repo_name: llsc-supercloud/teaching-examples


text:
struct GeneralUnstructuredMesh_2D <: UnstructuredMesh_2D
    name::String
    points::Vector{Point_2D}
    edges::Vector{<:SVector{L, UInt32} where {L}}
    materialized_edges::Vector{<:Edge_2D}
    faces::Vector{<:SArray{S, UInt32, 1, L} where {S<:Tuple, L}}
    materialized_faces::Vector{<:Face_2D}
    edge_face_connectivity::Vector{SVector{2, UInt32}}
    face_edge_connectivity::Vector{<:SArray{S, UInt32, 1, L} where {S<:Tuple, L}}
    boundary_edges::Vector{Vector{UInt32}}
    face_sets::Dict{String, Set{UInt32}}
end

function GeneralUnstructuredMesh_2D(;
        name::String = "DefaultMeshName",
        points::Vector{Point_2D} = Point_2D[],
        edges::Vector{<:SVector{L, UInt32} where {L}} = SVector{2, UInt32}[],
        materialized_edges::Vector{<:Edge_2D} = LineSegment_2D[],
        faces::Vector{<:SArray{S, UInt32, 1, L} where {S<:Tuple, L}} = SVector{4, UInt32}[],
        materialized_faces::Vector{<:Face_2D} = Triangle_2D[],
        edge_face_connectivity::Vector{SVector{2, UInt32}} = SVector{2, UInt32}[],
        face_edge_connectivity ::Vector{<:SArray{S, UInt32, 1, L} where {S<:Tuple, L}} = SVector{3, UInt32}[],
        boundary_edges::Vector{Vector{UInt32}} = Vector{UInt32}[],
        face_sets::Dict{String, Set{UInt32}} = Dict{String, Set{UInt32}}()
    ) where {F <: AbstractFloat, U <: Unsigned}
        return GeneralUnstructuredMesh_2D(name,
                                          points,
                                          edges,
                                          materialized_edges,
                                          faces,
                                          materialized_faces,
                                          edge_face_connectivity,
                                          face_edge_connectivity,
                                          boundary_edges,
                                          face_sets,
                                         )
end

Base.broadcastable(mesh::GeneralUnstructuredMesh_2D) = Ref(mesh)

# Return a mesh with boundary edges
function add_boundary_edges(mesh::M; boundary_shape="Unknown"
    ) where {M <: UnstructuredMesh_2D}
    if 0 === length(mesh.edge_face_connectivity)
        mesh = add_connectivity(mesh)
    end
    return M(name = mesh.name,
             points = mesh.points,
             edges = mesh.edges,
             materialized_edges = mesh.materialized_edges,
             faces = mesh.faces,
             materialized_faces = mesh.materialized_faces,
             edge_face_connectivity = mesh.edge_face_connectivity,
             face_edge_connectivity = mesh.face_edge_connectivity,
             boundary_edges = boundary_edges(mesh, boundary_shape),
             face_sets = mesh.face_sets
            )
end

# Return a mesh with face/edge connectivity and edge/face connectivity
function add_connectivity(mesh::UnstructuredMesh_2D)
    return add_edge_face_connectivity(mesh)
end

# Return a mesh with its edges
function add_edges(mesh::M) where {M <: UnstructuredMesh_2D}
    return M(name = mesh.name,
             points = mesh.points,
             edges = edges(mesh),
             materialized_edges = mesh.materialized_edges,
             faces = mesh.faces,
             materialized_faces = mesh.materialized_faces,
             edge_face_connectivity = mesh.edge_face_connectivity,
             face_edge_connectivity = mesh.face_edge_connectivity,
             boundary_edges = mesh.boundary_edges,
             face_sets = mesh.face_sets
            )
end

# Return a mesh with edge/face connectivity
function add_edge_face_connectivity(mesh::M) where {M <: UnstructuredMesh_2D}
    if 0 === length(mesh.face_edge_connectivity)
        mesh = add_face_edge_connectivity(mesh)
    end
    return M(name = mesh.name,
             points = mesh.points,
             edges = mesh.edges,
             materialized_edges = mesh.materialized_edges,
             faces = mesh.faces,
             materialized_faces = mesh.materialized_faces,
             edge_face_connectivity = edge_face_connectivity(mesh),
             face_edge_connectivity = mesh.face_edge_connectivity,
             boundary_edges = mesh.boundary_edges,
             face_sets = mesh.face_sets
            )
end

# Return a mesh with every field created
function add_everything(mesh::UnstructuredMesh_2D)
    return add_materialized_faces(
             add_materialized_edges(
               add_boundary_edges(mesh, boundary_shape = "Rectangle")))
end

# Return a mesh with face/edge connectivity
function add_face_edge_connectivity(mesh::M) where {M <: UnstructuredMesh_2D}
    if 0 === length(mesh.edges)
        mesh = add_edges(mesh)
    end
    return M(name = mesh.name,
             points = mesh.points,
             edges = mesh.edges,
             materialized_edges = mesh.materialized_edges,
             faces = mesh.faces,
             materialized_faces = mesh.materialized_faces,
             edge_face_connectivity = mesh.edge_face_connectivity,
             face_edge_connectivity = face_edge_connectivity(mesh),
             boundary_edges = mesh.boundary_edges,
             face_sets = mesh.face_sets
            )
end

# Return a mesh with materialized edges
function add_materialized_edges(mesh::M) where {M <: UnstructuredMesh_2D}
    if 0 === length(mesh.edges)
        mesh = add_edges(mesh)
    end
    return M(name = mesh.name,
             points = mesh.points,
             edges = mesh.edges,
             materialized_edges = materialize_edges(mesh),
             faces = mesh.faces,
             materialized_faces = mesh.materialized_faces,
             edge_face_connectivity = mesh.edge_face_connectivity,
             face_edge_connectivity = mesh.face_edge_connectivity,
             boundary_edges = mesh.boundary_edges,
             face_sets = mesh.face_sets
            )
end

# Return a mesh with materialized faces
function add_materialized_faces(mesh::M) where {M <: UnstructuredMesh_2D}
    return M(name = mesh.name,
             points = mesh.points,
             edges = mesh.edges,
             materialized_edges = mesh.materialized_edges,
             faces = mesh.faces,
             materialized_faces = materialize_faces(mesh),
             edge_face_connectivity = mesh.edge_face_connectivity,
             face_edge_connectivity = mesh.face_edge_connectivity,
             boundary_edges = mesh.boundary_edges,
             face_sets = mesh.face_sets
            )
end

# Return a vector of the faces adjacent to the face of ID face
function adjacent_faces(face::UInt32, mesh::UnstructuredMesh_2D)
    edges = mesh.face_edge_connectivity[face]
    the_adjacent_faces = UInt32[]
    for edge in edges
        faces = mesh.edge_face_connectivity[edge]
        for face_id in faces
            if face_id != face && face_id != 0
                push!(the_adjacent_faces, face_id)
            end
        end
    end
    return the_adjacent_faces
end

# Area of face
function area(face::SVector{N, UInt32}, points::Vector{Point_2D}) where {N}
    return area(materialize_face(face, points))
end

# Return the area of a face set
function area(mesh::UnstructuredMesh_2D, face_set::Set{UInt32})
    if 0 < length(mesh.materialized_faces)
        return mapreduce(x->area(mesh.materialized_faces[x]), +, face_set)
    else
        return mapreduce(x->area(mesh.faces[x], mesh.points), +, face_set)
    end
end

# Return the area of a face set by name
function area(mesh::UnstructuredMesh_2D, set_name::String)
    return area(mesh, mesh.face_sets[set_name])
end

# Bounding box of a vector of points
function boundingbox(points::Vector{Point_2D})
    x = getindex.(points, 1)
    y = getindex.(points, 2)
    return Rectangle_2D(minimum(x), minimum(y), maximum(x), maximum(y))
end

# Bounding box of a vector of points
function boundingbox(points::SVector{L, Point_2D}) where {L}
    x = getindex.(points, 1)
    y = getindex.(points, 2)
    return Rectangle_2D(minimum(x), minimum(y), maximum(x), maximum(y))
end

# Return a vector containing vectors of the edges in each side of the mesh's bounding shape, e.g.
# For a rectangular bounding shape the sides are North, East, South, West. Then the output would
# be [ [e1, e2, e3, ...], [e17, e18, e18, ...], ..., [e100, e101, ...]]
function boundary_edges(mesh::UnstructuredMesh_2D, boundary_shape::String)
    # edges which have face 0 in their edge_face connectivity are boundary edges
    the_boundary_edges = UInt32.(findall(x->x[1] === 0x00000000, mesh.edge_face_connectivity))
    if boundary_shape == "Rectangle"
        # Sort edges into NESW
        bb = boundingbox(mesh.points)
        y_north = bb.ymax
        x_east  = bb.xmax
        y_south = bb.ymin
        x_west  = bb.xmin
        p_NW = Point_2D(x_west, y_north)
        p_NE = bb.tr
        p_SE = Point_2D(x_east, y_south)
        p_SW = bb.bl
        edges_north = UInt32[]
        edges_east = UInt32[]
        edges_south = UInt32[]
        edges_west = UInt32[]
        # Insert edges so that indices move from NW -> NE -> SE -> SW -> NW
        for edge ∈  the_boundary_edges
            epoints = edgepoints(mesh.edges[edge], mesh.points)
            if all(x->abs(x[2] - y_north) < 1e-4, epoints)
                insert_boundary_edge!(edge, p_NW, edges_north, mesh)
            elseif all(x->abs(x[1] - x_east) < 1e-4, epoints)
                insert_boundary_edge!(edge, p_NE, edges_east, mesh)
            elseif all(x->abs(x[2] - y_south) < 1e-4, epoints)
                insert_boundary_edge!(edge, p_SE, edges_south, mesh)
            elseif all(x->abs(x[1] - x_west) < 1e-4, epoints)
                insert_boundary_edge!(edge, p_SW, edges_west, mesh)
            else
                @error "Edge $edge could not be classified as NSEW"
            end
        end
        return [ edges_north, edges_east, edges_south, edges_west ]
    else
        return [ convert(Vector{UInt32}, the_boundary_edges) ]
    end
end

# SVector of MVectors of point IDs representing the 3 edges of a triangle
function edges(face::SVector{3, UInt32})
    edges = SVector( MVector{2, UInt32}(face[1], face[2]),
                     MVector{2, UInt32}(face[2], face[3]),
                     MVector{2, UInt32}(face[3], face[1]) )
    # Order the linear edge vertices by ID
    for edge in edges
        if edge[2] < edge[1]
            e1 = edge[1]
            edge[1] = edge[2]
            edge[2] = e1
        end
    end
    return edges
end

# SVector of MVectors of point IDs representing the 4 edges of a quadrilateral
function edges(face::SVector{4, UInt32})
    edges = SVector( MVector{2, UInt32}(face[1], face[2]),
                     MVector{2, UInt32}(face[2], face[3]),
                     MVector{2, UInt32}(face[3], face[4]),
                     MVector{2, UInt32}(face[4], face[1]) )
    # Order the linear edge vertices by ID
    for edge in edges
        if edge[2] < edge[1]
            e1 = edge[1]
            edge[1] = edge[2]
            edge[2] = e1
        end
    end
    return edges
end

# SVector of MVectors of point IDs representing the 3 edges of a quadratic triangle
function edges(face::SVector{6, UInt32})
    edges = SVector( MVector{3, UInt32}(face[1], face[2], face[4]),
                     MVector{3, UInt32}(face[2], face[3], face[5]),
                     MVector{3, UInt32}(face[3], face[1], face[6]) )
    # Order the linear edge vertices by ID
    for edge in edges
        if edge[2] < edge[1]
            e1 = edge[1]
            edge[1] = edge[2]
            edge[2] = e1
        end
    end
    return edges
end

# SVector of MVectors of point IDs representing the 4 edges of a quadratic quadrilateral
function edges(face::SVector{8, UInt32})
    edges = SVector( MVector{3, UInt32}(face[1], face[2], face[5]),
                     MVector{3, UInt32}(face[2], face[3], face[6]),
                     MVector{3, UInt32}(face[3], face[4], face[7]),
                     MVector{3, UInt32}(face[4], face[1], face[8]) )
    # Order th linear edge vertices by ID
    for edge in edges
        if edge[2] < edge[1]
            e1 = edge[1]
            edge[1] = edge[2]
            edge[2] = e1
        end
    end
    return edges
end

# The unique edges from a vector of triangles or quadrilaterals represented by point IDs
function edges(mesh::UnstructuredMesh_2D)
    edges_filtered = sort(unique(reduce(vcat, edges.(mesh.faces))))
    return [ SVector(e.data) for e in edges_filtered ]
end

# A vector of length 2 SVectors, denoting the face ID each edge is connected to. If the edge
# is a boundary edge, face ID 0 is returned
function edge_face_connectivity(mesh::UnstructuredMesh_2D)
    # Each edge should only border 2 faces if it is an interior edge, and 1 face if it is
    # a boundary edge.
    # Loop through each face in the face_edge_connectivity vector and mark each edge with
    # the faces that it borders.
    if length(mesh.edges) === 0
        @error "Does not have edges!"
    end
    if length(mesh.face_edge_connectivity) === 0
        @error "Does not have face/edge connectivity!"
    end
    edge_face = [MVector{2, UInt32}(0, 0) for _ in eachindex(mesh.edges)]
    for (iface, edges) in enumerate(mesh.face_edge_connectivity)
        for iedge in edges
            # Add the face id in the first non-zero position of the edge_face conn. vec.
            if edge_face[iedge][1] == 0
                edge_face[iedge][1] = iface
            elseif edge_face[iedge][2] == 0
                edge_face[iedge][2] = iface
            else
                @error "Edge $iedge seems to have 3 faces associated with it!"
            end
        end
    end
    return [SVector(sort(two_faces).data) for two_faces in edge_face]
end

# Return an SVector of the points in the edge (Linear)
function edgepoints(edge::SVector{2, UInt32}, points::Vector{Point_2D})
    return SVector(points[edge[1]], points[edge[2]])
end

# Return an SVector of the points in the edge (Quadratic)
function edgepoints(edge::SVector{3, UInt32}, points::Vector{Point_2D})
    return SVector(points[edge[1]], points[edge[2]], points[edge[3]])
end

# Return an SVector of the points in the edge
function edgepoints(edge_id::UInt32, mesh::UnstructuredMesh_2D)
    return edgepoints(mesh.edges[edge_id], mesh.points)
end

# Return an SVector of the points in the face (Triangle)
function facepoints(face::SVector{3, UInt32}, points::Vector{Point_2D})
    return SVector(points[face[1]], points[face[2]], points[face[3]])
end

# Return an SVector of the points in the face (Quadrilateral)
function facepoints(face::SVector{4, UInt32}, points::Vector{Point_2D})
    return SVector(points[face[1]], points[face[2]], points[face[3]], points[face[4]])
end

# Return an SVector of the points in the face (Triangle6)
function facepoints(face::SVector{6, UInt32}, points::Vector{Point_2D})
    return SVector(points[face[1]], points[face[2]], points[face[3]],
                   points[face[4]], points[face[5]], points[face[6]])
end

# Return an SVector of the points in the face (Quadrilateral8)
function facepoints(face::SVector{8, UInt32}, points::Vector{Point_2D})
    return SVector(points[face[1]], points[face[2]], points[face[3]], points[face[4]],
                   points[face[5]], points[face[6]], points[face[7]], points[face[8]])
end

# Return an SVector of the points in the face
function facepoints(edge_id::UInt32, mesh::UnstructuredMesh_2D)
    return facepoints(mesh.faces[edge_id], mesh.points)
end

# Find the faces which share the vertex of ID v.
function faces_sharing_vertex(v::Integer, mesh::UnstructuredMesh_2D)
    shared_faces = UInt32[]
    for i ∈ 1:length(mesh.faces)
        if v ∈  mesh.faces[i]
            push!(shared_faces, UInt32(i))
        end
    end
    return shared_faces
end

# Return the face containing point p.
function findface(p::Point_2D, mesh::UnstructuredMesh_2D)
    if 0 < length(mesh.materialized_faces)
        return UInt32(findface_explicit(p, mesh.materialized_faces))
    else
        return UInt32(findface_implicit(p, mesh.faces, mesh.points))
    end
end

# Find the face containing the point p, with explicitly represented faces
function findface_explicit(p::Point_2D, faces::Vector{<:Face_2D})
    for i ∈ 1:length(faces)
        if p ∈  faces[i]
            return i
        end
    end
    return 0
end

# Return the face containing the point p, with implicitly represented faces
function findface_implicit(p::Point_2D,
                            faces::Vector{<:SArray{S, UInt32, 1, L} where {S<:Tuple, L}},
                            points::Vector{Point_2D})
    for i ∈ 1:length(faces)
        bool = p ∈  materialize_face(faces[i], points)
        if bool
            return i
        end
    end
    return 0
end

# Return the intersection algorithm that will be used for l ∩ mesh
function get_intersection_algorithm(mesh::UnstructuredMesh_2D)
    if length(mesh.materialized_edges) !== 0
        return "Edges - Explicit"
    elseif length(mesh.edges) !== 0
        return "Edges - Implicit"
    elseif length(mesh.materialized_faces) !== 0
        return "Faces - Explicit"
    else
        return "Faces - Implicit"
    end
end

# Insert the boundary edge into the correct place in the vector of edge indices, based on
# the distance from some reference point
function insert_boundary_edge!(edge_index::UInt32, p_ref::Point_2D, edge_indices::Vector{UInt32},
                               mesh::UnstructuredMesh_2D)
    # Compute the minimum distance from the edge to be inserted to the reference point
    insertion_distance = minimum(distance.(Ref(p_ref),
                                           edgepoints(mesh.edges[edge_index], mesh.points)))
    # Loop through the edge indices until an edge with greater distance from the reference point
    # is found, then insert
    nindices = length(edge_indices)
    for i ∈ 1:nindices
        epoints = edgepoints(mesh.edges[edge_indices[i]], mesh.points)
        edge_distance = minimum(distance.(Ref(p_ref), epoints))
        if insertion_distance < edge_distance
            insert!(edge_indices, i, edge_index)
            return nothing
        end
    end
    insert!(edge_indices, nindices+1, edge_index)
    return nothing
end

# Intersect a line with the mesh. Returns a vector of intersection points, sorted based
# upon distance from the line's start point
function intersect(l::LineSegment_2D, mesh::UnstructuredMesh_2D)
    # Edges are faster, so they are the default
    if length(mesh.edges) !== 0
        if 0 < length(mesh.materialized_edges)
            return intersect_edges_explicit(l, mesh.materialized_edges)
        else
            return intersect_edges_implicit(l, mesh.edges, mesh.points)
        end
    else
        if 0 < length(mesh.materialized_faces)
            return intersect_faces_explicit(l, mesh.materialized_faces)
        else
            return intersect_faces_implicit(l, mesh.faces, mesh.points)
        end
    end
end

# Intersect a line with an implicitly defined edge
function intersect_edge_implicit(l::LineSegment_2D,
                                 edge::SVector{L, UInt32} where {L},
                                 points::Vector{Point_2D})
    return l ∩ materialize_edge(edge, points)
end

# Intersect a line with materialized edges
function intersect_edges_explicit(l::LineSegment_2D, edges::Vector{LineSegment_2D})
    # A vector to hold all of the intersection points
    intersection_points = Point_2D[]
    for edge in edges
        npoints, point = l ∩ edge
        # If the intersections yields 1 or more points, push those points to the array of points
        if 0 < npoints
            push!(intersection_points, point)
        end
    end
    sort_intersection_points!(l[1], intersection_points)
    return intersection_points
end

# Intersect a line with implicitly defined edges
function intersect_edges_explicit(l::LineSegment_2D, edges::Vector{QuadraticSegment_2D})
    # A vector to hold all of the intersection points
    intersection_points = Point_2D[]
    for edge in edges
        npoints, points = l ∩ edge
        # If the intersections yields 1 or more points, push those points to the array of points
        if 0 < npoints
            append!(intersection_points, points[1:npoints])
        end
    end
    sort_intersection_points!(l[1], intersection_points)
    return intersection_points
end

# Intersect a line with a vector of implicitly defined linear edges
function intersect_edges_implicit(l::LineSegment_2D,
                                  edges::Vector{SVector{2, UInt32}},
                                  points::Vector{Point_2D})
    # A vector to hold all of the intersection points
    intersection_points = Point_2D[]
    # Intersect the line with each of the faces
    for edge in edges
        npoints, point = intersect_edge_implicit(l, edge, points)
        if 0 < npoints
            push!(intersection_points, point)
        end
    end
    sort_intersection_points!(l[1], intersection_points)
    return intersection_points
end

# Intersect a line with a vector of implicitly defined quadratic edges
function intersect_edges_implicit(l::LineSegment_2D,
                                  edges::Vector{SVector{3, UInt32}},
                                  points::Vector{Point_2D})
    # An array to hold all of the intersection points
    intersection_points = Point_2D[]
    # Intersect the line with each of the faces
    for edge in edges
        npoints, ipoints = intersect_edge_implicit(l, edge, points)
        if 0 < npoints
            append!(intersection_points, ipoints[1:npoints])
        end
    end
    sort_intersection_points!(l[1], intersection_points)
    return intersection_points
end

# Intersect a line with an implicitly defined face
function intersect_face_implicit(l::LineSegment_2D,
                                 face::SVector{L, UInt32} where {L},
                                 points::Vector{Point_2D})
    return l ∩ materialize_face(face, points)
end

# Intersect a line with explicitly defined linear faces
function intersect_faces_explicit(l::LineSegment_2D, faces::Vector{<:Face_2D} )
    # An array to hold all of the intersection points
    intersection_points = Point_2D[]
    for face in faces
        npoints, points = l ∩ face
        # If the intersections yields 1 or more points, push those points to the array of points
        if 0 < npoints
            append!(intersection_points, points[1:npoints])
        end
    end
    sort_intersection_points!(l[1], intersection_points)
    return intersection_points
end

# Intersect a line with implicitly defined faces
function intersect_faces_implicit(l::LineSegment_2D,
                                  faces::Vector{<:SArray{S, UInt32, 1, L} where {S<:Tuple, L}},
                                  points::Vector{Point_2D})
    # An array to hold all of the intersection points
    intersection_points = Point_2D[]
    # Intersect the line with each of the faces
    for face in faces
        npoints, ipoints = intersect_face_implicit(l, face, points)
        # If the intersections yields 1 or more points, push those points to the array of points
        if 0 < npoints
            append!(intersection_points, ipoints[1:npoints])
        end
    end
    sort_intersection_points!(l[1], intersection_points)
    return intersection_points
end

# If a point is a vertex
function isvertex(p::Point_2D, mesh::UnstructuredMesh_2D)
    for point in mesh.points
        if p ≈ point
            return true
        end
    end
    return false
end

# Return a LineSegment_2D from the point IDs in an edge
function materialize_edge(edge::SVector{2, UInt32}, points::Vector{Point_2D})
    return LineSegment_2D(edgepoints(edge, points))
end

# Return a QuadraticSegment_2D from the point IDs in an edge
function materialize_edge(edge::SVector{3, UInt32}, points::Vector{Point_2D})
    return QuadraticSegment_2D(edgepoints(edge, points))
end

# Return a LineSegment_2D or QuadraticSegment_2D
function materialize_edge(edge_id::UInt32, mesh::UnstructuredMesh_2D)
    return materialize_edge(mesh.edges[edge_id], mesh.points)
end

# Return a materialized edge for each edge in the mesh
function materialize_edges(mesh::UnstructuredMesh_2D)
    return materialize_edge.(mesh.edges, Ref(mesh.points))
end

# Return a Triangle_2D from the point IDs in a face
function materialize_face(face::SVector{3, UInt32}, points::Vector{Point_2D})
    return Triangle_2D(facepoints(face, points))
end

# Return a Quadrilateral_2D from the point IDs in a face
function materialize_face(face::SVector{4, UInt32}, points::Vector{Point_2D})
    return Quadrilateral_2D(facepoints(face, points))
end

# Return a Triangle6_2D from the point IDs in a face
function materialize_face(face::SVector{6, UInt32}, points::Vector{Point_2D})
    return Triangle6_2D(facepoints(face, points))
end

# Return a Quadrilateral8_2D from the point IDs in a face
function materialize_face(face::SVector{8, UInt32}, points::Vector{Point_2D})
    return Quadrilateral8_2D(facepoints(face, points))
end

# Return an SVector of the points in the edge
function materialize_face(face_id::UInt32, mesh::UnstructuredMesh_2D)
    return materialize_face(mesh.faces[face_id], mesh.points)
end

# Return a materialized face for each face in the mesh
function materialize_faces(mesh::UnstructuredMesh_2D)
    return materialize_face.(mesh.faces, Ref(mesh.points))
end

# Return the number of edges in a face
function num_edges(face::SVector{L, UInt32}) where {L}
    if L % 3 === 0
        return 0x00000003
    elseif L % 4 === 0
        return 0x00000004
    else
        # Error
        return 0x00000000
    end
end

function remap_points_to_hilbert(points::Vector{Point_2D})
    bb = boundingbox(points)
    npoints = length(points)
    # Generate a Hilbert curve
    hilbert_points = hilbert_curve(bb, npoints)
    nhilbert_points = length(hilbert_points)
    # For each point, get the index of the hilbert points that is closest
    point_indices = Vector{Int64}(undef, npoints)
    for i ∈ 1:npoints
        min_distance = 1.0e10
        for j ∈ 1:nhilbert_points
            pdistance = distance²(points[i], hilbert_points[j])
            if pdistance < min_distance
                min_distance = pdistance
                point_indices[i] = j
            end
        end
    end
    return sortperm(point_indices)
end

function reorder_points_to_hilbert!(mesh::UnstructuredMesh_2D)
    # Points
    # point_map     maps  new_points[i] == mesh.points[point_map[i]]
    # point_map_inv maps mesh.points[i] == new_points[point_map_inv[i]]
    point_map  = remap_points_to_hilbert(mesh.points)
    point_map_inv = UInt32.(sortperm(point_map))
    # reordered to resemble a hilbert curve
    permute!(mesh.points, point_map)

    # Adjust face indices
    # Point IDs have changed, so we need to change the point IDs referenced by the faces
    new_faces_vec = [ point_map_inv[face] for face in mesh.faces]
    for i in 1:length(mesh.faces)
        mesh.faces[i] = SVector(point_map_inv[mesh.faces[i]])
    end
    return nothing
end

function reorder_faces_to_hilbert!(mesh::UnstructuredMesh_2D)
    face_map  = UInt32.(remap_points_to_hilbert(centroid.(materialize_faces(mesh))))
    permute!(mesh.faces, face_map)
    return nothing
end

function reorder_to_hilbert!(mesh::UnstructuredMesh_2D)
    reorder_points_to_hilbert!(mesh)
    reorder_faces_to_hilbert!(mesh)
    return nothing
end

# Return the ID of the edge shared by two adjacent faces
function shared_edge(face1::UInt32, face2::UInt32, mesh::UnstructuredMesh_2D)
    for edge1 in mesh.face_edge_connectivity[face1]
        for edge2 in mesh.face_edge_connectivity[face2]
            if edge1 == edge2
                return edge1
            end
        end
    end
    return 0x00000000
end

# How to display a mesh in REPL
function Base.show(io::IO, mesh::UnstructuredMesh_2D)
    mesh_type = typeof(mesh)
    println(io, mesh_type)
    println(io, "  ├─ Name      : $(mesh.name)")
    size_MB = Base.summarysize(mesh)/1E6
    if size_MB < 1
        size_KB = size_MB*1000
        println(io, "  ├─ Size (KB) : $size_KB")
    else
        println(io, "  ├─ Size (MB) : $size_MB")
    end
    println(io, "  ├─ Points    : $(length(mesh.points))")
    nedges = length(mesh.edges)
    println(io, "  ├─ Edges     : $nedges")
    println(io, "  │  ├─ Linear         : $(count(x->x isa SVector{2, UInt32},  mesh.edges))")
    println(io, "  │  ├─ Quadratic      : $(count(x->x isa SVector{3, UInt32},  mesh.edges))")
    println(io, "  │  └─ Materialized?  : $(length(mesh.materialized_edges) !== 0)")
    println(io, "  ├─ Faces     : $(length(mesh.faces))")
    println(io, "  │  ├─ Triangle       : $(count(x->x isa SVector{3, UInt32},  mesh.faces))")
    println(io, "  │  ├─ Quadrilateral  : $(count(x->x isa SVector{4, UInt32},  mesh.faces))")
    println(io, "  │  ├─ Triangle6      : $(count(x->x isa SVector{6, UInt32},  mesh.faces))")
    println(io, "  │  ├─ Quadrilateral8 : $(count(x->x isa SVector{8, UInt32},  mesh.faces))")
    println(io, "  │  └─ Materialized?  : $(length(mesh.materialized_faces) !== 0)")
    println(io, "  ├─ Connectivity")
    println(io, "  │  ├─ Edge/Face : $(0 < length(mesh.edge_face_connectivity))")
    println(io, "  │  └─ Face/Edge : $(0 < length(mesh.face_edge_connectivity))")
    println(io, "  ├─ Boundary edges")
    nsides = length(mesh.boundary_edges)
    if nsides !== 0
        println(io, "  │  ├─ Edges : $(mapreduce(x->length(x), +, mesh.boundary_edges))")
    else
        println(io, "  │  ├─ Edges : 0")
    end
    println(io, "  │  └─ Sides : $nsides")
    println(io, "  └─ Face sets : $(length(keys(mesh.face_sets)))")
end

# Sort intersection points, deleting points that are less than minimum_segment_length apart
function sort_intersection_points!(p::Point_2D, points::Vector{Point_2D})
    if 2 <= length(points)
        sortpoints!(p, points)
        # Eliminate any points for which the distance between consecutive points
        # is less than the minimum segment length
        delete_ids = Int64[]
        id_start = 1
        for id_stop ∈ 2:length(points)
            if distance²(points[id_start], points[id_stop]) < minimum_segment_length^2
                push!(delete_ids, id_stop)
            else
                id_start = id_stop
            end
        end
        deleteat!(points, delete_ids)
    else
        return points
    end
end

# Return a mesh with name name, composed of the faces in the set face_ids
function submesh(name::String, mesh::M) where {M <: UnstructuredMesh_2D}
    # Setup faces and get all vertex ids
    face_ids = mesh.face_sets[name]
    submesh_faces = Vector{Vector{UInt32}}(undef, length(face_ids))
    vertex_ids = Set{UInt32}()
    for (i, face_id) in enumerate(face_ids)
        face_vec = collect(mesh.faces[face_id].data)
        submesh_faces[i] = face_vec
        union!(vertex_ids, Set(face_vec))
    end
    # Need to remap vertex ids in faces to new ids
    vertex_ids_sorted = sort(collect(vertex_ids))
    vertex_map = Dict{UInt32, UInt32}()
    for (i,v) in enumerate(vertex_ids_sorted)
        vertex_map[v] = i
    end
    submesh_points = Vector{Point_2D}(undef, length(vertex_ids_sorted))
    for (i, v) in enumerate(vertex_ids_sorted)
        submesh_points[i] = mesh.points[v]
    end
    # remap vertex ids in faces
    for face in submesh_faces
        for (i, v) in enumerate(face)
            face[i] = vertex_map[v]
        end
    end
    # At this point we have points, faces, & name.
    # Just need to get the face sets
    submesh_face_sets = Dict{String, Set{UInt32}}()
    for face_set_name in keys(mesh.face_sets)
        set_intersection = mesh.face_sets[face_set_name] ∩ face_ids
        if length(set_intersection) !== 0
            submesh_face_sets[face_set_name] = set_intersection
        end
    end
    # Need to remap face ids in face sets
    face_map = Dict{UInt32, UInt32}()
    for (i,f) in enumerate(face_ids)
        face_map[f] = i
    end
    for face_set_name in keys(submesh_face_sets)
        new_set = Set{UInt32}()
        for fid in submesh_face_sets[face_set_name]
            union!(new_set, face_map[fid])
        end
        submesh_face_sets[face_set_name] = new_set
    end
    return M(name = name,
             points = submesh_points,
             faces = [SVector{length(f), UInt32}(f) for f in submesh_faces],
             face_sets = submesh_face_sets
            )
end

# Plot
# -------------------------------------------------------------------------------------------------
if enable_visualization
    function convert_arguments(LS::Type{<:LineSegments}, mesh::UnstructuredMesh_2D)
        if 0 < length(mesh.materialized_edges)
            return convert_arguments(LS, mesh.materialized_edges)
        elseif 0 < length(mesh.edges)
            return convert_arguments(LS, materialize_edges(mesh))
        else
            return convert_arguments(LS, materialize_faces(mesh))
        end
    end

    function convert_arguments(P::Type{<:Mesh}, mesh::UnstructuredMesh_2D)
        if 0 < length(mesh.materialized_faces)
            return convert_arguments(P, mesh.materialized_faces)
        else
            return convert_arguments(P, materialize_faces(mesh))
        end
    end
end

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
function mimestring(mime::Type{<:MIME}, x)
    buf=IOBuffer()
    show(buf, mime(), x)
    String(take!(buf))
end
mimestring(x) = mimestring(MIME"text/plain", x)

struct MultiTerm <: AbstractTerm
    terms::StatsModels.TupleTerm
end
StatsModels.apply_schema(mt::MultiTerm, sch::StatsModels.Schema, Mod::Type) =
    apply_schema.(mt.terms, Ref(sch), Mod)

@testset "terms" begin

    using Statistics

    @testset "concrete_term" begin
        t = term(:aaa)
        @test string(t) == "aaa"
        @test mimestring(t) == "aaa(unknown)"

        t0 = concrete_term(t, [3, 2, 1])
        @test string(t0) == "aaa"
        @test mimestring(t0) == "aaa(continuous)"
        @test t0.mean == 2.0
        @test t0.var == var([1,2,3])
        @test t0.min == 1.0
        @test t0.max == 3.0

        t1 = concrete_term(t, [:a, :b, :c])
        @test t1.contrasts isa StatsModels.ContrastsMatrix{DummyCoding}
        @test string(t1) == "aaa"
        @test mimestring(t1) == "aaa(DummyCoding:3→2)"

        t3 = concrete_term(t, [:a, :b, :c], DummyCoding())
        @test t3.contrasts isa StatsModels.ContrastsMatrix{DummyCoding}
        @test string(t3) == "aaa"
        @test mimestring(t3) == "aaa(DummyCoding:3→2)"

        t2 = concrete_term(t, [:a, :a, :b], EffectsCoding())
        @test t2.contrasts isa StatsModels.ContrastsMatrix{EffectsCoding}
        @test mimestring(t2) == "aaa(EffectsCoding:2→1)"
        @test string(t2) == "aaa"

        t2full = concrete_term(t, [:a, :a, :b], StatsModels.FullDummyCoding())
        @test t2full.contrasts isa StatsModels.ContrastsMatrix{StatsModels.FullDummyCoding}
        @test mimestring(t2full) == "aaa(StatsModels.FullDummyCoding:2→2)"
        @test string(t2full) == "aaa"
    end

    @testset "term operators" begin
        a = term(:a)
        b = term(:b)
        @test a + b == (a, b)
        @test (a ~ b) == FormulaTerm(a, b)
        @test string(a~b) == "$a ~ $b"
        @test mimestring(a~b) ==
            """FormulaTerm
               Response:
                 a(unknown)
               Predictors:
                 b(unknown)"""
        @test mimestring(a ~ term(1) + b) ==
            """FormulaTerm
               Response:
                 a(unknown)
               Predictors:
                 1
                 b(unknown)"""
        @test a & b == InteractionTerm((a,b))
        @test string(a & b) == "$a & $b"
        @test mimestring(a & b) == "a(unknown) & b(unknown)"
        c = term(:c)
        @test (a+b)+c == (a,b,c)
        @test a+(b+c) == (a,b,c)
    end

    @testset "expand nested tuples of terms during apply_schema" begin
        sch = schema((a=rand(10), b=rand(10), c=rand(10)))

        # nested tuples of terms are expanded by apply_schema
        terms = (term(:a), (term(:b), term(:c)))
        terms2 = apply_schema(terms, sch, Nothing)
        @test terms2 isa NTuple{3, ContinuousTerm}
        @test terms2 == apply_schema(term.((:a, :b, :c)), sch, Nothing)

        # a term that generates multiple terms after apply_schema
        mterms = (terms[1], MultiTerm(terms[2]))
        terms3 = apply_schema(mterms, sch, Nothing)

        @test terms2 == terms3
    end

    @testset "Intercept and response traits" begin

        has_responses = [term(:y), term(1), InterceptTerm{true}(), term(:y)+term(:z),
                        term(:y) + term(0), term(:y) + InterceptTerm{false}()]
        no_responses = [term(0), InterceptTerm{false}()]

        has_intercepts = [term(1), InterceptTerm{true}()]
        omits_intercepts = [term(0), term(-1), InterceptTerm{false}()]

        using StatsModels: hasresponse, hasintercept, omitsintercept

        a = term(:a)

        for lhs in has_responses, rhs in has_intercepts
            @test hasresponse(lhs ~ rhs)
            @test hasintercept(lhs ~ rhs)
            @test !omitsintercept(lhs ~ rhs)

            @test hasresponse(lhs ~ rhs + a)
            @test hasintercept(lhs ~ rhs + a)
            @test !omitsintercept(lhs ~ rhs + a)

        end

        for lhs in no_responses, rhs in has_intercepts
            @test !hasresponse(lhs ~ rhs)
            @test hasintercept(lhs ~ rhs)
            @test !omitsintercept(lhs ~ rhs)

            @test !hasresponse(lhs ~ rhs + a)
            @test hasintercept(lhs ~ rhs + a)
            @test !omitsintercept(lhs ~ rhs + a)
        end

        for lhs in has_responses, rhs in omits_intercepts
            @test hasresponse(lhs ~ rhs)
            @test !hasintercept(lhs ~ rhs)
            @test omitsintercept(lhs ~ rhs)

            @test hasresponse(lhs ~ rhs + a)
            @test !hasintercept(lhs ~ rhs + a)
            @test omitsintercept(lhs ~ rhs + a)
        end

        for lhs in no_responses, rhs in omits_intercepts
            @test !hasresponse(lhs ~ rhs)
            @test !hasintercept(lhs ~ rhs)
            @test omitsintercept(lhs ~ rhs)

            @test !hasresponse(lhs ~ rhs + a)
            @test !hasintercept(lhs ~ rhs + a)
            @test omitsintercept(lhs ~ rhs + a)
        end

    end

end


text:
using ARGTools
using TreeKnit
using TreeTools
using Test

n = 15
ρ = 0.1
N = 10_000

f() = f(n, ρ, N)
function f(n, ρ, N)
	arg_ref = ARGTools.SimulateARG.simulate(N, ρ/N*n, n; simtype=:kingman)
	t1, t2 = ARGTools.trees_from_ARG(arg_ref)
	MCCs = ARGTools.MCCs_from_arg(arg_ref)

	arg = TreeKnit.SRG.arg_from_trees(t1, t2, MCCs)[1]
	tt1, tt2 = TreeKnit.SRG.trees_from_arg(arg)
	TreeTools.remove_internal_singletons!(tt1)
	TreeTools.remove_internal_singletons!(tt2)

	S1_ref = SplitList(t1)
	S1 = SplitList(tt1)
	S2_ref = SplitList(t2)
	S2 = SplitList(tt2)

	return S1 == S1_ref && S2 == S2_ref
end

@testset "Random args" begin
	@test mapreduce(i->f(), *, 1:50, init=true)
end

max_stars_repo_name: PierreBarrat/RecombTools

max_issues_repo_name: PierreBarrat/TreeKnit

max_forks_repo_name: PierreBarrat/TreeKnit


text:
module TestRanges

using Test
using MLJBase
using MLJTuning
using Random
import Distributions
const Dist = Distributions

# `in` for MLJType is overloaded to be `===` based. For purposed of
# testing here, we need `==` based:
function _in(x, itr)::Union{Bool,Missing}
    for y in itr
        ismissing(y) && return missing
        y == x && return true
    end
    return false
end
_issubset(itr1, itr2) = all(_in(x, itr2) for x in itr1)

@testset "boundedness traits" begin
    r1 = range(Float64, :K, lower=1, upper=10)
    r2 = range(Float64, :K, lower=-1, upper=Inf, origin=1, unit=1)
    r3 = range(Float64, :K, lower=0, upper=Inf, origin=1, unit=1)
    r4 = range(Float64, :K, lower=-Inf, upper=1, origin=0, unit=1)
    r5 = range(Float64, :K, lower=-Inf, upper=Inf, origin=1, unit=1)
    @test MLJTuning.boundedness(r1) == MLJTuning.Bounded
    @test MLJTuning.boundedness(r2) == MLJTuning.Other
    @test MLJTuning.boundedness(r3) == MLJTuning.PositiveUnbounded
    @test MLJTuning.boundedness(r4) == MLJTuning.Other
    @test MLJTuning.boundedness(r5) == MLJTuning.Other
end

mutable struct DummyModel <: Deterministic
    lambda::Float64
    metric::Float64
    kernel::Char
end

dummy_model = DummyModel(4, 9.5, 'k')

mutable struct SuperModel <: Deterministic
    K::Int64
    model1::DummyModel
    model2::DummyModel
end

dummy_model = DummyModel(1.2, 9.5, 'k')
super_model = SuperModel(4, dummy_model, deepcopy(dummy_model))

r1 = range(super_model, :(model1.kernel), values=['c', 'd'])
r2 = range(super_model, :K, lower=1, upper=10, scale=:log10)

@testset "models from cartesian range and resolutions" begin

    f1 = r1.field
    f2 = r2.field
    itr1 = iterator(r1, nothing)
    itr2 = iterator(r2, 7)

    # with method:
    m1 = MLJTuning.grid(super_model, [f1, f2], [itr1, itr2])
    m1r = MLJTuning.grid(MersenneTwister(123), super_model, [f1, f2],
                         [itr1, itr2])

    # generate all models by hand:
    models1 = [SuperModel(1, DummyModel(1.2, 9.5, 'c'), dummy_model),
               SuperModel(1, DummyModel(1.2, 9.5, 'd'), dummy_model),
               SuperModel(2, DummyModel(1.2, 9.5, 'c'), dummy_model),
               SuperModel(2, DummyModel(1.2, 9.5, 'd'), dummy_model),
               SuperModel(3, DummyModel(1.2, 9.5, 'c'), dummy_model),
               SuperModel(3, DummyModel(1.2, 9.5, 'd'), dummy_model),
               SuperModel(5, DummyModel(1.2, 9.5, 'c'), dummy_model),
               SuperModel(5, DummyModel(1.2, 9.5, 'd'), dummy_model),
               SuperModel(7, DummyModel(1.2, 9.5, 'c'), dummy_model),
               SuperModel(7, DummyModel(1.2, 9.5, 'd'), dummy_model),
               SuperModel(10, DummyModel(1.2, 9.5, 'c'), dummy_model),
               SuperModel(10, DummyModel(1.2, 9.5, 'd'), dummy_model)]

    @test _issubset(models1, m1) && _issubset(m1, models1)
    @test m1r != models1
    @test _issubset(models1, m1r) && _issubset(m1, models1)

    itr1 = iterator(r1, 1)

    # with method:
    m2 = MLJTuning.grid(super_model, [f1, f2], [itr1, itr2])

    # generate all models by hand:
    models2 = [SuperModel(1, DummyModel(1.2, 9.5, 'c'), dummy_model),
               SuperModel(2, DummyModel(1.2, 9.5, 'c'), dummy_model),
               SuperModel(3, DummyModel(1.2, 9.5, 'c'), dummy_model),
               SuperModel(5, DummyModel(1.2, 9.5, 'c'), dummy_model),
               SuperModel(7, DummyModel(1.2, 9.5, 'c'), dummy_model),
               SuperModel(10, DummyModel(1.2, 9.5, 'c'), dummy_model)]

    @test _issubset(models2, m2) && _issubset(m2, models2)

end

@testset "processing user specification of range in Grid" begin
    r1 = range(Int, :h1, lower=1, upper=10)
    r2 = range(Int, :h2, lower=20, upper=30)
    s = range(Char, :j1, values = ['x', 'y'])
    @test_throws ArgumentError MLJTuning.process_grid_range("junk", 42, 1)
    @test(@test_logs((:warn, r"Ignoring"),
                     MLJTuning.process_grid_range((s, 3), 42, 1)) ==
          ((s, ), (2, )))
    @test MLJTuning.process_grid_range(r1, 42, 1) == ((r1, ), (42, ))
    @test MLJTuning.process_grid_range((r1, 3), 42, 1) == ((r1, ), (3, ))
    @test MLJTuning.process_grid_range(s, 42, 1) == ((s, ), (2,))
    @test MLJTuning.process_grid_range([(r1, 3), r2, s], 42, 1) ==
        ((r1, r2, s), (3, 42, 2))
end

struct MySampler end
Base.rand(rng::AbstractRNG, ::MySampler) = rand(rng)

@testset "processing user specification of range in RandomSearch" begin
    r1 = range(Int, :h1, lower=1, upper=10, scale=exp)
    r2 = range(Int, :h2, lower=5, upper=Inf, origin=10, unit=5)
    r3 = range(Char, :j1, values = ['x', 'y'])
    s = MySampler()

    @test_throws(ArgumentError,
                 MLJTuning.process_random_range("junk",
                                                Dist.Uniform,
                                                Dist.Gamma,
                                                Dist.Cauchy))
    @test_throws(ArgumentError,
                 MLJTuning.process_random_range((r1, "junk"),
                                                Dist.Uniform,
                                                Dist.Gamma,
                                                Dist.Cauchy))
    @test_throws(ArgumentError,
                 MLJTuning.process_random_range((r3, "junk"),
                                                Dist.Uniform,
                                                Dist.Gamma,
                                                Dist.Cauchy))
    @test_throws(ArgumentError,
                 MLJTuning.process_random_range(("junk", s),
                                                Dist.Uniform,
                                                Dist.Gamma,
                                                Dist.Cauchy))

    # unpaired numeric range:
    pp = MLJTuning.process_random_range(r1,
                                         Dist.Uniform, # bounded
                                         Dist.Gamma,   # positive_unbounded
                                         Dist.Cauchy)  # other
    @test pp isa Tuple{Tuple{Symbol,MLJBase.NumericSampler}}
    p = first(pp)
    @test first(p) == :h1
    s = last(p)
    @test s.scale == r1.scale
    @test s.distribution == Dist.Uniform(1.0, 10.0)

    # unpaired nominal range:
    p = MLJTuning.process_random_range(r3,
                                       Dist.Uniform,
                                       Dist.Gamma,
                                       Dist.Cauchy)  |> first
    @test first(p) == :j1
    s = last(p)
    @test s.values == r3.values
    @test s.distribution.p == [0.5, 0.5]
    @test s.distribution.support == 1:2

    # (numeric range, distribution instance):
    p = MLJTuning.process_random_range((r2, Dist.Poisson(3)),
                                       Dist.Uniform,
                                       Dist.Gamma,
                                       Dist.Cauchy) |> first
    @test first(p) == :h2
    s = last(p)
    @test s.scale == r2.scale
    @test s.distribution == Dist.truncated(Dist.Poisson(3.0), 5.0, Inf)

    # (numeric range, distribution type):
    p = MLJTuning.process_random_range((r2, Dist.Poisson),
                                       Dist.Uniform,
                                       Dist.Gamma,
                                       Dist.Cauchy) |> first
    s = last(p)
    @test s.distribution == Dist.truncated(Dist.Poisson(r2.unit), 5.0, Inf)

end

end
true


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

using Markdown
using InteractiveUtils

# ╔═╡ 5168c034-1879-11eb-2df4-e74ba7880da1
using Pkg, DrWatson

# ╔═╡ 5168ed8e-1879-11eb-31cd-d3a727fff08d
begin
	@quickactivate "SR2TuringPluto"
	using Turing
	using StatisticalRethinking
	Turing.setprogress!(false);
end;

# ╔═╡ e6220984-1878-11eb-275c-db7fe829037a
md"## Intro clip 1_priors.jl"

# ╔═╡ 51695ecc-1879-11eb-2fd1-b1d0abfb4b8f
begin
	df = CSV.read(sr_datadir("Howell1.csv"), DataFrame)
	df = df[df.age .>= 18, :]
	x̄ = mean(df.weight)
	x = range(minimum(df.weight), maximum(df.weight), length = 100)
	Text(precis(df; io=String))
end

# ╔═╡ 5174d414-1879-11eb-28ed-b3f36412a25f
@model function ppl4_3(weights, heights)
    a ~ Normal(mean(weights), 50)
    b ~ Normal(0, 5)
    σ ~ Uniform(0, 20)
    μ = a .+ b .* weights
    for i in eachindex(heights)
        heights[i] ~ Normal(μ[i], σ)
    end
end

# ╔═╡ 11151d8c-1916-11eb-2d70-69180f7be854
begin
	m4_3t = ppl4_3(df.weight, df.height)
	priors4_3t = sample(m4_3t, Prior(), 50)
	Text(sprint(show, "text/plain", priors4_3t))
end

# ╔═╡ 1f9c6586-1916-11eb-055d-9dd844c9b0a7
begin
	priors4_3t_df = DataFrame(priors4_3t)
	Text(precis(priors4_3t_df; io=String))
end

# ╔═╡ f1943d70-1916-11eb-2cc5-5b6afaea5fbd
begin
	plot(xlims=(25, 70), ylims=(-10, 300), title="Possible prior regression", leg=false)
	for row in eachrow(priors4_3t_df)
		plot!(df.weight, row.a .+ row.b * df.weight)
	end
	scatter!(df.weight, df.height)
end

# ╔═╡ 0365d298-187f-11eb-1dc1-b7f61a13c7b3
md"## End of intro clip 1_priors.jl"

# ╔═╡ Cell order:
# ╟─e6220984-1878-11eb-275c-db7fe829037a
# ╠═5168c034-1879-11eb-2df4-e74ba7880da1
# ╠═5168ed8e-1879-11eb-31cd-d3a727fff08d
# ╠═51695ecc-1879-11eb-2fd1-b1d0abfb4b8f
# ╠═5174d414-1879-11eb-28ed-b3f36412a25f
# ╠═11151d8c-1916-11eb-2d70-69180f7be854
# ╠═1f9c6586-1916-11eb-055d-9dd844c9b0a7
# ╠═f1943d70-1916-11eb-2cc5-5b6afaea5fbd
# ╟─0365d298-187f-11eb-1dc1-b7f61a13c7b3