text:
function copy_variable!(toModel, variable)
    map(variable) do x
        JuMP.add_variable(toModel, JuMP.build_variable(error, (get_info(x))))
    end
end

function copy_variable!(toModel, variable, f)
    map(variable) do x
        JuMP.add_variable(toModel, JuMP.build_variable(error, f(get_info(x))))
    end
end

function copy_variable!(toModel, variable::JuMP.VariableRef)
    JuMP.add_variable(toModel, JuMP.build_variable(error, get_info(variable)))
end

function copy_variable!(toModel, variable::JuMP.VariableRef, f)
    JuMP.add_variable(toModel, JuMP.build_variable(error, f(get_info(variable))), "exp")
end

# constuct variable info object for a single variable
function get_info(x::VariableRef)
    has_lb_local = false
    lb_local = NaN
    has_ub_local = false
    ub_local = NaN
    is_fixed_local = false
    fixed_value_local = NaN
    has_start_local = false
    start_value_local = NaN
    is_binary_local = false
    is_integer_local = false

    if has_lower_bound(x)
        has_lb_local = true
        lb_local = lower_bound(x)
    end

    if has_upper_bound(x)
        has_ub_local = true
        ub_local = upper_bound(x)
    end

    if is_fixed(x)
        is_fixed_local = true
        fixed_value_local = fix_value(x)
    end

    if start_value(x) != nothing
        has_start_local = true
        start_value_local = start_value(x)
    end

    if is_binary(x)
        is_binary_local = true
    end

    if is_integer(x)
        is_integer_local = true
    end

    VariableInfo(
        has_lb_local,
        lb_local,
        has_ub_local,
        ub_local,
        is_fixed_local,
        fixed_value_local,
        has_start_local,
        start_value_local,
        is_binary_local,
        is_integer_local,
    )
end

function relaxbinary(x::VariableInfo)
    VariableInfo(
        true,
        0.0,
        true,
        1.0,
        x.has_fix,
        x.fixed_value,
        x.has_start,
        x.start,
        false,
        false,
    )
end

function relaxinteger(x::VariableInfo)
    VariableInfo(
        true,
        x.lower_bound,
        true,
        x.upper_bound,
        x.has_fix,
        x.fixed_value,
        x.has_start,
        x.start,
        false,
        false,
    )
end

function UnitIntervalInformation(;UB::Float64=1.0)
    VariableInfo(true, 0.0, true, UB, false, NaN, false, NaN, false, false)
end

function objcoef(x::JuMP.VariableRef)
    affine_expression = objective_function(owner_model(x))
    if x in keys(affine_expression.terms)
        affine_expression.terms[x]
    else
        0.0
    end
end

# function coef(aff, x::JuMP.VariableRef)
#     if x in keys(aff.terms)
#         aff.terms[x]
#     else
#         0.0
#     end
# end

function coef(aff, x::JuMP.VariableRef)
    if x in keys(aff)
        aff[x]
    else
        0.0
    end
end

function densekey_to_tuple(key::Any)
    if typeof(key) <: JuMP.Containers.DenseAxisArrayKey
        return length(key.I)==1 ? key.I[1] : key.I
    else
        return key
    end
end

function unpack_expansions(a::Dict{AbstractTree,Dict{Symbol,Any}})
   assign=Array{Expr,1}()
   found=Dict{Symbol,Bool}()
   for (i,j) in a
       for (k,l) in j
           if k in keys(found)
               push!(assign,:($k[$i] = $l))
           else
               if isdefined(Main,k)
                   m=Symbol("→"*string(k))
                   push!(assign,:($m=$k))
               end
               push!(assign,:($k = Dict{AbstractTree,Any}()))
               push!(assign,:($k[$i] = $l))
               found[k]=true
           end
       end
   end
   assign
end

function clear_expansions(a::Dict{AbstractTree,Dict{Symbol,Any}})
   assign=Array{Expr,1}()
   for (i,j) in a
       for (k,l) in j
           m=Symbol("→"*string(k))
           if isdefined(Main,m)
               push!(assign,:($k=$m))
               push!(assign,:($m=nothing))
           else
               push!(assign,:($k=nothing))
           end
       end
       break
   end
   assign
end

function compute_objval(scenarios::Dict{Leaf,Float64}, probabilities::Dict{AbstractTree,Float64}, risk::Union{Risk,Array{Risk,1}})
    scenario_objs=Array{Tuple{Float64,Float64,Leaf},1}()
	EV_weight=1.0
	EV=0.0
	for (leaf,val) in scenarios
		push!(scenario_objs,(probabilities[leaf],val,leaf))
		EV+=probabilities[leaf]*val
	end

	obj=0.0

    if typeof(risk)==Risk
		if risk.α==1.0 || risk.λ==0.0
			risk=[]
		else
			risk=[risk]
		end
	end
    for i in 1:length(risk)
        EV_weight-=risk[i].λ
        so=scenario_objs
        if risk[i].offset!=nothing
            for j in 1:length(so)
                so[j]=(so[j][1],so[j][2]-risk[i].offset[so[j][3]],so[j][3])
            end
        end
        sort!(so,by=i->i[2],rev=true)
        beta=risk[i].α
        for scen in so
            if scen[1]>beta
                obj+=scen[2]*risk[i].λ*beta/risk[i].α
                beta=0
            else
                obj+=scen[2]*risk[i].λ*scen[1]/risk[i].α
                beta-=scen[1]
            end
        end
    end

	obj+EV_weight*EV
end

function solution_to_dictionary(jmodel::JuDGEModel;prefix="")
    function helper(jmodel::JuDGEModel,node::AbstractTree,solution::T where {T <: Dict})
        vars=all_variables(jmodel.sub_problems[node])
        for v in vars
            temp=string(v)
            i=findfirst('[',temp)
            if i==nothing
                if Symbol(prefix*temp) ∉ keys(solution)
                    solution[Symbol(prefix*temp)]=Dict{AbstractTree,Float64}()
                end
                solution[Symbol(prefix*temp)][node]=JuMP.value(v)
            else
                if Symbol(temp[1:i-1]) ∉ keys(solution)
                    solution[Symbol(prefix*temp[1:i-1])]=Dict{String,Dict{AbstractTree,Float64}}()
                end
                if temp[i+1:length(temp)-1] ∉ keys(solution[Symbol(prefix*temp[1:i-1])])
                    solution[Symbol(prefix*temp[1:i-1])][temp[i+1:length(temp)-1]]=Dict{AbstractTree,Float64}()
                end
                solution[Symbol(prefix*temp[1:i-1])][temp[i+1:length(temp)-1]][node]=JuMP.value(v)
            end
        end

        for (x,var) in jmodel.master_problem.ext[:expansions][node]
             if typeof(var) <: AbstractArray
                 if Symbol(prefix*string(x)*"_master") ∉ keys(solution)
                     solution[Symbol(prefix*string(x)*"_master")]=Dict{String,Dict{AbstractTree,Float64}}()
                 end
                 val=JuMP.value.(var)
                 if typeof(var) <: JuMP.Containers.SparseAxisArray
                     val=val.data
                 end
                 for key in keys(val)
                     temp=""
                     if typeof(val) <: Array
                         strkey=string(key)
                         strkey=replace(strkey,"CartesianIndex("=>"")
                         strkey=replace(strkey,")"=>"")
                         strkey=replace(strkey,", "=>",")
                         temp*=strkey
                     elseif typeof(val) <: Dict
                         strkey=string(key)
                         strkey=replace(strkey,")"=>"")
                         strkey=replace(strkey,"("=>"")
                         strkey=replace(strkey,", "=>",")
                         temp*=strkey
                     else
                         for i in 1:length(val.axes)-1
                            temp*=string(key[i])*","
                         end
                         temp*=string(key[length(val.axes)])
                     end
                     if temp ∉ keys(solution[Symbol(prefix*string(x)*"_master")])
                         solution[Symbol(prefix*string(x)*"_master")][temp]=Dict{AbstractTree,Float64}()
                     end
                     solution[Symbol(prefix*string(x)*"_master")][temp][node]=val[key]
                 end
             else
                 if Symbol(prefix*string(x)*"_master") ∉ keys(solution)
                     solution[Symbol(prefix*string(x)*"_master")]=Dict{AbstractTree,Float64}()
                 end
                 solution[Symbol(prefix*string(x)*"_master")][node]=JuMP.value(var)
             end
        end

        if typeof(node)==Tree
            for child in node.children
                helper(jmodel,child,solution)
            end
        else
            if Symbol(prefix*"scenario_obj") ∉ keys(solution)
                solution[Symbol(prefix*"scenario_obj")]=Dict{AbstractTree,Float64}()
            end
            solution[Symbol(prefix*"scenario_obj")][node]=JuMP.value(jmodel.master_problem.ext[:scenprofit_var][node])
        end
        solution
    end

    if termination_status(jmodel.master_problem) != MathOptInterface.OPTIMAL && termination_status(jmodel.master_problem) != MathOptInterface.INTERRUPTED && termination_status(jmodel.master_problem) != MathOptInterface.TIME_LIMIT  && termination_status(jmodel.master_problem) != MathOptInterface.LOCALLY_SOLVED
        error("You need to first solve the decomposed model.")
    end

    if prefix!=""
        prefix*="_"
    end

    solution=Dict{Symbol,Any}()
    helper(jmodel, jmodel.tree, solution)
end

function set_starting_solution!(deteq::DetEqModel,jmodel::JuDGEModel)
	for node in collect(jmodel.tree)
	   for (name,exps) in jmodel.master_problem.ext[:expansions][node]
	      for index in keys(jmodel.master_problem.ext[:expansions][node][name])
	         key=JuDGE.densekey_to_tuple(index)
	         set_start_value(deteq.problem.ext[:master_vars][node][name][key],JuMP.value(jmodel.master_problem.ext[:expansions][node][name][index]))
	      end
	   end
	end

	for (node,sp) in jmodel.sub_problems
		temp=Dict{String,VariableRef}()
		for variable in all_variables(sp)
			temp[string(variable)]=variable
		end
		for variable in keys(deteq.problem.ext[:vars][node])
			set_start_value(deteq.problem.ext[:vars][node][variable], JuMP.value(temp[string(variable)]))
		end
	end
end

function get_active_columns(jmodel::JuDGEModel;inttol=10^-7)
    active=Dict{AbstractTree,Array{Any,1}}()
    for node in collect(jmodel.tree)
        active[node]=[]
        for col in jmodel.master_problem.ext[:columns][node]
            if JuMP.value(col.var)>inttol
                push!(active[node],(col,JuMP.value(col.var)))
            end
        end
    end

    active
end

function overprint(str)
    print("\e[2K")
    print("\e[1G")
    print(str)
end

function printleft(str)
    print("\u1b[40G")
    print("\u1b[1K")
    print("\u1b[1G")
    print(str)
end

function printright(str)
    print("\u1b[41G")
    print("\u1b[0K")
    print(str)
end


text:
using JSON, JLD2

CMEParams = Dict()
for i in 1:26
    @load pwd()*"/src/CMEParamsData/CMEParams"*string(i)*".jld2" tempDict
    for k in keys(tempDict)
        if k < 50
            CMEParams[k] = tempDict[k]
        end
    end
end

open("src/CMEParamsData/CMEParams.json","w") do f
    JSON.print(f, CMEParams)
end

@save "src/CMEParamsData/CMEParams.jld2" CMEParams

ext: jl

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
#############################################################################
# geom_mean_hypocone.jl
# Constrains T to
#   A #_t B ⪰ T
# where:
#   * A #_t B is the t-weighted geometric mean of A and B:
#     A^{1/2} (A^{-1/2} B A^{-1/2})^t A^{1/2}
#     Parameter t should be in [0,1].
#   * Constraints A ⪰ 0, B ⪰ 0 are added.
#
# Note on parameter fullhyp:
#   In many applications one doesn't need the full hypograph
#     hyp_t = {(A,B,T) : A #_t B ⪰ T}
#   but rather it is enough to work with a convex set C_t that satisfies
#     (A,B,A #_t B) \in C_t
#     (A,B,T) \in C_t  =>  A #_t B ⪰ T
#   In this case one should set fullhyp = false. The SDP description will be
#   (slightly) smaller. (By default fullhyp is set to true).
#
# All expressions and atoms are subtypes of AbstractExpr.
# Please read expressions.jl first.
#
#REFERENCE
#   Ported from CVXQUAD which is based on the paper: "Lieb's concavity
#   theorem, matrix geometric means and semidefinite optimization" by Hamza
#   Fawzi and James Saunderson (arXiv:1512.03401)
#############################################################################

struct GeomMeanHypoCone
    A::AbstractExpr
    B::AbstractExpr
    t::Rational
    size::Tuple{Int,Int}
    fullhyp::Bool

    function GeomMeanHypoCone(
        A::AbstractExpr,
        B::AbstractExpr,
        t::Rational,
        fullhyp::Bool = true,
    )
        if size(A) != size(B)
            throw(DimensionMismatch("A and B must be the same size"))
        end
        n = size(A)[1]
        if size(A) != (n, n)
            throw(DimensionMismatch("A and B must be square"))
        end
        if t < 0 || t > 1
            throw(DomainError(t, "t must be in the range [0, 1]"))
        end
        return new(A, B, t, (n, n), fullhyp)
    end

    function GeomMeanHypoCone(
        A::Value,
        B::AbstractExpr,
        t::Rational,
        fullhyp::Bool = true,
    )
        return GeomMeanHypoCone(Constant(A), B, t, fullhyp)
    end
    function GeomMeanHypoCone(
        A::AbstractExpr,
        B::Value,
        t::Rational,
        fullhyp::Bool = true,
    )
        return GeomMeanHypoCone(A, Constant(B), t, fullhyp)
    end
    function GeomMeanHypoCone(
        A::Value,
        B::Value,
        t::Rational,
        fullhyp::Bool = true,
    )
        return GeomMeanHypoCone(Constant(A), Constant(B), t, fullhyp)
    end

    function GeomMeanHypoCone(
        A::AbstractExprOrValue,
        B::AbstractExprOrValue,
        t::Integer,
        fullhyp::Bool = true,
    )
        return GeomMeanHypoCone(A, B, t // 1, fullhyp)
    end
end

struct GeomMeanHypoConeConstraint <: Constraint
    head::Symbol
    id_hash::UInt64
    T::AbstractExpr
    cone::GeomMeanHypoCone

    function GeomMeanHypoConeConstraint(T::AbstractExpr, cone::GeomMeanHypoCone)
        if size(T) != cone.size
            throw(DimensionMismatch("T must be size $(cone.size)"))
        end
        id_hash = hash((cone.A, cone.B, cone.t, :GeomMeanHypoCone))
        return new(:GeomMeanHypoCone, id_hash, T, cone)
    end

    function GeomMeanHypoConeConstraint(T::Value, cone::GeomMeanHypoCone)
        return GeomMeanHypoConeConstraint(Constant(T), cone)
    end
end

in(T, cone::GeomMeanHypoCone) = GeomMeanHypoConeConstraint(T, cone)

function AbstractTrees.children(constraint::GeomMeanHypoConeConstraint)
    return (
        constraint.T,
        constraint.cone.A,
        constraint.cone.B,
        "t=$(constraint.cone.t)",
    )
end

# For t ∈ [0,1] the t-weighted matrix geometric mean is matrix concave (arxiv:1512.03401).
# So if A and B are convex sets, then T ⪯ A #_t B will be a convex set.
function vexity(constraint::GeomMeanHypoConeConstraint)
    A = vexity(constraint.cone.A)
    B = vexity(constraint.cone.B)
    T = vexity(constraint.T)

    # NOTE: can't say A == NotDcp() because the NotDcp constructor prints a warning message.
    if typeof(A) == ConvexVexity || typeof(A) == NotDcp
        return NotDcp()
    end
    if typeof(B) == ConvexVexity || typeof(B) == NotDcp
        return NotDcp()
    end
    # Copied from vexity(c::GtConstraint)
    vex = -ConcaveVexity() + T
    if vex == ConcaveVexity()
        vex = NotDcp()
    end
    return vex
end

function conic_form!(
    constraint::GeomMeanHypoConeConstraint,
    unique_conic_forms::UniqueConicForms,
)
    if !has_conic_form(unique_conic_forms, constraint)
        A = constraint.cone.A
        B = constraint.cone.B
        t = constraint.cone.t
        T = constraint.T
        fullhyp = constraint.cone.fullhyp

        is_complex =
            sign(A) == ComplexSign() ||
            sign(B) == ComplexSign() ||
            sign(T) == ComplexSign()
        if is_complex
            make_temporary = () -> HermitianSemidefinite(size(A)[1])
        else
            make_temporary = () -> Semidefinite(size(A)[1])
        end

        if fullhyp && t != 0 && t != 1
            W = make_temporary()
            conic_form!(
                W in GeomMeanHypoCone(A, B, t, false),
                unique_conic_forms,
            )
            conic_form!(W ⪰ T, unique_conic_forms)
            cache_conic_form!(
                unique_conic_forms,
                constraint,
                Array{Convex.ConicConstr,1}(),
            )
        else
            p = t.num
            q = t.den

            if t == 0
                #println("geom_mean_hypocone p=$p q=$q t==0")
                conic_form!(A ⪰ 0, unique_conic_forms)
                conic_form!(B ⪰ 0, unique_conic_forms)
                conic_form!(A ⪰ T, unique_conic_forms)
            elseif t == 1
                #println("geom_mean_hypocone p=$p q=$q t==1")
                conic_form!(A ⪰ 0, unique_conic_forms)
                conic_form!(B ⪰ 0, unique_conic_forms)
                conic_form!(B ⪰ T, unique_conic_forms)
            elseif t == 1 // 2
                #println("geom_mean_hypocone p=$p q=$q t==1/2")
                conic_form!([A T; T' B] ⪰ 0, unique_conic_forms)
            elseif ispow2(q)
                #println("geom_mean_hypocone p=$p q=$q ispow2(q)")
                Z = make_temporary()
                if t < 1 / 2
                    conic_form!(
                        Z in GeomMeanHypoCone(A, B, 2 * t, false),
                        unique_conic_forms,
                    )
                    conic_form!([A T; T' Z] ⪰ 0, unique_conic_forms)
                else
                    conic_form!(
                        Z in GeomMeanHypoCone(A, B, 2 * t - 1, false),
                        unique_conic_forms,
                    )
                    conic_form!([B T; T' Z] ⪰ 0, unique_conic_forms)
                end
            elseif ispow2(p) && t > 1 // 2
                #println("geom_mean_hypocone p=$p q=$q ispow2(p) && t>1/2")
                Z = make_temporary()
                conic_form!(
                    Z in GeomMeanHypoCone(A, T, (2 * p - q) // p, false),
                    unique_conic_forms,
                )
                conic_form!([Z T; T B] ⪰ 0, unique_conic_forms)
            elseif t < 1 / 2
                #println("geom_mean_hypocone p=$p q=$q t<1/2")
                X = make_temporary()
                # Decompose t = (p/2^l) * (2^l/q) where l=floor(log2(q))
                l = floor(Int, log2(q))
                conic_form!(
                    X in GeomMeanHypoCone(A, B, p // (2^l), false),
                    unique_conic_forms,
                )
                conic_form!(
                    T in GeomMeanHypoCone(A, X, (2^l) // q, false),
                    unique_conic_forms,
                )
            else
                #println("geom_mean_hypocone p=$p q=$q else")
                conic_form!(
                    T in GeomMeanHypoCone(B, A, 1 - t, false),
                    unique_conic_forms,
                )
            end

            cache_conic_form!(
                unique_conic_forms,
                constraint,
                Array{Convex.ConicConstr,1}(),
            )
            #cache_conic_form!(unique_conic_forms, constraint, ConicConstr([A,B,T], :geom_mean_hypocone, [size(A), size(B), size(T)]))
            #cache_conic_form!(unique_conic_forms, conic_form!(T, unique_conic_forms))
        end
    end
    return get_conic_form(unique_conic_forms, constraint)
end


text:
# Energy calculation

export total_energy,
    kinetic_energy,
    temperature,
    potential_energy

"""
    total_energy(s, neighbors=nothing)

Compute the total energy of the system.
If the interactions use neighbor lists, the neighbors should be computed
first and passed to the function.
"""
total_energy(s, neighbors=nothing) = kinetic_energy(s) + potential_energy(s, neighbors)

kinetic_energy_noconvert(s) = sum(mass.(s.atoms) .* sum.(abs2, s.velocities)) / 2

"""
    kinetic_energy(s)

Compute the kinetic energy of the system.
"""
function kinetic_energy(s::System{D, G, T}) where {D, G, T}
    ke = kinetic_energy_noconvert(s)
    # Convert energy to per mol if required
    if dimension(s.energy_units) == u"𝐋^2 * 𝐌 * 𝐍^-1 * 𝐓^-2"
        return T(uconvert(s.energy_units, ke * Unitful.Na))
    else
        return T(uconvert(s.energy_units, ke))
    end
end

const temp_conversion_factor = ustrip(u"nm^2 * u * K^-1 * ps^-2", Unitful.k)

"""
    temperature(system)

Calculate the temperature of a system from the kinetic energy of the atoms.
"""
function temperature(s::System{D, G, T}) where {D, G, T}
    ke = kinetic_energy_noconvert(s)
    df = 3 * length(s) - 3
    if unit(ke) == NoUnits
        k = T(temp_conversion_factor)
    else
        k = T(uconvert(u"K^-1" * unit(ke), Unitful.k))
    end
    return 2 * ke / (df * k)
end

"""
    potential_energy(s, neighbors=nothing)

Compute the potential energy of the system.
If the interactions use neighbor lists, the neighbors should be computed
first and passed to the function.

    potential_energy(inter::PairwiseInteraction, vec_ij, coord_i, coord_j,
                     atom_i, atom_j, box_size)
    potential_energy(inter::SpecificInteraction, coords_i, coords_j,
                     box_size)
    potential_energy(inter::SpecificInteraction, coords_i, coords_j,
                     coords_k, box_size)
    potential_energy(inter::SpecificInteraction, coords_i, coords_j,
                     coords_k, coords_l, box_size)
    potential_energy(inter, system, neighbors=nothing)

Calculate the potential energy due to a given interation type.
Custom interaction types should implement this function.
"""
function potential_energy(s::System{D, G, T}, neighbors=nothing) where {D, G, T}
    n_atoms = length(s)
    potential = zero(T) * s.energy_units

    for inter in values(s.pairwise_inters)
        if inter.nl_only
            if isnothing(neighbors)
                error("An interaction uses the neighbor list but neighbors is nothing")
            end
            @inbounds for ni in 1:neighbors.n
                i, j, weight_14 = neighbors.list[ni]
                dr = vector(s.coords[i], s.coords[j], s.box_size)
                if weight_14
                    potential += potential_energy(inter, dr, s.coords[i], s.coords[j], s.atoms[i],
                                                    s.atoms[j], s.box_size, true)
                else
                    potential += potential_energy(inter, dr, s.coords[i], s.coords[j], s.atoms[i],
                                                    s.atoms[j], s.box_size)
                end
            end
        else
            for i in 1:n_atoms
                for j in (i + 1):n_atoms
                    dr = vector(s.coords[i], s.coords[j], s.box_size)
                    potential += potential_energy(inter, dr, s.coords[i], s.coords[j], s.atoms[i],
                                                    s.atoms[j], s.box_size)
                end
            end
        end
    end

    for inter_list in values(s.specific_inter_lists)
        potential += potential_energy(inter_list, s.coords, s.box_size)
    end

    for inter in values(s.general_inters)
        potential += potential_energy(inter, s, neighbors)
    end

    return uconvert(s.energy_units, potential)
end

@views function potential_energy(inter_list::InteractionList2Atoms, coords, box_size)
    return sum(potential_energy.(inter_list.inters, coords[inter_list.is], coords[inter_list.js],
                                    (box_size,)))
end

@views function potential_energy(inter_list::InteractionList3Atoms, coords, box_size)
    return sum(potential_energy.(inter_list.inters, coords[inter_list.is], coords[inter_list.js],
                                    coords[inter_list.ks], (box_size,)))
end

@views function potential_energy(inter_list::InteractionList4Atoms, coords, box_size)
    return sum(potential_energy.(inter_list.inters, coords[inter_list.is], coords[inter_list.js],
                                    coords[inter_list.ks], coords[inter_list.ls], (box_size,)))
end


text:
"""
Plot the normalized pump curves for the pumps in a network. Displays the plot on the screen
by default. Use keywords `screen` and `savepath` to control whether to display and/or save
to file.
"""
function plot_pumps(pumps::Dict{String,Any}; normalized=true, screen=true, reuse=true,
                    savepath=nothing)
    m = 50
    Qhat = Array(range(0.0, 2.0, length=m))
    etahat = -Qhat.^2 .+ 2*Qhat
    Ghat = -1/3*Qhat.^2 .+ 4/3
    Phat = 1/3*Qhat .+ 2/3

    p1 = Plots.plot()
    p2 = Plots.plot()
    p3 = Plots.plot()
    plotargs = Dict(:lt=>:scatter, :shape=>:circle, :ms=>6, :msc=>:auto)
    for (i,(key,pump)) in enumerate(pumps)
        Qbep = pump["q_bep"]
        Gbep = pump["g_bep"]
        Pbep = pump["p_bep"]
        etabep = _WM._DENSITY * _WM._GRAVITY * inv(Pbep) * Gbep * Qbep
        head = array_from_tuples(pump["head_curve"])
        eff_is_curve = haskey(pump, "efficiency_curve")
        if eff_is_curve
            eff = array_from_tuples(pump["efficiency_curve"])
            qmin = max(head[1,1], eff[1,1])/Qbep
            qmax = min(head[end,1], eff[end,1])/Qbep
        else
            eff = pump["efficiency"]
            qmin = head[1,1]/Qbep
            qmax = head[end,1]/Qbep
        end

        # create power curves by linear interpolation of head and eff curves
        if size(head)[1] > 1
            headitp = _ITP.LinearInterpolation(head[:,1]/Qbep, head[:,2]/Gbep)
        else
            headitp(x::Float64) = head[1,2]
            headitp(x::Array{Float64}) = head[1,2]*ones(size(x))
        end
        if eff_is_curve
            effitp = _ITP.LinearInterpolation(eff[:,1]/Qbep, eff[:,2]/etabep)
        else
            effitp(x::Float64) = eff
            effitp(x::Array{Float64}) = eff*ones(size(x))
        end
        if qmin==qmax
            q = [1]
            pow = [1]
        else
            q = Array(range(qmin, stop=qmax, length=7)) # can change number of points here;
                                                        # make this a user option?
            pow = headitp(q).*q./effitp(q)
        end

        # plot head curve
        if normalized
            Plots.plot!(p1, head[:,1]/Qbep, head[:,2]/Gbep, mc=i; plotargs...)
        else
            Plots.plot!(p1, head[:,1], head[:,2], mc=i; plotargs...)
            Plots.plot!(p1, Qbep*Qhat, Gbep*Ghat, lc=i, lw=2, label="") # fit curve
        end
        #plot efficiency curve
        if eff_is_curve
            if normalized
                Plots.plot!(p2, eff[:,1]/Qbep, eff[:,2]/etabep, label=pump["name"], mc=i;
                            plotargs...)
            else
                Plots.plot!(p2, eff[:,1], eff[:,2], label=pump["name"], mc=i; plotargs...)
            end
        else # plot the single point
            if normalized
                Plots.plot!(p2, (1, 1), label=pump["name"], mc=i; plotargs...)
            else
                Plots.plot!(p2, (Qbep, etabep), label=pump["name"], mc=i; plotargs...)
            end
        end
        if !normalized # fit curve
            Plots.plot!(p2, Qbep*Qhat, etabep*etahat, lc=i, lw=2, label="")
        end
        # plot power curve
        if normalized
            Plots.plot!(p3, q, pow, ylims=(0,2), mc=i; plotargs...)
        else
            Plots.plot!(p3, q*Qbep, pow*Pbep, mc=i; plotargs...)
            Plots.plot!(p3, Qbep*Qhat, Pbep*Phat, lc=i, lw=2, label="") # fit curve
        end
    end
    if normalized
        # plot normalized fit curves and add axes labels and legend
        Plots.plot!(p1, Qhat, Ghat, lc=:black, lw=2, label="", ylabel=L"\hat{G}",
                    legend=:none)
        Plots.plot!(p2, Qhat, etahat, lc=:black, lw=2, label="", ylabel=L"\hat{\eta}",
                    legend=:bottom)
        Plots.plot!(p3, Qhat, Phat, lc=:black, lw=2, label="", xlabel=L"\hat{Q}",
                    ylabel=L"\hat{P}", legend=:none)
    else
        # add axes labels and legend (not-normalized fit curves already plotted)
        Plots.plot!(p1, ylabel=L"G", legend=:none)
        Plots.plot!(p2, ylabel=L"\eta", legend=:bottom)
        Plots.plot!(p3, xlabel=L"Q", ylabel=L"P", legend=:none)
    end

    width = 400
    height = width*4/6*3
    pfig = Plots.plot(p1, p2, p3, size=(width,height), layout=(3,1), reuse=reuse);
    # display on the screen and/or save
    screen ? display(pfig) : nothing
    savepath != nothing ? Plots.savefig(pfig, savepath) : nothing
end


text:
"""
GSL fileformat

Describes how elementary sectors are built from airblocks (GAR)

See EUROCONTROL NEST Manual for Gsl fileformat description

Output:
Sectors dictionary with the following elements:
* ID        Sector ID           String
* Name      Sector name         String
* Category  Airspace category   String  Single character
* Type      Sectory typeassert  String  FIR=Flight Information Region,
                                        ERSA=Elementary Restricted Airspace,
                                        ES=Elementary Sector
                                        ERAS=Elementary Regulated Airspace
                                        UNK=Unknown type
* Airblocks                      Vector{Airblock}

Airblocks struct with the following elements:
* Name      Airblock name
* LowerFL   Lower flight level  Int64
* UpperFL   Upper flight level  Int64

GSL fileformat example
4
S;BALTIC;;1;_;ES
A;BALTIC;+;0;999
S;BGGLFIR;SONDRESTROM FIR;1;_;FIR
A;002BG;+;0;999
S;BIFAROER;FAEROER CTA;1;_;FIR
A;011BI;+;0;200
S;BIRD;_;9;_;FIR
A;001BI;+;0;999
A;002BI;+;0;999
A;003BI;+;0;999
A;004BI;+;0;999
A;011BI;+;200;999
A;100BI;+;0;999
A;104BI;+;0;999
A;107BI;+;0;999
A;108BI;+;0;999

"""
module Gsl

export read

struct Airblock
    Name::String
    LowerFL::Int64
    UpperFL::Int64
end

struct Sector
    ID::String
    Name::String
    Category::String
    Type::String
    Airblocks::Vector{Airblock}
end

function read(file)
    sectors = Dict{String, Sector}()
    id = ""
    sname = ""
    category = 0
    type = ""
    airblocks = Vector{Airblock}()
    for line in eachline(file)
        line_elements = split(line, ';')
        if line_elements[1] == "S"
            if id != ""
                sectors[id] = Sector(id, sname, category, type, airblocks)
            end
            id = line_elements[2]
            sname = line_elements[3]
            category = line_elements[5]
            type = line_elements[6]
            airblocks = Vector{Airblock}()
        elseif line_elements[1] == "A"
            aname = line_elements[2]
            lowerfl = parse(Int64, line_elements[4])
            upperfl = parse(Int64, line_elements[5])
            airblock = Airblock(aname, lowerfl, upperfl)
            airblocks = vcat(airblocks, airblock)
        end
    end
    sectors[id] = Sector(id, sname, category, type, airblocks)
    return sectors
end

end # module


text:
# create pde structure and mesh structure
pde[1], mesh[1] = Preprocessing.initializeexasim();

# Define PDE model: governing equations, initial solutions, and boundary conditions
pde[1].model = "ModelD";            # ModelC, ModelD, ModelW
pde[1].modelfile = "pdemodel1";    # name of a file defining the PDE model
include(pde[1].modelfile * ".jl"); # include the PDE model file
pde[1].modelnumber = 1;

# Set discretization parameters, physical parameters, and solver parameters
pde[1].porder = 3;                  # polynomial degree
pde[1].physicsparam = [1.0 0.0];    # thermal conductivity and boundary value
pde[1].tau = [1.0];                 # DG stabilization parameter

# Choose computing platform and set number of processors
#pde[1].platform = "gpu";           # choose this option if NVIDIA GPUs are available
pde[1].mpiprocs = 1;                # number of MPI processors

# create a linear mesh for a square domain
mesh[1].p, mesh[1].t = Mesh.SquareMesh(8,8,1); # a mesh of 8 by 8 quadrilaterals
# expressions for disjoint boundaries
mesh[1].boundaryexpr = [p -> (p[2,:] .< 1e-3), p -> (p[1,:] .> 1-1e-3), p -> (p[2,:] .> 1-1e-3), p -> (p[1,:] .< 1e-3)];
mesh[1].boundarycondition = [1 1 1 1]; # Set boundary condition for each disjoint boundary


text:
struct NamingError <: Exception
    invalid_name::AbstractString
    name_of::String
    reason::String
end

struct FieldRequired <: Exception
    field::String
    structname::String
end

NamingError(a::AbstractString, nameof::String) = NamingError(a, nameof, "it may be too long or contain invalid characters.")

Base.showerror(io::IO, e::NamingError) = print(io, "Name $(e.invalid_name) is an invalid name for `$(e.name_of)` because `$(e.reason)`.")
Base.showerror(io::IO, e::FieldRequired) = print(io, "Field `$(e.field)` is required to construct a `$(e.structname)`.")


text:
export GuildCreate,
    GuildUpdate,
    GuildDelete,
    GuildBanAdd,
    GuildBanRemove,
    GuildEmojisUpdate,
    GuildIntegrationsUpdate,
    GuildMemberAdd,
    GuildMemberRemove,
    GuildMemberUpdate,
    GuildMembersChunk,
    GuildRoleCreate,
    GuildRoleUpdate,
    GuildRoleDelete

"""
Sent when a new [`Guild`](@ref) is created.
"""
struct GuildCreate <: AbstractEvent
    guild::Guild
end
@boilerplate GuildCreate :docs :mock
GuildCreate(; kwargs...) = GuildCreate(Guild(; kwargs...))
GuildCreate(d::Dict{Symbol, Any}) = GuildCreate(; d...)

"""
Sent when a [`Guild`](@ref) is updated.
"""
struct GuildUpdate <: AbstractEvent
    guild::Guild
end
@boilerplate GuildUpdate :docs :mock
GuildUpdate(; kwargs...) = GuildUpdate(Guild(; kwargs...))
GuildUpdate(d::Dict{Symbol, Any}) = GuildUpdate(; d...)

"""
Sent when a [`Guild`](@ref) is deleted.
"""
struct GuildDelete <: AbstractEvent
    guild::AbstractGuild  # Supposed to be an UnavailableGuild.
end
@boilerplate GuildDelete :docs :mock
GuildDelete(; kwargs...) = GuildDelete(AbstractGuild(; kwargs...))
GuildDelete(d::Dict{Symbol, Any}) = GuildDelete(; d...)

"""
Sent when a [`User`](@ref) is banned from a [`Guild`](@ref).
"""
struct GuildBanAdd <: AbstractEvent
    guild_id::Snowflake
    user::User
end
@boilerplate GuildBanAdd :constructors :docs :mock

"""
Sent when a [`User`](@ref) is unbanned from a [`Guild`](@ref).
"""
struct GuildBanRemove <: AbstractEvent
    guild_id::Snowflake
    user::User
end
@boilerplate GuildBanRemove :constructors :docs :mock

"""
Sent when a [`Guild`](@ref) has its [`Emoji`](@ref)s updated.
"""
struct GuildEmojisUpdate <: AbstractEvent
    guild_id::Snowflake
    emojis::Vector{Emoji}
end
@boilerplate GuildEmojisUpdate :constructors :docs :mock

"""
Sent when a [`Guild`](@ref) has its [`Integration`](@ref)s updated.
"""
struct GuildIntegrationsUpdate <: AbstractEvent
    guild_id::Snowflake
end
@boilerplate GuildIntegrationsUpdate :constructors :docs :mock

"""
Sent when a [`Member`](@ref) is added to a [`Guild`](@ref).
"""
struct GuildMemberAdd <: AbstractEvent
    guild_id::Snowflake
    member::Member
end
@boilerplate GuildMemberAdd :docs :mock
GuildMemberAdd(d::Dict{Symbol, Any}) = GuildMemberAdd(; d...)
function GuildMemberAdd(; kwargs...)
    return GuildMemberAdd(snowflake(kwargs[:guild_id]), Member(; kwargs...))
end

"""
Sent when a [`Member`](@ref) is removed from a [`Guild`](@ref).
"""
struct GuildMemberRemove <: AbstractEvent
    guild_id::Snowflake
    user::User
end
@boilerplate GuildMemberRemove :constructors :docs :mock

"""
Sent when a [`Member`](@ref) is updated in a [`Guild`](@ref).
"""
struct GuildMemberUpdate <: AbstractEvent
    guild_id::Snowflake
    roles::Vector{Snowflake}
    user::User
    nick::OptionalNullable{String}
end
@boilerplate GuildMemberUpdate :constructors :docs :mock

"""
Sent when the [`Client`](@ref) requests guild members with [`request_guild_members`](@ref).
"""
struct GuildMembersChunk <: AbstractEvent
    guild_id::Snowflake
    members::Vector{Member}
end
@boilerplate GuildMembersChunk :constructors :docs :mock

"""
Sent when a new [`Role`](@ref) is created in a [`Guild`](@ref).
"""
struct GuildRoleCreate <: AbstractEvent
    guild_id::Snowflake
    role::Role
end
@boilerplate GuildRoleCreate :constructors :docs :mock

"""
Sent when a [`Role`](@ref) is updated in a [`Guild`](@ref).
"""
struct GuildRoleUpdate <: AbstractEvent
    guild_id::Snowflake
    role::Role
end
@boilerplate GuildRoleUpdate :constructors :docs :mock

"""
Sent when a [`Role`](@ref) is deleted from a [`Guild`](@ref).
"""
struct GuildRoleDelete <: AbstractEvent
    guild_id::Snowflake
    role_id::Snowflake
end
@boilerplate GuildRoleDelete :constructors :docs :mock


text:
module _Int

include("pi_iterative_gauss-legendre_int_simple.jl")
include("pi_iterative_gauss-legendre_int_inplace.jl")

end


text:
using NonconvexMultistart
function convert_common_kwargs(opt::NonconvexMultistart.HyperoptAlg, opt_kwargs;
    callback=nothing,
    maxiters=nothing,
    maxtime=nothing,
    abstol=nothing,
    reltol=nothing)

    conv_opt_kwargs = (; opt_kwargs...)

    if !isnothing(callback)
        @warn "common callback argument is currently not used by $(opt)"
    end

    if !isnothing(maxiters)
        @info "common maxiters argument refers to how many of the potential starting points will be evaluated by $(opt)"
        conv_opt_kwargs = (; conv_opt_kwargs..., iters=maxiters)
    end

    if !isnothing(maxtime)
        @warn "common maxtime argument is currently not used by $(opt)"
    end

    if !isnothing(abstol)
        @warn "common abstol argument is currently not used by $(opt)"
    end

    if !isnothing(reltol)
        @warn "common reltol argument is currently not used by $(opt)"
    end

    return conv_opt_kwargs
end

function __create_options(opt::NonconvexMultistart.HyperoptAlg;
    opt_kwargs=nothing)

    options = !isnothing(opt_kwargs) ? NonconvexMultistart.HyperoptOptions(;opt_kwargs...) : NonconvexMultistart.HyperoptOptions()

    if isa(options.sampler, NonconvexMultistart.Hyperopt.Hyperband)
        error("$(options.sampler) is currently not support by GalacticOptim")
    end

    return options
end

function _create_options(opt::NonconvexMultistart.HyperoptAlg;
    opt_kwargs=nothing,
    sub_options=nothing,
    convergence_criteria=nothing)

    options = (; options = !isnothing(opt_kwargs) ? NonconvexMultistart.HyperoptOptions(;sub_options= __create_options(opt.sub_alg, opt_kwargs=sub_options) ,opt_kwargs...) : NonconvexMultistart.HyperoptOptions(;sub_options= __create_options(opt.sub_alg, opt_kwargs=sub_options)))

    if isa(options.options.sampler, NonconvexMultistart.Hyperopt.Hyperband)
        error("$(options.options.sampler) is currently not support by GalacticOptim")
    end

    return options
end

check_optimizer_backend(opt::NonconvexMultistart.HyperoptAlg) = false


text:
"""
```
standard_spec!(m, vint, directory_path; cdvt = vint,
    dsid = data_id(m), cdid = cond_id(m), fcast-date = Dates.lastdayofquarter(Dates.today()),
    fcast_end = DSGE.quartertodate("2030-Q1"))
```
adds the standard settings to the model object `m` for
the Brookings Phillips Curve project.
"""
function standard_spec!(m::AbstractDSGEModel, vint::String,
                        directory_path::String = dirname(@__FILE__); cdvt::String = vint,
                        dsid::Int = data_id(m),
                        cdid::Int = cond_id(m),
                        fcast_date::Dates.Date = Dates.lastdayofquarter(Dates.today()),
                        fcast_end::Date = DSGE.quartertodate("2024-Q4"),
                        four_folders_down::Bool = false,
                        save_orig::Bool = false)
    m <= Setting(:data_vintage, vint)
    if four_folders_down
        m <= Setting(:saveroot, joinpath(directory_path, "../../../../save/"))
        m <= Setting(:dataroot, joinpath(directory_path, "../../../../save/input_data"))
    else
        m <= Setting(:saveroot, joinpath(directory_path, "../../../save/"))
        m <= Setting(:dataroot, joinpath(directory_path, "../../../save/input_data"))
    end
    if save_orig
        if four_folders_down
            m <= Setting(:saveroot, joinpath(directory_path, "../../../../save_orig/"))
            m <= Setting(:dataroot, joinpath(directory_path, "../../../../save/input_data"))
        else
            m <= Setting(:saveroot, joinpath(directory_path, "../../../save_orig/"))
            m <= Setting(:dataroot, joinpath(directory_path, "../../../save/input_data"))
        end
    end
    m <= Setting(:data_id, dsid)
    m <= Setting(:cond_id, cdid)
    m <= Setting(:hours_first_observable, true)

    m <= Setting(:cond_vintage, cdvt)
    m <= Setting(:date_forecast_start,  fcast_date)
    m <= Setting(:forecast_horizons, subtract_quarters(fcast_end, fcast_date) + 1)
    m <= Setting(:impulse_response_horizons, 20)
    m <= Setting(:date_conditional_end, fcast_date)
    m <= Setting(:use_population_forecast, true)
    m <= Setting(:forecast_uncertainty_override, Nullable{Bool}(false))

    m <= Setting(:sampling_method, :SMC)
    m <= Setting(:forecast_block_size, 500)
    m <= Setting(:forecast_jstep, 1)

    if cdid == 1
        m <= Setting(:cond_semi_names, [])
        m <= Setting(:cond_full_names, [:obs_hours])
    end
end

function standard_spec200129!(m::AbstractDSGEModel, vint::String,
                        directory_path::String = dirname(@__FILE__); cdvt::String = vint,
                        dsid::Int = data_id(m),
                        cdid::Int = cond_id(m),
                        fcast_date::Dates.Date = Dates.lastdayofquarter(Dates.today()),
                        fcast_end::Date = DSGE.quartertodate("2024-Q4"),
                        four_folders_down::Bool = false)
    m <= Setting(:data_vintage, vint)
    if four_folders_down
        m <= Setting(:saveroot, joinpath(directory_path, "../../../../save/200129"))
        m <= Setting(:dataroot, joinpath(directory_path, "../../../../save/input_data"))
    else
        m <= Setting(:saveroot, joinpath(directory_path, "../../../save/200129"))
        m <= Setting(:dataroot, joinpath(directory_path, "../../../save/input_data"))
    end
    m <= Setting(:data_id, dsid)
    m <= Setting(:cond_id, cdid)
    m <= Setting(:hours_first_observable, true)

    m <= Setting(:cond_vintage, cdvt)
    m <= Setting(:date_forecast_start,  fcast_date)
    m <= Setting(:forecast_horizons, subtract_quarters(fcast_end, fcast_date) + 1)
    m <= Setting(:impulse_response_horizons, 20)
    m <= Setting(:date_conditional_end, fcast_date)
    m <= Setting(:use_population_forecast, true)
    m <= Setting(:forecast_uncertainty_override, Nullable{Bool}(false))

    m <= Setting(:sampling_method, :SMC)
    m <= Setting(:forecast_block_size, 500)
    m <= Setting(:forecast_jstep, 1)

    m <= Setting(:date_regime2_start_text, "900331", true, "reg2start",
                 "The text version to be saved of when regime 2 starts")
    m <= Setting(:preZLB, "false", true, "preZLB", "")

    if cdid == 1
        m <= Setting(:cond_semi_names, [])
        m <= Setting(:cond_full_names, [:obs_hours])
    end
end

function standard_spec200204_dsgevar!(m::AbstractDSGEModel, vint::String,
                        directory_path::String = dirname(@__FILE__); cdvt::String = vint,
                        dsid::Int = data_id(m),
                        cdid::Int = cond_id(m),
                        fcast_date::Dates.Date = Dates.lastdayofquarter(Dates.today()),
                        fcast_end::Date = DSGE.quartertodate("2024-Q4"),
                        four_folders_down::Bool = false)
    m <= Setting(:data_vintage, vint)
    if four_folders_down
        m <= Setting(:saveroot, joinpath(directory_path, "../../../../save/"))
        m <= Setting(:dataroot, joinpath(directory_path, "../../../../save/input_data"))
    else
        m <= Setting(:saveroot, joinpath(directory_path, "../../../save/"))
        m <= Setting(:dataroot, joinpath(directory_path, "../../../save/input_data"))
    end
    m <= Setting(:data_id, dsid)
    m <= Setting(:cond_id, cdid)
    m <= Setting(:hours_first_observable, false)

    m <= Setting(:cond_vintage, cdvt)
    m <= Setting(:date_forecast_start,  fcast_date)
    m <= Setting(:forecast_horizons, subtract_quarters(fcast_end, fcast_date) + 1)
    m <= Setting(:impulse_response_horizons, 20)
    m <= Setting(:date_conditional_end, fcast_date)
    m <= Setting(:use_population_forecast, false)
    m <= Setting(:forecast_uncertainty_override, Nullable{Bool}(false))

    m <= Setting(:sampling_method, :SMC)
    m <= Setting(:forecast_block_size, 500)
    m <= Setting(:forecast_jstep, 1)

    m <= Setting(:date_regime2_start_text, "900331", true, "reg2start",
                 "The text version to be saved of when regime 2 starts")
    m <= Setting(:preZLB, "false", true, "preZLB", "")
    m <= Setting(:dsgevar, "true", true, "dsgevar", "indicates output from DSGEVAR estimation")
end

"""
```
get_κp(m::AbstractDSGEModel, θ::Matrix{Float64})
```
computes the slope of the price Phillips curve, given a
`n_draws × n_parameters` matrix.
"""
function get_κp(m::AbstractDSGEModel, θ::Matrix{Float64})
    κ = Vector{Float64}(undef, size(θ, 1))
    for i = 1:size(θ, 1)
        κ[i] = get_κp(m, θ[i, :])
    end
    return κ
end

function get_κp(m::AbstractDSGEModel, θ::Vector{Float64})
    DSGE.update!(m, θ)
    return ((1 - m[:ζ_p]*m[:β]*exp((1 - m[:σ_c])*m[:z_star]))*
            (1 - m[:ζ_p]))/(m[:ζ_p]*((m[:Φ]- 1)*m[:ϵ_p] + 1))/(1 + m[:ι_p]*m[:β]*exp((1 - m[:σ_c])*m[:z_star]))
end

"""
```
get_κw(m::AbstractDSGEModel, θ::Matrix{Float64})
```
computes the slope of the wage Phillips curve, given a
`n_draws × n_parameters` matrix.
"""
function get_κw(m::AbstractDSGEModel, θ::Matrix{Float64})
    κ = Vector{Float64}(undef, size(θ, 1))
    for i = 1:size(θ, 1)
        κ[i] = get_κw(m, θ[i, :])
    end
    return κ
end

function get_κw(m::AbstractDSGEModel, θ::Vector{Float64})
    DSGE.update!(m, θ)
    return (1 - m[:ζ_w]*m[:β]*exp((1 - m[:σ_c])*m[:z_star]))*
    (1 - m[:ζ_w])/(m[:ζ_w]*((m[:λ_w] - 1)*m[:ϵ_w] + 1))/(1 + m[:β]*exp((1 - m[:σ_c])*m[:z_star]))
end

"""
```
draw_κ(m::AbstractDSGEModel; N::Int = 100000,
       kernel_density::Bool = true) where {S <: Real}
```
draws from the implied prior for κp and κw, given priors
for the parameters which specify these slopes (e.g. ζ_p).
"""
function draw_κ(m::AbstractDSGEModel; N::Int = 100000,
                kernel_density::Bool = true) where {S <: Real}
    paras = zeros(N, length(m.parameters))
    for i = 1:N
        success = false
        while !success
            try
                paras[i, :] = rand(m.parameters, 1)
            catch e
                continue
            end
            success = true
        end
    end
    κp = get_κp(m, paras)
    κw = get_κw(m, paras)

    if kernel_density
        kdep = kde(κp)
        kdew = kde(κw)

        return kdep, kdew
    else

        return κp, κw
    end
end


text:
export  setMRSconfiguration!,
        setMRSruns!,
        setMRSdatabaseName!,
        setMRSmaxThreads!

function setMRSconfiguration!( mrs::MRS, mpSim::MPsim )
    mrs.mpSim = deepcopy( mpSim )
    return true
end  # setMRSconfiguration!( mrs, sim )

function setMRSruns!( mrs::MRS, nRuns::Integer )
    if nRuns < 0
        return false
    end  # if nRuns < 0

    mrs.nRuns = nRuns
    return true
end  # setMRSruns!( mrs, nruns )

function setMRSdatabaseName!( mrs::MRS, dbName::AbstractString="" )
    if dbName == ""
        mrs.resultsDBname = ""
        mrs.resultsDB = SQLite.DB()
        return true
    end  # if dbName == ""

    # Ensure the folder exists.
    if !ispath( dirname( dbName ) )
        mkpath( dirname( dbName ) )
    end  # if !ispath( dirname( dbName ) )

    mrs.resultsDBname = string( dbName )
    mrs.resultsDB = SQLite.DB( mrs.resultsDBname )
    return true
end  # setMRSdatabaseName!( mrs, dbName )

function setMRSmaxThreads!( mrs::MRS, maxThreads::Integer )
    if maxThreads <= 0
        return false
    end  # if maxThreads <= 0

    mrs.maxThreads = maxThreads
    return true
end  # setMRSmaxThreads!( mrs, maxThreads )

verifySimulation!( mrs::MRS ) = verifySimulation!( mrs.mpSim )


text:
@testset "targetmachine" begin

host_triple = triple()
host_t = Target(host_triple)

let
    tm = TargetMachine(host_t, host_triple)
    dispose(tm)
end

TargetMachine(host_t, host_triple) do tm
    @test target(tm) == host_t
    @test triple(tm) == host_triple
    @test cpu(tm) == ""
    @test features(tm) == ""
    asm_verbosity!(tm, true)

    # emission
    Context() do ctx
    Builder(ctx) do builder
    LLVM.Module("SomeModule", ctx) do mod
        ft = LLVM.FunctionType(LLVM.VoidType(ctx))
        fn = LLVM.Function(mod, "SomeFunction", ft)

        entry = BasicBlock(fn, "entry")
        position!(builder, entry)

        ret!(builder)

        asm = String(convert(Vector{UInt8}, emit(tm, mod, LLVM.API.LLVMAssemblyFile)))

        mktemp() do path, io
            emit(tm, mod, LLVM.API.LLVMAssemblyFile, path)
            @test asm == read(path, String)
        end

        @test_throws LLVMException emit(tm, mod, LLVM.API.LLVMAssemblyFile, "/")
    end
    end
    end

    Context() do ctx
    LLVM.Module("SomeModule", ctx) do mod
        FunctionPassManager(mod) do fpm
            add_transform_info!(fpm, tm)
            add_library_info!(fpm, triple(tm))
        end
        ModulePassManager() do mpm
            add_transform_info!(mpm, tm)
            add_library_info!(mpm, triple(tm))
        end
    end
    end

    DataLayout(tm)
end

end


text:
using SPICE

file = "gm_de431.tpc"

ids = Int[]
gms = Float64[]
re = r"BODY(?<id>[0-9]+)_GM\s+= \(\s+(?<gm>[0-9\.\+\-E]+)\s+\)"
lines = open(readlines, file)
for line in lines
    s = string(line)
    if occursin(re, s)
        m = match(re, s)
        push!(ids, parse(Int, m["id"]))
        push!(gms, parse(Float64, m["gm"]))
    end
end

open("gm.jl", "w") do f
    write(f, "# Automatically generated by $(basename(@__FILE__)), do not edit!\n\n")
    write(f, "grav_param(::Float64, ::SolarSystemBarycenter) = $(gms[findfirst(ids .== 10)])\n")
    write(f, "grav_param(::SolarSystemBarycenter) = grav_param(Float64, SolarSystemBarycenter())\n")
    for (id, gm) in zip(ids, gms)
        try
            body = replace(titlecase(bodc2n(id)), r"\s"=>"")
            body = body == "Moon" ? "Luna" : body
            write(f, "grav_param(::Type{Float64}, ::$body) = $gm\n")
            write(f, "grav_param(::$body) = grav_param(Float64, $body())\n")
        catch
            @show id
        end
    end
end


text:
using Test, RangeEnclosures
using RangeEnclosures: relative_precision, available_solvers

using DynamicPolynomials: @polyvar

include("univariate.jl")
include("multivariate.jl")


text:
function generate_posit_h(io, posit_defs)
  #generates "posit.h" based on the posit_definitions
  write(io, """
    #ifndef __POSIT_H
    #define __POSIT_H

    #include <stdint.h>
    #include <math.h>
    #include <stdbool.h>
    """)

  #define a data structure that holds the environment settings and extern it.
  write(io, """

  typedef struct {
    bool nanmode;
    bool underflows;
  } environment_t;

  extern "C" environment_t POSIT_ENV;

  extern "C" void set_nanmode(bool nanmode);
  extern "C" void set_underflow(bool underflows);
  extern "C" bool get_nanmode();
  extern "C" bool get_underflow();

  """)

  #create literal constants.
  for n in sort(collect(keys(posit_defs)))
    T = top_bits(n)
    write(io,"""
    #define P$(n)INF      $(c_literal(T))
    #define P$(n)MAXREAL  $(c_literal(T - one(T)))
    #define P$(n)MINREAL  $(c_literal(T | one(T)))
    #define P$(n)POSSMALL $(c_literal(one(T)))
    #define P$(n)NEGSMALL $(c_literal(-one(T)))
    #define P$(n)ZER      $(c_literal(zero(T)))
    """)
  end
  write(io, "\n")

  for n in sort(collect(keys(posit_defs)))
    for es in posit_defs[n]
      write(io, "typedef union {uint$(n)_t udata; int$(n)_t sdata;} p$(n)e$(es)_t ;\n")
    end
    write(io, "\n")
  end

  for n in sort(collect(keys(posit_defs)))
    for es in posit_defs[n]
      write(io, "typedef struct {uint64_t *cache; uint64_t exponent; uint64_t leading_ones;} fdp_cache$(n)e$(es)_t;\n")
    end
    write(io, "\n")
  end

  write(io, "#endif")  #to close out the header guard #define

end


text:
""" "Each Tuple corresponds to (con_name, internal_index, moi_index)"""
function get_all_constraint_index(op_problem::OperationsProblem)
    con_index = Vector{Tuple{Symbol, Int64, Int64}}()
    for (key, value) in op_problem.canonical.constraints
        for (idx, constraint) in enumerate(value)
            moi_index = JuMP.optimizer_index(constraint);
            push!(con_index, (key, idx, moi_index.value))
        end
    end

    return con_index
end

""" "Each Tuple corresponds to (con_name, internal_index, moi_index)"""
function get_all_var_index(op_problem::OperationsProblem)
    var_index = Vector{Tuple{Symbol, Int64, Int64}}()
    for (key, value) in op_problem.canonical.variables
        for (idx, variable) in enumerate(value)
            moi_index = JuMP.optimizer_index(variable);
            push!(var_index, (key, idx, moi_index.value))
        end
    end

    return var_index
end

function get_con_index(op_problem::OperationsProblem, index::Int64)

    for i in get_all_constraint_index(op_problem::OperationsProblem)
        if i[3] == index
            return op_problem.canonical.constraints[i[1]].data[i[2]]
        end
    end

    @info "Index not found"

    return

end

function get_var_index(op_problem::OperationsProblem, index::Int64)

    for i in get_all_var_index(op_problem::OperationsProblem)
        if i[3] == index
            return op_problem.canonical.variables[i[1]].data[i[2]]
        end
    end

    @info "Index not found"

    return

end


text:
# Loads SMB results from Matthias' GloGEM model
#
# Format is Elevation bands vs year total-balance & winter-balance
# Elev  1980  1981  1982  1983  1984  1985  1986  1987  1988  1989  1990  1991  1992  1993  1994  1995  1996  1997  1998  1999  2000  2001  2002  2003  2004  2005  2006  2007  2008  2009  2010  2011  2012  1980  1981  1982  1983  1984  1985  1986  1987  1988  1989  1990  1991  1992  1993  1994  1995  1996  1997  1998  1999  2000  2001  2002  2003  2004  2005  2006  2007  2008  2009  2010  2011  2012
# 2495  -0.23  -0.55  -0.79  -0.61  -0.21  -0.90  -1.19  -0.82  -1.17  -0.79  -1.20  -0.87  -1.45  -0.91  -1.30  -0.35  -1.18  -0.97  -1.44  -1.30  -1.41  -0.28  -1.72  -2.36  -0.59  -1.93  -1.43  -2.05  -1.33  -1.73  -0.98  -2.06  -1.68   0.72   0.71   0.76   0.72   0.55   0.53   0.68   0.50   0.64   0.44   0.53   0.52   0.70   0.70   0.71   0.65   0.06   0.62   0.47   0.78   0.65   0.93   0.32   0.62   0.71   0.23   0.47   0.11   0.56   0.59   0.53   0.48   0.54
# 2505  -0.21  -0.52  -0.76  -0.59  -0.18  -0.87  -1.16  -0.80  -1.13  -0.76  -1.17  -0.85  -1.42  -0.88  -1.27  -0.32  -1.14  -0.94  -1.41  -0.91  -1.38  -0.26  -1.68  -2.33  -0.57  -1.52  -1.40  -2.01  -1.30  -1.70  -0.95  -1.73  -1.65   0.72   0.71   0.76   0.73   0.56   0.54   0.69   0.51   0.65   0.45   0.53   0.53   0.70   0.71   0.72   0.66   0.07   0.63   0.48   0.78   0.65   0.94   0.32   0.62   0.72   0.24   0.47   0.12   0.56   0.59   0.53   0.48   0.54
# ...

dir2reg  = Dict(
    :Alaska=>1,
    :WesternCanada=>2,
    :ArcticCanadaNorth=>3,
    :ArcticCanadaSouth=>4,
    :Greenland=>5,
    :Iceland=>6,
    :Svalbard=>7, # and Jan Mayen
    :Scandinavia=>8,
    :RussianArctic=>9,
    :NorthAsia=>10,
    :CentralEurope=>11,
    :CaucasusMiddleEast=>12,
    :CentralAsia=>13,
    :SouthAsiaWest=>14,
    :SouthAsiaEast=>15,
    :LowLatitudes=>16,
    :SouthernAndes=>17,
    :NewZealand=>18,
    :SubAntarctic=>19
)

reg2dir = Dict(r=>n for (n,r) in dir2reg)

function load!(ds::DataSet{BdotData, GloGEMLoader}, dem, glaciermask, target_proj,
               outline, cropbox, pp::Phys)
    @unpack files, opts, daterange = ds
    @unpack sigma = opts
    make2D = get(opts, :make2D, true)

    rgi = getrgi(outline.gid)
    vers,reg,nr = parse_rgi(rgi)
    if vers!=v"6"
        error("GloGEM runs on RGI version 6 but here it is version $vers")
    end
    # make filename
    #base = "/home/mauro/netfs/vaw/root/scratch_net/sermeq/mhuss/r6_global_results"
    #ff = joinpath(base, string(reg2dir[reg]), "PAST/mb_elevation", "belev_$(VAWTools.int2str5(nr)).dat")
    ff = ds.files[1]
    #
    out = get_cache_ram!(ds, :bdot, ()->readdlm(ff, F, header=true))
    # only keep first half as second half is winter-balance:
    years = map(x->parse(Int,x),out[2][2:end÷2])
    bdots = out[1][:,2:end÷2] # values at elevation bands and years
    for i in eachindex(bdots)
        if bdots[i]==-99
            bdots[i]=FILL
        end
    end

    ele_ = out[1][:,1] # middle of elevation band
    # convert to range
    st = diff(ele_)[1]
    ele = ele_[1]:st:ele_[end]
    @assert collect(ele)==ele_

    # average over some years
    if daterange==(Date(),Date()) # use all
        bdot = squeeze(mean(bdots,2),2)
    else
        inds = Dates.year(daterange[1]).<=years.<=Dates.year(daterange[2])
        bdot = squeeze(mean(bdots[:,inds],2),2)
    end

    # scale to m ice / a
    scale_m_ice!(bdot, pp)
    vals = Gridded1d(ele,bdot)

    # fill gaps
    fillgaps!(vals, 3)

    if make2D
        # make a 2D bdot field:
        bands, bandi = VAWTools.bin_grid(dem, vals.x, glaciermask) #, min_bands=0)
        bdot_2D = VAWTools.map_back_to_2D(size(dem.v), bandi, vals.v)
        return SomeData(vals=Gridded(dem.x, dem.y, bdot_2D), sigma=sigma)
    else
        return SomeData(vals=vals, sigma=sigma)
    end
end


text:
function set_default_values(cur_player::AbstractPlayer)

  cur_player.grid.height = 0
  fill!(cur_player.grid.table, "")

  empty!(cur_player.bag.pieces)

  if cur_player.round.is_keeping_score
    Round(cur_player)
  end

  cur_player.hold = Nullable{AbstractPiece}()
  cur_player.piece = Nullable{AbstractPiece}()
  cur_player.shadow = Piece(cur_player)

  cur_player.level = 1
  cur_player.score = 0

  cur_player.lines = 0
  cur_player.combo = 0

  cur_player.glues = 0

end


text:
module Basis

include("abstractomega.jl")

end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Blosc_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Blosc")
JLLWrappers.@generate_main_file("Blosc", UUID("0b7ba130-8d10-5ba8-a3d6-c5182647fed9"))
end  # module Blosc_jll


text:
using MyGraph

cap = -1
T = [9 1 4 9 0 4 8 9 0 1] # map(parse, split(readline(STDIN), r"(\s+)"))
G = DictGraph()

# search for the capital vertex
# T[Q] = Q then Q is the capital

for Q = 1:lastindex(T)
	Q - 1 == T[Q] && (global cap = T[Q])
end

if cap == -1
    error("no capital exists")
else
    println("The capital is $cap")
    for i = 1:lastindex(T)
        addEdge!(G, i - 1, T[i])
    end
end

D = DFS(G)
B = BFS(G, cap)

R = [0 for _ = 1:B.count]

for k in keys(B)
    R[B[k][2] + 1] += 1
end

println("R = $(R[2:end])")


text:
"""
Analysis data of Monte Carlo simulation
"""
@with_kw mutable struct MCAnalysis
    acc_rate::Float64 = 0.
    prop_local::Int = 0
    acc_local::Int = 0
    acc_rate_global::Float64 = 0.
    prop_global::Int = 0
    acc_global::Int = 0
end

"""
Parameters of Monte Carlo
"""
@with_kw struct MCParameters
    global_moves::Bool = false
    global_rate::Int = 5
    thermalization::Int = 0 # number of thermalization sweeps
    sweeps::Int = 1000 # number of sweeps (after thermalization)
    measure_rate::Int = 1 # measure at every nth sweep (after thermalization)
    print_rate::Int = 1000

    beta::Float64
end

"""
Monte Carlo simulation
"""
mutable struct MC{M<:Model, C, RT<:AbstractRecorder} <: MonteCarloFlavor
    model::M
    conf::C
    configs::RT
    last_sweep::Int

    thermalization_measurements::Dict{Symbol, AbstractMeasurement}
    measurements::Dict{Symbol, AbstractMeasurement}
    p::MCParameters
    a::MCAnalysis

    MC{M,C,RT}() where {M,C,RT} = new()
    function MC(m::M, c::C, cs::RT, ls, tm, mm, p, a) where {M, C, RT}
        new{M, C, RT}(m, c, cs, ls, tm, mm, p, a)
    end
end

"""
    MC(m::M; kwargs...) where M<:Model

Create a Monte Carlo simulation for model `m` with keyword parameters `kwargs`.
"""
function MC(m::M;
        seed::Int=-1,
        thermalization_measurements = Dict{Symbol, AbstractMeasurement}(),
        measurements = :default,
        last_sweep = 0,
        measure_rate = 1,
        recording_rate = measure_rate,
        recorder = Discarder(MC, M, recording_rate),
        kwargs...
    ) where M<:Model

    conf = rand(MC, m)
    mc = MC{M, typeof(conf), typeof(recorder)}()
    mc.conf = conf
    mc.model = m
    kwdict = Dict(kwargs)
    if :T in keys(kwargs)
        kwdict[:beta] = 1/kwargs[:T]
        delete!(kwdict, :T)
    end
    mc.p = MCParameters(measure_rate = measure_rate; kwdict...)
    mc.last_sweep = last_sweep
    mc.configs = recorder

    init!(
        mc, seed = seed, conf = conf,
        thermalization_measurements = thermalization_measurements,
        measurements = measurements
    )
    return mc
end

"""
    MC(m::M, params::Dict)
    MC(m::M, params::NamedTuple)

Create a Monte Carlo simulation for model `m` with (keyword) parameters
as specified in the dictionary/named tuple `params`.
"""
MC(m::Model, params::Dict{Symbol, T}) where T = MC(m; params...)
MC(m::Model, params::NamedTuple) = MC(m; params...)

# convenience
@inline beta(mc::MC) = mc.p.beta
@inline model(mc::MC) = mc.model
@inline conf(mc::MC) = mc.conf
@inline last_sweep(mc::MC) = mc.last_sweep
@inline configurations(mc::MC) = mc.configs
@inline lattice(mc::MC) = lattice(mc.model)
@inline parameters(mc::MC) = merge(parameters(mc.p), parameters(mc.model))
@inline function parameters(p::MCParameters)
    (beta = p.beta, thermalization = p.thermalization, sweeps = p.sweeps)
end

# cosmetics
Base.summary(mc::MC) = "MC simulation of $(summary(mc.model))"
function Base.show(io::IO, mc::MC)
    print(io, "Monte Carlo simulation\n")
    print(io, "Model: ", mc.model, "\n")
    print(io, "Beta: ", round(beta(mc), sigdigits=3), " (T ≈ $(round(1/beta(mc), sigdigits=3)))\n")
    N_th_meas = length(mc.thermalization_measurements)
    N_me_meas = length(mc.measurements)
    print(io, "Measurements: ", N_th_meas + N_me_meas, " ($N_th_meas + $N_me_meas)")
end
Base.show(io::IO, m::MIME"text/plain", mc::MC) = print(io, mc)

function ConfigRecorder(mc::MC, model::Model, rate = 10)
    ConfigRecorder{typeof(compress(mc, model, conf(mc)))}(rate)
end

# implement MonteCarloFlavor interface
"""
    init!(mc::MC[; seed::Real=-1])

Initialize the Monte Carlo simulation `mc`.
If `seed !=- 1` the random generator will be initialized with `Random.seed!(seed)`.
"""
function init!(mc::MC;
        seed::Real=-1,
        conf=rand(MC, mc.model),
        thermalization_measurements = Dict{Symbol, AbstractMeasurement}(),
        measurements = :default
    )
    seed == -1 || Random.seed!(seed)
    mc.conf = conf
    mc.a = MCAnalysis()

    mc.thermalization_measurements = thermalization_measurements
    if measurements isa Dict{Symbol, AbstractMeasurement}
        mc.measurements = measurements
    elseif measurements == :default
        mc.measurements = default_measurements(mc, mc.model)
    else
        @warn(
            "`measurements` should be of type Dict{Symbol, AbstractMeasurement}" *
            ", but are $(typeof(measurements)). No measurements have been set."
        )
        mc.measurements = Dict{Symbol, AbstractMeasurement}()
    end
    init!(mc, mc.model)

    nothing
end

function resume_init!(mc::MC)
    mc.a = MCAnalysis()
    init!(mc, mc.model)
    nothing
end

"""
    run!(mc::MC[; kwargs...])

Runs the given Monte Carlo simulation `mc`. Returns true if the run finished and
false if it cancelled early to generate a resumable save-file.

### Keyword Arguments:
- `verbose = true`: If true, print progress messaged to stdout.
- `thermalization`: Number of thermalization sweeps. Uses the value passed to
`DQMC` by default.
- `sweeps`: Number of measurement sweeps. Uses the value passed to `DQMC` by
default.
- `safe_every::TimePeriod`: Set the interval for regularly scheduled saves.
- `safe_before::Date`: If this date is passed, `run!` will generate a resumable
save file and exit
- `grace_period = Minute(5)`: Buffer between the current time and `safe_before`.
The time required to generate a save file should be included here.
- `resumable_filename`: Name of the resumable save file. The default is based on
`safe_before`.
- `overwrite = false`: If set to true a file with the same name as
`resumable_filename` will be overwritten. (This will create a temporary backup)

See also: [`resume!`](@ref)
"""
@bm function run!(
        mc::MC;
        verbose::Bool = true,
        sweeps::Int = mc.p.sweeps,
        thermalization = mc.p.thermalization,
        safe_before::TimeType = now() + Year(100),
        safe_every::TimePeriod = Hour(10000),
        grace_period::TimePeriod = Minute(5),
        resumable_filename::String = "resumable_$(Dates.format(safe_before, "d_u_yyyy-HH_MM")).jld",
        overwrite = false
    )

    do_th_measurements = !isempty(mc.thermalization_measurements)
    do_me_measurements = !isempty(mc.measurements)
    !do_me_measurements && @warn(
        "There are no measurements set up for this simulation!"
    )

    # Update number of sweeps
    if (mc.p.thermalization != thermalization) || (mc.p.sweeps != sweeps)
        verbose && println("Rebuilding DQMCParameters with new number of sweeps.")
        p = MCParameters(
            mc.p.global_moves,
            mc.p.global_rate,
            thermalization,
            sweeps,
            mc.p.measure_rate,
            mc.p.print_rate,
            mc.p.beta
        )
        mc.p = p
    end
    total_sweeps = sweeps + thermalization

    start_time = now()
    last_checkpoint = now()
    max_sweep_duration = 0.0
    verbose && println("Started: ", Dates.format(start_time, "d.u yyyy HH:MM"))

    _time = time()
    do_th_measurements && prepare!(mc.thermalization_measurements, mc, mc.model)
    for i in mc.last_sweep+1:total_sweeps
        sweep(mc)
        mc.last_sweep = i

        if mc.p.global_moves && mod(i, mc.p.global_rate) == 0
            mc.a.prop_global += 1
            mc.a.acc_global += global_move(mc, mc.model, conf(mc))
        end

        # For optimal performance whatever is most likely to fail should be
        # checked first.
        if i <= thermalization && iszero(mod(i, mc.p.measure_rate)) && do_th_measurements
            measure!(mc.thermalization_measurements, mc, mc.model, i)
        end
        if (i == thermalization+1)
            do_th_measurements && finish!(mc.thermalization_measurements, mc, mc.model)
            do_me_measurements && prepare!(mc.measurements, mc, mc.model)
        end
        if i > thermalization && iszero(mod(i, mc.p.measure_rate)) && do_me_measurements
            measure!(mc.measurements, mc, mc.model, i)
        end
        (i > thermalization) && push!(mc.configs, mc, mc.model, i)

        print_rate = mc.p.print_rate
        if print_rate != 0 && iszero(mod(i, print_rate))
            mc.a.acc_rate /= print_rate
            mc.a.acc_rate_global /= print_rate / mc.p.global_rate
            sweep_dur = (time() - _time)/print_rate
            max_sweep_duration = max(max_sweep_duration, sweep_dur)
            if verbose
                println("\t", i)
                @printf("\t\tsweep dur: %.3fs\n", sweep_dur)
                @printf("\t\tacc rate (local) : %.1f%%\n", mc.a.acc_rate*100)
                if mc.p.global_moves
                    @printf("\t\tacc rate (global): %.1f%%\n", mc.a.acc_rate_global*100)
                    @printf(
                        "\t\tacc rate (global, overall): %.1f%%\n",
                        mc.a.acc_global/mc.a.prop_global*100
                    )
                end
            end

            mc.a.acc_rate = 0.0
            mc.a.acc_rate_global = 0.0
            flush(stdout)
            _time = time()
        end

        if safe_before - now() < Millisecond(grace_period) +
                Millisecond(round(Int, 2e3max_sweep_duration))

            println("Early save initiated for sweep #$i.\n")
            verbose && println("Current time: ", Dates.format(now(), "d.u yyyy HH:MM"))
            verbose && println("Target time:  ", Dates.format(safe_before, "d.u yyyy HH:MM"))
            resumable_filename = save(resumable_filename, mc, overwrite=overwrite)
            verbose && println("\nEarly save finished")

            return false
        elseif (now() - last_checkpoint) > safe_every
            verbose && println("Performing scheduled save.")
            last_checkpoint = now()
            save(resumable_filename, mc, overwrite = overwrite, rename = false)
        end
    end
    do_me_measurements && finish!(mc.measurements, mc, mc.model)

    mc.a.acc_rate = mc.a.acc_local / mc.a.prop_local
    mc.a.acc_rate_global = mc.a.acc_global / mc.a.prop_global

    end_time = now()
    if verbose
        println("Ended: ", Dates.format(end_time, "d.u yyyy HH:MM"))
        @printf("Duration: %.2f minutes", (end_time - start_time).value/1000. /60.)
        println()
    end

    return true
end

"""
    sweep(mc::MC)

Performs a sweep of local moves.
"""
@inline @bm function sweep(mc::MC)
    c = conf(mc)
    m = model(mc)
    β = beta(mc)
    @inbounds for i in eachindex(c)
        ΔE, passthrough = propose_local(mc, m, i, c)
        mc.a.prop_local += 1
        # Metropolis
        if ΔE <= 0 || rand() < exp(- β*ΔE)
            accept_local!(mc, m, i, c, ΔE, passthrough)
            mc.a.acc_rate += 1.0
            mc.a.acc_local += 1
        end
    end
    mc.a.acc_rate /= length(c)
    nothing
end

"""
    replay(mc::MC[, configurations::Iterable = mc.configs; kwargs...])

Replays previously generated configurations and measures observables along the
way.

### Keyword Arguments:
- `verbose = true`: If true, print progress messaged to stdout.
- `safe_every::TimePeriod`: Set the interval for regularly scheduled saves.
- `safe_before::Date`: If this date is passed, `replay!` will generate a
resumable save file and exit
- `grace_period = Minute(5)`: Buffer between the current time and `safe_before`.
The time required to generate a save file should be included here.
- `filename`: Name of the save file. The default is based on `safe_before`.
- `measure_rate = 1`: Rate at which measurements are taken. Note that this is
based on the recorded configurations, not actual sweeps.
"""
function replay!(
        mc::MC, configurations = mc.configs;
        verbose::Bool=true,
        safe_before::TimeType = now() + Year(100),
        safe_every::TimePeriod = Hour(10000),
        grace_period::TimePeriod = Minute(5),
        filename::String = "resumable_$(Dates.format(safe_before, "d_u_yyyy-HH_MM")).jld",
        overwrite = false,
        measure_rate = 1
    )
    if isempty(configurations)
        println("Nothin to replay (configurations empty). Exiting")
        return true
    end

    start_time = now()
    last_checkpoint = now()
    max_sweep_duration = 0.0
    verbose && println("Started: ", Dates.format(start_time, "d.u yyyy HH:MM"))

    # Check for measurements
    !isempty(mc.thermalization_measurements) && @warn(
        "There is no thermalization process in a replayed simulation."
    )
    isempty(mc.measurements) && @warn(
        "There are no measurements set up for this simulation!"
    )

    mc.p = MCParameters(
        global_moves = mc.p.global_moves,
        global_rate = mc.p.global_rate,
        thermalization = mc.p.thermalization,
        sweeps = mc.p.sweeps,
        measure_rate = measure_rate,
        print_rate = mc.p.print_rate,
        beta = mc.p.beta
    )
    mc.conf = copy(decompress(mc, mc.model, configurations[1]))
    resume_init!(mc)

    _time = time()
    verbose && println("\n\nReplaying measurement stage - ", length(configurations))
    prepare!(mc.measurements, mc, mc.model)
    for i in mc.last_sweep+1:mc.p.measure_rate:length(configurations)
        mc.conf .= decompress(mc, mc.model, configurations[i])
        mc.last_sweep = i

        energy(mc, mc.model, mc.conf)
        measure!(mc.measurements, mc, mc.model, i)

        if mc.p.print_rate != 0 && iszero(mod(i, mc.p.print_rate))
            sweep_dur = (time() - _time)/print_rate
            max_sweep_duration = max(max_sweep_duration, sweep_dur)
            if verbose
                println("\t", i)
                @printf("\t\tsweep dur: %.3fs\n", sweep_dur)
            end
            flush(stdout)
            _time = time()
        end

        if safe_before - now() < Millisecond(grace_period) +
                Millisecond(round(Int, 2e3max_sweep_duration))

            println("Early save initiated for sweep #$i.\n")
            verbose && println("Current time: ", Dates.format(now(), "d.u yyyy HH:MM"))
            verbose && println("Target time:  ", Dates.format(safe_before, "d.u yyyy HH:MM"))
            filename = save(filename, mc, overwrite = overwrite)
            verbose && println("\nEarly save finished")

            return false
        elseif (now() - last_checkpoint) > safe_every
            verbose && println("Performing scheduled save.")
            last_checkpoint = now()
            save(resumable_filename, mc, overwrite = overwrite, rename = false)
        end
    end
    finish!(mc.measurements, mc, mc.model)

    end_time = now()
    verbose && println("Ended: ", Dates.format(end_time, "d.u yyyy HH:MM"))
    verbose && @printf("Duration: %.2f minutes", (end_time - start_time).value/1000. /60.)

    return true
end

#     save_mc(filename, mc, entryname)
#
# Saves (minimal) information necessary to reconstruct a given `mc::MC` to a
# JLD-file `filename` under group `entryname`.
#
# When saving a simulation the default `entryname` is `MC`
function save_mc(file::JLDFile, mc::MC, entryname::String="MC")
    write(file, entryname * "/VERSION", 1)
    write(file, entryname * "/tag", "MC")
    write(file, entryname * "/type", typeof(mc))
    save_parameters(file, mc.p, entryname * "/parameters")
    write(file, entryname * "/conf", mc.conf)
    _save(file, mc.configs, entryname * "/configs")
    write(file, entryname * "/last_sweep", mc.last_sweep)
    save_measurements(file, mc, entryname * "/Measurements")
    save_model(file, mc.model, entryname * "/Model")
    nothing
end

#     load_mc(data, ::Type{<: MC})
#
# Loads a MC from a given `data` dictionary produced by `JLD.load(filename)`.
function _load(data, ::Val{:MC})
    if !(data["VERSION"] == 1)
        throw(ErrorException("Failed to load $T version $(data["VERSION"])"))
    end
    p = _load(data["parameters"], Val(:MCParameters))
    conf = data["conf"]
    configs = _load(data["configs"], to_tag(data["configs"]))
    last_sweep = data["last_sweep"]
    model = _load(data["Model"], to_tag(data["Model"]))

    measurements = _load(data["Measurements"], Val(:Measurements))

    MC(
        model, conf, configs, last_sweep,
        measurements[:TH], measurements[:ME],
        p, MCAnalysis()
    )
end

#   save_parameters(file::JLDFile, p::MCParameters, entryname="Parameters")
#
# Saves (minimal) information necessary to reconstruct a given
# `p::MCParameters` to a JLD-file `filename` under group `entryname`.
#
# When saving a simulation the default `entryname` is `MC/Parameters`
function save_parameters(file::JLDFile, p::MCParameters, entryname::String="Parameters")
    write(file, entryname * "/VERSION", 1)
    write(file, entryname * "/tag", "MCParameters")

    write(file, entryname * "/global_moves", Int(p.global_moves))
    write(file, entryname * "/global_rate", p.global_rate)
    write(file, entryname * "/thermalization", p.thermalization)
    write(file, entryname * "/sweeps", p.sweeps)
    write(file, entryname * "/measure_rate", p.measure_rate)
    write(file, entryname * "/print_rate", p.print_rate)
    write(file, entryname * "/beta", p.beta)

    nothing
end

#     load_parameters(data, ::Type{<: MCParameters})
#
# Loads a MCParameters object from a given `data` dictionary produced by
# `JLD.load(filename)`.
function _load(data, ::Val{:MCParameters})
    if !(data["VERSION"] == 1)
        throw(ErrorException("Failed to load $T version $(data["VERSION"])"))
    end

    MCParameters(
        Bool(data["global_moves"]),
        data["global_rate"],
        data["thermalization"],
        data["sweeps"],
        data["measure_rate"],
        data["print_rate"],
        data["beta"]
    )
end

include("MC_mandatory.jl")
include("MC_optional.jl")


text:
function DiffEqBase.__solve(prob::AbstractNoiseProblem,args...;dt=0.0,kwargs...)
  if dt == 0.0 || dt == nothing
    error("dt must be provided to simulate a noise process. Please pass dt=...")
  end
  W = deepcopy(prob.noise)
  if typeof(W) <: NoiseProcess
    if prob.seed != 0
      Random.seed!(W.rng,prob.seed)
    else
      Random.seed!(W.rng,rand(UInt64))
    end
  end
  W.curt = prob.tspan[1]
  W.dt = dt
  setup_next_step!(W,nothing,nothing)
  tType = typeof(W.curt)
  while W.curt < prob.tspan[2]
    if tType <: AbstractFloat && abs(W.curt + dt - prob.tspan[2]) < 100*eps(dt) # Correct the end due to floating point error
      dt = prob.tspan[2]- W.curt
      accept_step!(W,dt,nothing,nothing)
      W.curt = prob.tspan[2]
    else
      accept_step!(W,dt,nothing,nothing)
    end
  end
  W
end


text:
"""
    gapfill!(x t, fs; m::Bool=true)
    y = gapfill(x, t, fs; m::Bool=true)

Fill gaps in `x`, sampled at `fs`, with gap indices given by `t[:,1]` and
gap lengths in μs given by `t[:,2]`.

Specify `m=false` to fill with NaNs; else, fill with the mean of non-NaN
values in `x`.
"""
function gapfill!(x::Array{T,1}, t::Array{Int64,2}, fs::Float64; m::Bool=true) where T<: Real
  (fs == 0.0 || isempty(x)) && (return x)
  nt = size(t,1)
  ng = nt - (t[nt,2] == 0 ? 2 : 1)
  (ng == 0) && (return x)

  timegaps = t[2:end,2]                       # number of time gaps
  mx = m ? mean(skipmissing(x)) : T(NaN)      # mean or NaN
  Δ = round(Int64, 1.0/(fs*μs))               # sampling interval in μs

  # this always yields rows in w whose indices are correct output assignments
  w = t_win(t, Δ)
  nw = size(w,1)
  broadcast!(-, w, w, w[1,1])
  broadcast!(div, w, w, Δ)
  broadcast!(+, w, w, 1)
  nx = maximum(w)-minimum(w)+1

  if minimum(timegaps) ≥ 0
    resize!(x, nx)
    for i = nw:-1:1
      N = w[i,2]-w[i,1]+1
      j = t[i,1]
      copyto!(x, w[i,1], x, j, N)
      if i > 1
        fill_s = w[i-1,2]+1
        fill_e = w[i,1]-1
        x[fill_s:fill_e] .= mx
      end
    end
  else
    # Fix for issue #29
    x1 = Array{T,1}(undef, nx)

    # with negative gaps, source start indices are in t[:,1]
    for i in 1:nw
      N = w[i,2]-w[i,1]+1
      j = t[i,1]
      copyto!(x1, w[i,1], x, j, N)
      if i > 1
        fill_s = w[i-1,2]+1
        fill_e = w[i,1]-1
        (fill_e > fill_s) && (x1[fill_s:fill_e] .= mx)
      end
    end
    resize!(x, nx)
    copyto!(x, 1, x1, 1, nx)
  end
  return nothing
end


text:
circles = html_div() do
    html_div(id="fields") do
        html_div(className="container-slider") do
            html_label("Samples: "),
            html_div(
                dcc_slider(
                        id = "n_samples",
                        min = 1,
                        max = 10000,
                        marks = Dict(i => "$i" for i = 0:1000:10000),
                        value = 100,
                ),
            style = (width = "70%", display = "inline-block"),
            className = "slider")
        end,
        html_div(className="container-slider") do
            html_label("Noise: "),
            html_div(
                dcc_slider(
                        id = "noise",
                        min = 0.0,
                        max = 1.0,
                        step = 0.01,
                        marks = Dict(i => "$i" for i = 0.0:0.1:1.0),
                        value = 0.1,
                ),
            style = (width = "70%", display = "inline-block"),
            className = "slider")
        end,
        html_div(className="container-slider") do
            html_label("Factor: "),
            html_div(
                dcc_slider(
                        id = "factor",
                        min = 0.0,
                        max = 1.0,
                        step = 0.01,
                        marks = Dict(i => "$i" for i = 0.0:0.1:1.0),
                        value = 0.1,
                ),
            style = (width = "70%", display = "inline-block"),
            className = "slider")
        end
    end,
    html_div(className="container-graph") do
        dcc_graph(
            id = "plot_circles",
            figure = gen_circles(1000, 0.25,0.25)
        )
    end
end

callback!(app, Output("plot_circles", "figure"),[Input("n_samples", "value"), Input("noise", "value"),Input("factor", "value")]) do n_samples, noise, factor
    gen_circles(n_samples,noise,factor)
end


text:
# Autogenerated wrapper script for Libiconv_jll for x86_64-apple-darwin
export libcharset, libiconv

JLLWrappers.@generate_wrapper_header("Libiconv")
JLLWrappers.@declare_library_product(libcharset, "@rpath/libcharset.1.dylib")
JLLWrappers.@declare_library_product(libiconv, "@rpath/libiconv.2.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libcharset,
        "lib/libcharset.1.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libiconv,
        "lib/libiconv.2.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using JSON, DataFrames, JSON2DataFrame, Test

data = JSON.parsefile("data/colbase-0.json")
data = data["results"]["bindings"]

d = deepcopy(data)
flat = JSON2DataFrame.flatten.(d)
@test flat[1]["c_type"] == "uri"
df = JSON2DataFrame.dicts_to_df(flat)
@test df[1,:c_type] == "uri"
df = JSON2DataFrame.flatten_to_df(d)
@test df[1,:c_type] == "uri"


text:
# Test cmap.jl

# Tests could be more comprehensive (of course).  In some cases all
# that is done is to execute the function in some way just to make
# sure no exceptions are thrown as a result of something breaking with
# a new version of Julia.

println("testing cmap")

# cmap: Test cmap by calling it with no arguments. This generates a
# list of all the colour maps and in the process generates them all.
cmap()

# Test that all options such as shifting do not crash
testmap = cmap("C1", shift=0.25)
testmap = cmap("C2", N=24, shift=0.25)
testmap = cmap("C2", reverse=true)
testmap = cmap("C3", shift=0.25, chromaK = 0.5)

# equalisecolourmap: Generate a colour map in Lab space with an uneven
# ramp in lightness and check that this is corrected
rgblab = "LAB"
labmap = zeros(256,3)
labmap[1:127,1] = range(0, stop=20, length=127)
labmap[128:256,1] = range(20, stop=100, length=129)
formula = "CIE76"
W = [1,0,0]
sigma = 1
rgbmap = equalisecolourmap(rgblab, labmap, formula, W, sigma)
# Convert to Nx3 array and then back to lab space. Then check that dL
# is roughly constant
labmap2 = srgb2lab(rgbmap)
#dL = gradient(labmap2[:,1])  # gradient dramas
dL = labmap2[2:end,1] - labmap2[1:end-1,1]
@test maximum(dL[2:end-1]) - minimum(dL[2:end-1]) < 1e-1

# linearrgbmap
rgbmap = linearrgbmap([1,1,1],99)
rgbmap2 = PerceptualColourMaps.RGBA2FloatArray(rgbmap)
# check middle colour is [.5 .5 .5]
@test maximum(abs.(rgbmap2[50:50,:] - [0.5 0.5 0.5])) < 1e-6

# UInt32colormap
ui32 = RGBA2UInt32(cmap("L3"))

# srgb2lab and lab2srgb
rgb = [1 0 0
       0 1 0
       0 0 1]

# lab values of rgb values above
lab = [ 53.2377   80.1029   67.2144
        87.7346  -86.1606   83.1862
        32.2966   79.2009 -107.8615]

clab = srgb2lab(rgb)
crgb = lab2srgb(clab)
@test maximum(abs.(rgb-crgb)) < 1e-1
@test maximum(abs.(lab-clab)) < 1e-1


text:
# test for comparisons and hashing

let
    x = HTMLText("test")
    y = HTMLText("test")
    @test x == y
    @test hash(x) == hash(y)
    x.text *= "1"
    @test x != y
    @test hash(x) != hash(y)
    y.text *= "1"
    @test x == y
    @test hash(x) == hash(y)
end

let
    x = HTMLElement(:div)
    y = HTMLElement(:div)
    @test x == y
    @test hash(x) == hash(y)
    push!(x, HTMLElement(:p))
    @test x != y
    @test hash(x) != hash(y)
    push!(y, HTMLElement(:p))
    @test x == y
    @test hash(x) == hash(y)
    setattr!(x,"class","test")
    @test x != y
    @test hash(x) != hash(y)
    setattr!(y,"class","test")
    @test x == y
    @test hash(x) == hash(y)
end

let
    x = HTMLDocument("html", HTMLElement(:html))
    y = HTMLDocument("html", HTMLElement(:html))
    @test x == y
    @test hash(x) == hash(y)
    x.doctype = ""
    @test x != y
    @test hash(x) != hash(y)
    y.doctype = ""
    @test x == y
    @test hash(x) == hash(y)
end


text:
seed = [23324113542]
m = 2 ^ 31 - 1
a = 16807
c = 17

"""
randint(a, b) returns a random integer in [a, b[.
"""
function randint(min :: Int64, max :: Int64) :: Int64
    value = (a *  seed[0] +  c) %  m
    seed[0] = value;
    value = (value % (max - min)) + min
    return value
end

"""Returns a random FP value in [0;1[. !! Experimental"""
function rand()
    return div(randint(0, 2 ^ 31) * (2 ^ 31) + randint(0, 2 ^ 31) + 0., (2 ^ 62) - 1.)
end

function init_seed()
    println("Please enter a new value for the seed")
    seed[0] = input_int()
end