text:
# # Noisy nonlinear regression
#
# First, let us load the relevant packages and define our problem.
# We want to interpolate a noisy version of `sin(x²/2)`.

using KernelMachines, Statistics, Plots, Random

func(x) = sin(x^2 / 2)
Random.seed!(1234) # for reproducibility
xs = 8 .* rand(100)
ys = @. func(xs) + rand()
scatter(xs, ys, color="black", xlabel="", ylabel="", primary=false)

# Now, let us train a Kernel Machine on the problem. We use the additive Gaussian
# kernel. Other than the kernel, we can choose the dimensionality of the hidden
# spaces, in this case, `(3, 3, 3)`, and the regularization coefficient, `cost=0.0005`.

krm = KernelMachineRegression(
    xs, ys;
    kernel=additivegaussiankernel,
    dims=(3, 3, 3), cost=0.0005
)
fit!(krm)
us = range(extrema(xs)...; step = 0.01)
pred_krm = predict(krm, us)
mean(abs2, pred_krm .- func.(us) .- 0.5)

# Let us compare with Kernel Ridge regression. Other than the kernel, we can choose the
# regularization coefficient, `cost=0.0005`.

kr = KernelRegression(xs, ys; cost=0.0005, kernel=gaussiankernel)
fit!(kr)
pred_kr = predict(kr, us)
mean(abs2, pred_kr .- func.(us) .- 0.5)

# We can then visualize the results. As the Gaussian kernel has only a fixed
# resolution, it works well on the part of the data of comparable spatial resolution,
# but is unable to fit the part at higher frequency.

theme(:wong)

plt = scatter(xs, ys, color="black", xlabel="", ylabel="", primary=false,
    legend=:bottomleft)

plot!(plt, us, pred_krm, label="KM", linewidth=2)
plot!(plt, us, pred_kr, label="KR", linewidth=2)
plot!(plt, us, func.(us) .+ 0.5, label="Truth", linewidth=2, color="black")
plt

# Note that the fitting procedure of `KernelMachineRegression` is done via
# [Optim.jl](https://github.com/JuliaNLSolvers/Optim.jl). One can therefore
# pass both different optimization mathods and different options, such as
# number of iterations, or required accuracy. The default method is
# `ConjugateGradient`.
#
# !!! note
#
#     The API to pass optimization method and options to Optim.jl may change
#     in the future.
#
# Here, we try a different optimization method (`GradientDescent`) and
# different options (maximum number of iterations and required relative
# precision).
#

using Optim: Options, GradientDescent, ConjugateGradient

krm_descent = KernelMachineRegression(
    xs, ys;
    kernel=additivegaussiankernel,
    dims=(3, 3, 3), cost=0.0005
)
krm_cg_approx = KernelMachineRegression(
    xs, ys;
    kernel=additivegaussiankernel,
    dims=(3, 3, 3), cost=0.0005
)
fit!(krm_descent, GradientDescent(), Options(iterations=1000))
fit!(krm_cg_approx, ConjugateGradient(), Options(f_reltol=1e-3))
pred_krm_descent = predict(krm_descent, us)
pred_krm_cg_approx = predict(krm_cg_approx, us);

# In this simple test, `ConjugateGradient` with approximate convergence
# conditions is significantly faster and gives a reasonably good solution.

plt = scatter(xs, ys, color="black", xlabel="", ylabel="", primary=false,
    legend=:bottomleft)

plot!(plt, us, pred_krm, label="CG", linewidth=2)
plot!(plt, us, pred_krm_descent, label="Descent", linewidth=2)
plot!(plt, us, pred_krm_cg_approx, label="CG approx", linewidth=2)
plot!(plt, us, func.(us) .+ 0.5, label="Truth", linewidth=2, color="black")
plt

ext: jl


text:
using Interpolations: LinearInterpolation, Throw

"""
    setup_ionization_energies([filename])

Parses the table of ionization energies and returns it as a dictionary mapping elements to
their ionization energies, `[χ₁, χ₂, χ₃]` in eV.
"""
function setup_ionization_energies(fname=joinpath(_data_dir,
                                                  "BarklemCollet2016-ionization_energies.dat"))
    open(fname, "r") do f
        d = Dict{UInt8, Vector{Float64}}()
        for line in eachline(f)
            if line[1] != '#'
                toks = split(strip(line))
                Z = parse(UInt8, toks[1])
                #the second token is the atomic symbol, which we ignore
                d[Z] = parse.(Float64, toks[3:end])
            end
        end
        d
    end
end

"""
    setup_partition_funcs()

Returns a Dict holding the default (Barklem & Collet 2016) partition functions.

The partition function is linearly interpolated from values tabulated for the paper. The partition
functions can only be evaluated at temperatures between 1e-5 K and 1e4 K. The tabulated values were
NOT computed with a cutoff principle quantum number. This can lead to a couple percent error for
some neutral species in solar type stars (I think the problem becomes worse in hotter stars). See
the paper for more details.

Note that the paper seems to suggest that we could actually use cubic interpolation. That should be
revisited in the future.
"""
function setup_partition_funcs(atoms=joinpath(_data_dir, "BarklemCollet2016-atomic_partition.dat"),
                              mols=joinpath(_data_dir, "BarklemCollet2016-molecular_partition.dat"))
    merge(read_partition_funcs.([atoms, mols])...)
end

"""
    setup_equilibrium_constants()

Returns a Dict holding the default (Barklem & Collet 2016) log equilibrium constants, which are in
pressure units.

In equation 7 in this paper, m is the reduced mass, m₁m₂/(m₁ + m₂).  This can be verified by
re-deriving the equilibrium constants from the partition functions and dissolution constants
provided by the paper.
"""
setup_equilibrium_constants() = read_partition_funcs(
                                 joinpath(_data_dir, "BarklemCollet2016-equilibrium_constants.dat"),
                                 transform=x->x+1) #convert from log(mks) to log(cgs)

"""
    functon read_atomic_partition_funcs([transform=identity])

Constructs a Dict holding tables containing partition function or equilibrium constant values across
temperature.

The optional argument, `transform` is applied to each value. It is used to perform unit conversions,
when loading equilibrium constants

"""
function read_partition_funcs(fname; transform=identity)
    temperatures = Vector{Float64}()
    data_pairs = Vector{Tuple{Species,Vector{Float64}}}()
    open(fname, "r") do f
        for line in eachline(f)
            if (length(line)>=9) && contains(line, "T [K]")
                append!(temperatures, parse.(Float64,split(strip(line[10:length(line)]))))
            elseif line[1] == '#'
                continue
            else # add entries to the dictionary
                row_entries = split(strip(line))
                species_code = popfirst!(row_entries)
                #ignore deuterium and ionized molecules
                if species_code[1:2] != "D_" && species_code[end] != '+' && species_code[end] != '-'
                    push!(data_pairs, (Species(species_code),
                                       transform.(parse.(Float64, row_entries))))
                end
            end
        end
    end

    map(data_pairs) do (species, vals)
        species, LinearInterpolation(temperatures, vals, extrapolation_bc=Throw())
    end |> Dict
end


text:
using Zygote
n_comp = 2
function g(w, s=4.0)
	return 1.0 - s*w*(1-w)
end
function loss(w, s=4.0)
	fw = w
	for n = 1:n_comp
		fw = g(fw, s)
	end
	return fw
end
grad_loss(w, s=4.0) = gradient(w -> loss(w, s), w)[1]
next(w, η, s=4.0) = (w - η*grad_loss(w, s) + 1.0) % 1.0
dnext(w, η, s=4.0) = gradient(w -> next(w, η, s), w)[1]
d2next(w, η, s=4.0) = gradient(w -> dnext(w, η, s), w)[1]

max_stars_repo_name: nishaChandramoorthy/neuralOMET

max_issues_repo_name: nishaChandramoorthy/neuralOMET

max_forks_repo_name: nishaChandramoorthy/neuralOMET


text:
using Pkg
using Metatheory
using Metatheory.EGraphs

to_sexpr_pattern(e::QuoteNode) = e.value
to_sexpr_pattern(e::Symbol) = "?$e"
function to_sexpr_pattern(e::Expr)
    @assert e.head == :call
    e1 = join([e.args[1] ;  to_sexpr_pattern.(e.args[2:end])], ' ')
    "($e1)"
end

to_sexpr(e::Symbol) = e
to_sexpr(e::Expr) = "($(join(to_sexpr.(e.args),' ')))"

function eggify(rules)
    egg_rules = []
    for rule in rules
        l = to_sexpr_pattern(rule.left)
        r = to_sexpr_pattern(rule.right)
        if rule.mode == :symbolic
            push!(egg_rules,"\tvec![rw!( \"$(rule.expr)\" ; \"$l\" => \"$r\" )]")
        elseif rule.mode == :equational
            push!(egg_rules,"\trw!( \"$(rule.expr)\" ; \"$l\" <=> \"$r\" )")
        else
            println("Unsupported Rewrite Mode")
            @assert false
        end

    end
    return join(egg_rules, ",\n")
end
##########################################
# REPLACE WITH YOUR THEORY AND QUERY HERE

theory = @theory begin
    :a * b => :a
    a * b == b * a
end

query = :(a * (b * c))

###########################################

G = EGraph( query )
@time saturate!(G, theory)
extractor = addanalysis!(G, ExtractionAnalysis, astsize)
ex = extract!(G, extractor)
println( "Best found: $ex")

rust_code =
"""
use egg::{*, rewrite as rw};
//use std::time::Duration;

fn main() {
    let rules : &[Rewrite<SymbolLang, ()>] = &vec![
    $(eggify(theory))
    ].concat();

    let start = "$(to_sexpr(query))".parse().unwrap();
    let runner = Runner::default().with_expr(&start)
        // More options here https://docs.rs/egg/0.6.0/egg/struct.Runner.html
        //.with_iter_limit(10)
        //.with_node_limit(1_000_000)
        //.with_time_limit(Duration::new(30,0))
        .run(rules);
    runner.print_report();
    let mut extractor = Extractor::new(&runner.egraph, AstSize);
    let (best_cost, best_expr) = extractor.find_best(runner.roots[0]);
    println!("best cost: {}, best expr {}", best_cost, best_expr);
}
"""
open("src/main.rs", "w") do f
    write(f, rust_code)
end


text:
using LinearAlgebra
include("approx_inequality.jl")
# export ProximableFunction
export prox, prox!

# abstract type ProximableFunction end
export IndNonnegativeBallL2

struct IndNonnegativeBallL2{R <: Real}
    r::R
    function IndNonnegativeBallL2{R}(r::R) where {R <: Real}
        if r <= 0
            error("parameter r must be positive")
        else
            new(r)
        end
    end
end

is_convex(f::IndNonnegativeBallL2) = true
is_set(f::IndNonnegativeBallL2) = true

IndNonnegativeBallL2(r::R=1.0) where {R <: Real} = IndNonnegativeBallL2{R}(r)

function (f::IndNonnegativeBallL2)(x::AbstractArray{T}) where {R <: Real, T <: RealOrComplex{R}}
    if isapprox_le(norm(x), f.r, atol=eps(R), rtol=sqrt(eps(R)))
        for k in eachindex(x)
            if x[k] < 0
                return R(Inf)
            end
        end
        return R(0)
    end
    return R(Inf)
end

function prox!(y::AbstractArray{T}, f::IndNonnegativeBallL2, x::AbstractArray{T}, gamma::R=R(1)) where {R <: Real, T <: RealOrComplex{R}}
    for k in eachindex(x)
        if x[k] < 0
            y[k] = R(0)
        else
            y[k] = x[k]
        end
    end

    scal = f.r/norm(y)
    if scal > 1
        return R(0)
    end
    for k in eachindex(y)
        y[k] = scal*y[k]
    end
    return R(0)
end

function prox(f::IndNonnegativeBallL2, x::AbstractArray{T}, gamma::R=R(1)) where {R <: Real, T <: RealOrComplex{R}}
    y = copy(x)
    for k in eachindex(x)
        if x[k] < 0
            y[k] = R(0)
        else
            y[k] = x[k]
        end
    end

    scal = f.r/norm(y)
    if scal > 1
        return y, R(0)
    end
    for k in eachindex(y)
        y[k] = scal*y[k]
    end
    return y, R(0)
end

fun_name(f::IndNonnegativeBallL2) = "indicator of an L2 norm ball and nonnegative orthant"
fun_dom(f::IndNonnegativeBallL2) = "AbstractArray{Real}, AbstractArray{Complex}"
fun_expr(f::IndNonnegativeBallL2) = "x ↦ 0 if ||x|| ⩽ r & all(0 ⩽ x), +∞ otherwise"
fun_params(f::IndNonnegativeBallL2) = "r = $(f.r)"

function prox_naive(f::IndNonnegativeBallL2, x::AbstractArray{T}, gamma::R=R(1)) where {R <: Real, T <: RealOrComplex{R}}
    y = copy(x)
    for k in eachindex(x)
        if x[k] < 0
            y[k] = R(0)
        end
    end

    normx = norm(y)
    if normx > f.r
        y = (f.r/normx)*y
    end
    return y, R(0)
end


text:
using Lattices

ups(::Type{ST}, ltc::Lattices.AbstractLattice) where ST = ups(ST, size(ltc))
downs(::Type{ST}, ltc::Lattices.AbstractLattice) where ST = downs(ST, size(ltc))
Random.rand(::Type{ST}, ltc::Lattices.AbstractLattice) where ST = rand(ST, size(ltc))

ups(ltc::Lattices.AbstractLattice) = ups(Bit{Int}, ltc)
downs(ltc::Lattices.AbstractLattice) = downs(Bit{Int}, ltc)
Random.rand(ltc::Lattices.AbstractLattice) = rand(Bit{Int}, ltc)

HilbertSpace{L}(ltc::Lattices.AbstractLattice) where L = HilbertSpace{L}(size(ltc)...)
HilbertSpace(ltc::Lattices.AbstractLattice) = HilbertSpace{Bit{Int}}(ltc)


text:
"""
    SolidFraction computes the solid fraction (= (1-phi) as a function of T

"""
function SolidFraction(T::Array, Phi_o::Array,dt::Float64; P_bar=[], MeltFrac=[], PhaseRat=[],PhaseDiagrams=[])

    Phi_new =   Phi_o.*0.0;
    dPhi_dt =   Phi_o.*0.0;
    if isempty(P_bar)
        # in case we use a parameterized method only (default)
        SolidFraction_Parameterized!(T, Phi_o, Phi_new, dPhi_dt, dt);

    else
        # melt fraction is computed in a phase-wise manner
        #   which can be:
        #               "none"              -   no melting
        #               "parameterized"     -   parameterized mekting model (aka, our default)
        #               "PD"                -   we use a phase diagram to compute the melt fraction, as specified in PhaseDiagramData
        numPhases   =   size(PhaseRat)[end];
        dim         =   length(size(PhaseRat))-1;
        Phi_melt    =   Phi_o.*0.0;
        for iPhase  =   1:numPhases
            Phi_melt_av = Phi_o.*0.0;
            if  MeltFrac[iPhase]=="parameterized"
                Phi_melt_av_temp = Phi_o.*0.0;

                # Use parameterized melting diagram
                SolidFraction_Parameterized!(T, Phi_o, Phi_melt_av_temp, dPhi_dt, dt);
                Phi_melt_av .= 1.0 .- Phi_melt_av_temp;

            elseif MeltFrac[iPhase]=="PD"
                # Interpolate melt fraction from the phase diagram
                interp_meltWt       =   PhaseDiagrams[iPhase].meltWt;

                for i=eachindex(P_bar)
                    meltWt          =   interp_meltWt(T[i],P_bar[i]);

                    Phi_melt_av[i]  =  meltWt;
                end

            elseif MeltFrac[iPhase]=="none"
                # remains zero

            else
                error("Unknown melt fraction type for phase $iPhase: namely $(MeltFrac[iPhase]). Choose from: [none, PD, parameterized] ")
            end

            if dim==2
                Phi_melt +=  Phi_melt_av.*PhaseRat[:,:,iPhase]
            elseif dim==3
                Phi_melt +=  Phi_melt_av.*PhaseRat[:,:,:,iPhase]
            end

        end
        Phi_new   .=   1.0 .- Phi_melt;    # Phi=solid fraction
        dPhi_dt    =   (Phi_new .- Phi_o)./dt;

    end

    return Phi_new, dPhi_dt
end

function SolidFraction_Parameterized!(T::Array, Phi_o::Array, Phi::Array, dPhi_dt::Array, dt::Float64)
   # Compute the melt fraction of the domain, assuming T=Celcius
    # Taken from L.Caricchi (pers. comm.)

    # Also compute dPhi/dt, which is used to compute latent heat

    #Theta      =   (800.0 .- T)./23.0;
    #Phi        =   1.0 .- 1.0./(1.0 .+ exp.(Theta));

    Phi       .=   1.0 .- 1.0./(1.0 .+ exp.((800.0 .- T)./23.0));

    dPhi_dt   .=   (Phi .- Phi_o)./dt;
    Phi_o     .=   Phi;

end

"""
    This computes the lithostatic pressure from a given densit matrix

        General form:

        P = ComputeLithostaticPressure(Rho, Grid);

        with:
            Grid:   2D or 3D arrays that describe the 1D grid coordinates
                    2D - (x,z)
                    3D - (x,y,z)

            Rho:    2D or 3D matrix with density distribution

        out:
            P:      2D or 3D matrix with pressure [in bar!]
"""
function ComputeLithostaticPressure(Rho, Grid);
    # This computes the lithostatic pressure [in bar] from a given density matrix
    g   =   9.81;                   # m/s2
    z   =   Grid[end]               # z coordinates
    dim =   length(size(Rho))

    dz  =   z[2]-z[1];              # note that we assume a constant spacing in z
    P   =   Rho*g*dz;

    if dim==2;
        P[end,:]    .= 0.0;
    elseif dim==3
        P[end,:,:]  .= 0.0;
    end
    P   =   reverse(P,dims=dim);        # reverse array as we go from top->bottom
    P   =   cumsum( P/1e5, dims=dim);   # sum
    P   =   reverse(P,dims=dim);        # same

    return P;
end

struct PhaseDiagram
    meltRho
    meltWt
    rockRho
    rockVp
    rockVs
    rockVpVs
    meltVp
    meltVs
    meltVpVs
end

"""
    This preloads phase diagram data from disk

        General form:

        PhaseDiagramData = LoadPhaseDiagrams(PhaseDiagramNames, PlotDiagrams=false);

        with:
            PhaseDiagramNames:  Array with names that are either [""] or contain the name & directory of the phase diagram
            PlotDiagrams:       Plot the diagrams [warning; to be removed]

        out:
            PhaseDiagramData:   Array with interpolation objects that describe the phase diagrams

"""
function LoadPhaseDiagrams(PhaseDiagramNames, PlotDiagrams=false);
    # This pre-loads phase diagram and creates the interpolation objects to
    #   efficiently query them @ a later stage.
    #
    #   The phase diagram input format is the LaMEM input format, which is
    #   to a large extend similar to what Perple_X gives as an output
    #   (with more comments @ the beginning)

    PhaseData   =   [PhaseDiagram([],[],[],  [], [], [], [], [], []) for i=1:length(PhaseDiagramNames)];    # initialize
    i           =   0;
    for PhaseDiagramName in PhaseDiagramNames
        i   +=   1;

        if length(PhaseDiagramName)>0

            # open diagram and read info about the size of the diagram
            fid             =   open(PhaseDiagramName,"r")
            line            =   readline(fid)
            numFields       =   parse(Int64,line)                       # the # of fields encoded in the diagram
            for i=1:48; line=readline(fid);  end                        # skip to line 50, where the first "real" input starts
            Tmin            =   parse(Float64,readline(fid))-273.15     # minimum T in C
            dT              =   parse(Float64,readline(fid))            # step-size in T [C]
            numT            =   parse(Int64,readline(fid))              # number of T points
            Pmin            =   parse(Float64,readline(fid))            # minimum P in bar
            dP              =   parse(Float64,readline(fid))            # step-size in P [bar]
            numP            =   parse(Int64,readline(fid))              # number of P points
            close(fid)

            # Read the full data from the file
            data_diagram    =   CSV.File(PhaseDiagramName, skipto=56, header=false);    # read all diagram data

            Pvec            =   Pmin:dP:dP*(numP-1)+Pmin;
            Tvec            =   Tmin:dT:dT*(numT-1)+Tmin;
            meltRho         =   reshape(data_diagram.Column1, (numT,numP));              # density of melt
            meltWt          =   reshape(data_diagram.Column2, (numT,numP));              # melt fraction
            rockRho         =   reshape(data_diagram.Column3, (numT,numP));              # density of solid rock

            intp_meltRho    =   LinearInterpolation((Tvec, Pvec), meltRho,  extrapolation_bc = Line());
            intp_meltWt     =   LinearInterpolation((Tvec, Pvec), meltWt,   extrapolation_bc = Line());
            intp_rockRho    =   LinearInterpolation((Tvec, Pvec), rockRho,  extrapolation_bc = Line());
            PhaseData[i]    =   PhaseDiagram(intp_meltRho,intp_meltWt,intp_rockRho, [], [], [], [], [], []);

            if numFields>5
                # Seismic velocities are encoded as well
                rockVp          =   reshape(data_diagram.Column6, (numT,numP));              # Vp velocity of rock [km/s]
                rockVs          =   reshape(data_diagram.Column7, (numT,numP));              # Vs velocity of rock [km/s]
                rockVpVs        =   reshape(data_diagram.Column8, (numT,numP));              # Vp/Vs velocity of rock [-]

                meltVp          =   reshape(data_diagram.Column9, (numT,numP));              # Vp velocity of melt [km/s]
                meltVs          =   reshape(data_diagram.Column10,(numT,numP));              # Vs velocity of melt [km/s]
                meltVpVs        =   reshape(data_diagram.Column11,(numT,numP));              # Vp/Vs velocity of melt [-]

                # add interpolation objects to PhaseDiagramData struct
                intp_rockVp     =   LinearInterpolation((Tvec, Pvec), rockVp,       extrapolation_bc = Line());
                intp_rockVs     =   LinearInterpolation((Tvec, Pvec), rockVs,       extrapolation_bc = Line());
                intp_rockVpVs   =   LinearInterpolation((Tvec, Pvec), rockRho,      extrapolation_bc = Line());

                intp_meltVp     =   LinearInterpolation((Tvec, Pvec), meltVp,       extrapolation_bc = Line());
                intp_meltVs     =   LinearInterpolation((Tvec, Pvec), meltVs,       extrapolation_bc = Line());
                intp_meltVpVs   =   LinearInterpolation((Tvec, Pvec), meltVpVs,     extrapolation_bc = Line());

                # store data in a single structure
                PhaseData[i]    =   PhaseDiagram(intp_meltRho,intp_meltWt,intp_rockRho, intp_rockVp, intp_rockVs, intp_rockVpVs, intp_meltVp, intp_meltVs, intp_meltVpVs);

            elseif numFields==5
                # nothing to be done

            elseif numFields<5
                error("The phase diagram should at least have the collumns: meltRho | meltWt | rockRho | T [K] | P [bar] ")

            end

            # plot phase diagram if requested
            if PlotDiagrams # note that this requires the Plots package to be installed - move this to a separate routine?
                Rho = meltWt.*meltRho .+ (1.0 .- meltWt).*rockRho

                # density
                p2 = heatmap(Tvec,Pvec/1e3,meltWt',  title=PhaseDiagramName,  c=:batlow,  xlabel="T [C]",   ylabel="P [kbar]", dpi=400, fontsize=6, colorbar_title="Melt fraction")
                p3 = heatmap(Tvec,Pvec/1e3,meltRho',                          c=:lajolla, xlabel="T [C]",   ylabel="P [kbar]", dpi=400, fontsize=6, colorbar_title="ρ melt [kg/m3]")
                p1 = heatmap(Tvec,Pvec/1e3,rockRho',                          c=:lajolla, xlabel="T [C]",   ylabel="P [kbar]", dpi=400, fontsize=6, colorbar_title="ρ rock [kg/m3]")
                p4 = heatmap(Tvec,Pvec/1e3,Rho',                              c=:lajolla, xlabel="T [C]",   ylabel="P [kbar]", dpi=400, fontsize=6, colorbar_title="ρ combined[kg/m3]")
                plot(p1,p2,p3,p4); png("PhaseDiagram_DensitiesMeltfraction")

                if numFields>5
                    # Seismic velocities for rock
                    p1 = heatmap(Tvec,Pvec/1e3,rockVp',                         c=:batlow,  xlabel="T [C]",   ylabel="P [kbar]", dpi=400, fontsize=6, colorbar_title="rock Vp [km/s]")
                    p2 = heatmap(Tvec,Pvec/1e3,rockVs',  title=PhaseDiagramName,c=:batlow,  xlabel="T [C]",   ylabel="P [kbar]", dpi=400, fontsize=6, colorbar_title="rock Vs [km/s]")
                    p3 = heatmap(Tvec,Pvec/1e3,rockVpVs',                       c=:batlow,  xlabel="T [C]",   ylabel="P [kbar]", dpi=400, fontsize=6, colorbar_title="rock Vp/Vs []")
                    plot(p1,p2,p3); png("PhaseDiagram_RockSeismicVelocities")

                    # Seismic velocities for melt
                    p1 = heatmap(Tvec,Pvec/1e3,meltVp',                         c=:batlow,  xlabel="T [C]",   ylabel="P [kbar]", dpi=400, fontsize=6, colorbar_title="melt Vp [km/s]")
                    p2 = heatmap(Tvec,Pvec/1e3,meltVs',  title=PhaseDiagramName,c=:batlow,  xlabel="T [C]",   ylabel="P [kbar]", dpi=400, fontsize=6, colorbar_title="melt Vs [km/s]")
                    p3 = heatmap(Tvec,Pvec/1e3,meltVpVs',                       c=:batlow,  xlabel="T [C]",   ylabel="P [kbar]", dpi=400, fontsize=6, colorbar_title="melt Vp/Vs []")
                    plot(p1,p2,p3); png("PhaseDiagram_MeltSeismicVelocities")

                    # Combined velocities
                    Vp      =  meltWt.*meltVp   .+ (1.0 .- meltWt).*rockVp;
                    Vs      =  meltWt.*meltVs   .+ (1.0 .- meltWt).*rockVs;
                    VpVs    =  meltWt.*meltVpVs .+ (1.0 .- meltWt).*rockVpVs;
                    p1 = heatmap(Tvec,Pvec/1e3,Vp',                         c=:batlow,  xlabel="T [C]",   ylabel="P [kbar]", dpi=400, fontsize=6, colorbar_title="melt Vp [km/s]")
                    p2 = heatmap(Tvec,Pvec/1e3,Vs',  title=PhaseDiagramName,c=:batlow,  xlabel="T [C]",   ylabel="P [kbar]", dpi=400, fontsize=6, colorbar_title="melt Vs [km/s]")
                    p3 = heatmap(Tvec,Pvec/1e3,VpVs',                       c=:batlow,  xlabel="T [C]",   ylabel="P [kbar]", dpi=400, fontsize=6, colorbar_title="melt Vp/Vs []")
                    plot(p1,p2,p3); png("PhaseDiagram_CombinedSeismicVelocities")

                end

            end

        else
            PhaseData[i]    =   PhaseDiagram([],[],[],  [], [], [], [], [], []);
        end

    end

    println("Preloaded phase diagram data from disk ...")

    return PhaseData # Send back the interpolation objects

end

"""
    This computes the density distribution and pressure in the domain

        General form:

        Rho_new, P_bar = ComputeDensityAndPressure(Rho, T, FullGrid, Grid, Tracers, ρ, PhaseDiagramData);

        with:

        out:

"""
function ComputeDensityAndPressure(Rho, T, FullGrid, Grid, Tracers, ρ, PhaseDiagramData);
    # This computes density and pressure for the domain.
    # Nonlinear iterations are employed, as (lithostatic) pressure depends on density

    dim         =   length(Grid);

    # compute phase ratio @ every point
    PhaseRatio  =   PhaseRatioFromTracers(FullGrid, Grid, Tracers, "DistanceWeighted");
    numPhases   =   size(PhaseRatio)[end];

    if length(PhaseDiagramData) < numPhases
        error("The array PhaseDiagrams has less entries than the maximum detected phase")
    end

    # Compute pressure
    P_bar_new   =   ComputeLithostaticPressure(Rho, Grid);  #
    Error       =   1.0;
    P_bar       =   P_bar_new;
    Rho_new     =   P_bar .* 0.0;
    it          =   1;
    while (Error>1e-3)  & (it<20)           # density changes P, which changes density (if using phase diagrams), which is why we use iterations
        it          +=  it;
        P_bar       =   P_bar_new;

        Rho_new     =   P_bar .* 0.0;
        Rho_Phase   =   ones(size(Rho_new));
        for iPhase  =   1:numPhases

            if isa(ρ[iPhase], Number)
                # We have a constant density
                Rho_Phase = Rho_Phase*0.0 .+ ρ[iPhase];

            elseif ρ[iPhase]=="PD"
                # we interpolate density from the phase diagram
                interp_meltRho      =   PhaseDiagramData[iPhase].meltRho;
                interp_rockRho      =   PhaseDiagramData[iPhase].rockRho;
                interp_meltWt       =   PhaseDiagramData[iPhase].meltWt;

                for i=eachindex(P_bar)
                    rockRho         =   interp_rockRho(T[i],P_bar[i]);
                    meltRho         =   interp_meltRho(T[i],P_bar[i]);
                    meltWt          =   interp_meltWt(T[i],P_bar[i]);

                    Rho_Phase[i]    =   (1.0 - meltWt)*rockRho + meltRho*meltWt;
                end
            end

            if dim==2
                Rho_new +=  Rho_Phase.*PhaseRatio[:,:,iPhase]
            elseif dim==3
                Rho_new +=  Rho_Phase.*PhaseRatio[:,:,:,iPhase]
            end

        end

        P_bar_new     =   ComputeLithostaticPressure(Rho_new, Grid);  #

        Error = maximum( abs.(P_bar_new-P_bar) )/maximum(abs.(P_bar_new))
        #@show Error
    end

    return Rho_new, P_bar, PhaseRatio

end

"""
    This computes the average of a certain property

        General form:

            PhaseRatioAverage!(Average, prop_vec,  PhaseRatio)

        with:
            PhaseRatio - PhaseRatio matrix for the current grid

"""
function PhaseRatioAverage!(Average::Array, prop_vec,  PhaseRatio)

    numPhases   =   size(PhaseRatio)[end];
    dim         =   length(size(PhaseRatio))-1;

    # Catch errors
    if numPhases>length(prop_vec)
        error("you did not define sufficient properties.")
    end

    Average     =   Average.*0.0;
    for iPhase  =   1:numPhases

        if      dim==2
            Average +=  prop_vec[iPhase].*PhaseRatio[:,:,iPhase]
        elseif  dim==3
            Average +=  prop_vec[iPhase].*PhaseRatio[:,:,:,iPhase]
        end

    end

    return Average
end

"""
    This computes seismic velocities

        General form:

            Vp,Vs,VpVs = ComputeSeismicVelocities(Grid, T, P_bar, PhaseRatio, PhaseDiagramData)

        with:

"""
function ComputeSeismicVelocities(Grid, T, P_bar, PhaseRatio, PhaseDiagramData)
    # This computes seismic velocities on the given grid, taking melt fraction into account
    # This requires phase diagrams that list Vp,Vs etc. as a function of P and T

    dim         =   length(Grid);
    Vp          =   zeros(size(T));
    Vs          =   zeros(size(T));
    VpVs        =   zeros(size(T));

    numPhases   =   size(PhaseRatio)[end];
    for iPhase  =   1:numPhases

            Vp_av               =   zeros(size(T));
            Vs_av               =   zeros(size(T));
            VpVs_av             =   zeros(size(T));

            interp_meltVp       =   PhaseDiagramData[iPhase].meltVp;
            interp_meltVs       =   PhaseDiagramData[iPhase].meltVs;
            interp_meltVpVs     =   PhaseDiagramData[iPhase].meltVpVs;

            interp_rockVp       =   PhaseDiagramData[iPhase].rockVp;
            interp_rockVs       =   PhaseDiagramData[iPhase].rockVs;
            interp_rockVpVs     =   PhaseDiagramData[iPhase].rockVpVs;

            interp_meltWt       =   PhaseDiagramData[iPhase].meltWt;

            if !isempty(interp_meltVp)
                for i=eachindex(P_bar)
                    Phi_melt        =   interp_meltWt(  T[i],   P_bar[i]);

                    meltVp          =   interp_meltVp(  T[i],   P_bar[i]);
                    rockVp          =   interp_rockVp(  T[i],   P_bar[i]);

                    meltVs          =   interp_meltVs(  T[i],   P_bar[i]);
                    rockVs          =   interp_rockVs(  T[i],   P_bar[i]);

                    rockVpVs        =   interp_rockVpVs(T[i],   P_bar[i]);
                    meltVpVs        =   interp_meltVpVs(T[i],   P_bar[i]);

                    Vp_av[i]        =   rockVp*(1.0 - Phi_melt)     + meltVp*Phi_melt;
                    Vs_av[i]        =   rockVs*(1.0 - Phi_melt)     + meltVs*Phi_melt;
                    VpVs_av[i]      =   rockVpVs*(1.0 - Phi_melt)   + meltVpVs*Phi_melt;

                end

            end

            if dim==2
                Vp      += Vp_av  .*PhaseRatio[:,:,  iPhase];
                Vs      += Vs_av  .*PhaseRatio[:,:,  iPhase];
                VpVs    += VpVs_av.*PhaseRatio[:,:,  iPhase];
            else
                Vp      += Vp_av  .*PhaseRatio[:,:,:,iPhase];
                Vs      += Vs_av  .*PhaseRatio[:,:,:,iPhase];
                VpVs    += VpVs_av.*PhaseRatio[:,:,:,iPhase];
            end
    end

    return Vp,Vs,VpVs

end


text:
module GenPrintf
export printf, @fmt_str

abstract type FormatSpecifier end

struct PrintArgument <: FormatSpecifier end
struct PrintValue{V} <: FormatSpecifier end

struct PrintBoth{A,B} <: FormatSpecifier
  a::A
  b::B
end

const ARGN_ERR =
  "wrong number of arguments sent to printf"

include("impl.jl")
include("meta.jl")
include("cons.jl")
include("list.jl")

"""
Macro to create formatting strings.
{} are where arguments will go.

examples:
  fmt"ab cd"
  fmt"hello, {}!"
"""
macro fmt_str(s)
  values = split(unescape_string(s), "{}")
  print_list(values)
end

"""
Prints formatted output. Use the fmt macro for it
to generate code at compile time.

examples:
  printf(fmt"ab cd")
  printf(fmt"hello, {}!", "world")
  printf("ab{}{}", "c", "d")
"""
@generated function printf(fmt::FormatSpecifier, args...)
  printf_impl(fmt, length(args))
end

function printf(fmt::String, args...)
  i = 1
  values = split(fmt, "{}")
  @assert length(values) == length(args)+1 ARGN_ERR

  for value in values
    print(value)
    if i != length(args) + 1
      print(args[i])
      i += 1
    end
  end
end

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Minuit2"
version = v"6.22.6"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/root-project/root/archive/v6-22-06.tar.gz", "81fe6403a3cf51bb1c411f240d9c233473a833e5738b3abf68ed55d0d27ce1cd")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/root-*/math/minuit2

sed -i '/^add_library(Minuit2$/a SHARED' src/CMakeLists.txt
sed -i '/^add_library(Minuit2Math$/a SHARED' src/math/CMakeLists.txt
mkdir build && cd build
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release -Dminuit2_standalone=ON ..
make -j ${nproc}
make install

if [[ "${target}" == *-mingw* ]]; then
    # We have to manually build the shared library
    mkdir -p "${libdir}"
    cd "${prefix}/lib"
    ar x libMinuit2.dll.a
    cc -shared -o "${libdir}/libMinuit2.dll" *.o
    rm *.o
    ar x libMinuit2Math.dll.a
    cc -shared -o "${libdir}/libMinuit2Math.dll" *.o
    rm *.o
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libMinuit2", :libMinuit2)
    LibraryProduct("libMinuit2Math", :libMinuit2Math)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# This file is a part of ROOTFramework.jl, licensed under the MIT License (MIT).

using Cxx

cxxinclude("TAxis.h")

export TAxisPtr
const TAxisPtr = pcpp"TAxis"

axis_spec(edge::AbstractRange{<:Real}) = begin
    nbins = length(edge) - 1
    (nbins < 1) && error("Number of bins must be at least 1")
    (nbins, Float64(first(edge)), Float64((last(edge))))
end

axis_spec(edges::NTuple{1,AbstractRange{<:Real}}) = (axis_spec(edges[1])...,)
axis_spec(edges::NTuple{2,AbstractRange{<:Real}}) = (axis_spec(edges[1])..., axis_spec(edges[2])...)
axis_spec(edges::NTuple{3,AbstractRange{<:Real}}) = (axis_spec(edges[1])..., axis_spec(edges[2])..., axis_spec(edges[3])...)

axis_AbstractRange(nbins, low, up) = low:((up - low) / nbins):up

Base.AbstractRange{<:Real}(axis::TAxisPtr) = axis_AbstractRange((@cxx axis->GetNbins()), (@cxx axis->GetXmin()), (@cxx axis->GetXmax()))


text:
# export PARAMETER_DEFAULT, PARAMETER_UNSTABLE_BUT_FAST, PARAMETER_STABLE_BUT_SLOW

using MathOptInterface
MOI = MathOptInterface
const MOIU = MOI.Utilities
const AFFEQ{T} = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}

abstract type AbstractBlockMatrix{T} <: AbstractMatrix{T} end
function nblocks end
function block end

function Base.size(bm::AbstractBlockMatrix)
    n = mapreduce(blk -> LinearAlgebra.checksquare(block(bm, blk)),
                  +, 1:nblocks(bm), init=0)
    return (n, n)
end
function Base.getindex(bm::AbstractBlockMatrix, i::Integer, j::Integer)
    (i < 0 || j < 0) && throw(BoundsError(i, j))
    for k in 1:nblocks(bm)
        blk = block(bm, k)
        n = size(blk, 1)
        if i <= n && j <= n
            return blk[i, j]
        elseif i <= n || j <= n
            return 0
        else
            i -= n
            j -= n
        end
    end
    i, j = (i, j) .+ size(bm)
    throw(BoundsError(i, j))
end
Base.getindex(A::AbstractBlockMatrix, I::Tuple) = getindex(A, I...)

abstract type BlockSolution{T} <: AbstractBlockMatrix{T} end
struct PrimalSolution{T} <: BlockSolution{T}
    blocks::Vector{Matrix{T}}
end
# getptr(X::PrimalSolution, blk) = getResultYMat(X.problem, blk)
struct VarDualSolution{T} <: BlockSolution{T}
    blocks::Vector{Matrix{T}}
end
# getptr(X::VarDualSolution, blk) = getResultXMat(X.problem, blk)
nblocks(X::BlockSolution) = length(X.blocks)
function block(X::BlockSolution, blk::Integer)
    return X.blocks[blk]
end
# Needed by MPB_wrapper
function Base.getindex(A::BlockSolution, i::Integer)
    block(A, i)
end

mutable struct Optimizer{T} <: MOI.AbstractOptimizer
    objconstant::T
    objsign::Int
    blockdims::Vector{Int}
    varmap::Vector{Tuple{Int, Int, Int}} # Variable Index vi -> blk, i, j
    b::Vector{T}
    solve_time::Float64
    silent::Bool
    options::Dict{Symbol, Any}
    y::Vector{T}
    X::PrimalSolution{T}
    Z::VarDualSolution{T}
    primalobj::T
    dualobj::T
    phasevalue::Symbol
    tempfile::String
    elemdata::Vector{Any}
	presolve::Bool
    function Optimizer{T}(; presolve::Bool = true) where T
		optimizer = new(
            zero(T), 1, Int[], Tuple{Int, Int, Int}[], T[],
            NaN, false, Dict{Symbol, Any}(), T[], PrimalSolution{T}(Matrix{T}[]), VarDualSolution{T}(Matrix{T}[]), zero(T), zero(T), :noINFO, mktempdir(), [], true)
		if !presolve
			optimizer.presolve = false
		end
		return optimizer
    end
end

varmap(optimizer::Optimizer, vi::MOI.VariableIndex) = optimizer.varmap[vi.value]

function MOI.supports(optimizer::Optimizer, param::MOI.RawParameter)
	return param.name in keys(SET_PARAM)
end
function MOI.set(optimizer::Optimizer, param::MOI.RawParameter, value)
	if !MOI.supports(optimizer, param)
		throw(MOI.UnsupportedAttribute(param))
	end
	optimizer.options[param.name] = value
end
function MOI.get(optimizer::Optimizer, param::MOI.RawParameter)
	# TODO: This gives a poor error message if the name of the parameter is invalid.
	return optimizer.options[param.name]
end

MOI.supports(::Optimizer, ::MOI.Silent) = true
function MOI.set(optimizer::Optimizer, ::MOI.Silent, value::Bool)
	optimizer.silent = value
end
MOI.get(optimizer::Optimizer, ::MOI.Silent) = optimizer.silent

MOI.get(::Optimizer, ::MOI.SolverName) = "SDPA-GMP"

# See https://www.researchgate.net/publication/247456489_SDPA_SemiDefinite_Programming_Algorithm_User's_Manual_-_Version_600
# "SDPA (SemiDefinite Programming Algorithm) User's Manual — Version 6.00" Section 6.2
const RAW_STATUS = Dict(
    :noINFO        => "The iteration has exceeded the maxIteration and stopped with no informationon the primal feasibility and the dual feasibility.",
    :pdOPT => "The normal termination yielding both primal and dual approximate optimal solutions.",
    :pFEAS => "The primal problem got feasible but the iteration has exceeded the maxIteration and stopped.",
    :dFEAS => "The dual problem got feasible but the iteration has exceeded the maxIteration and stopped.",
    :pdFEAS => "Both primal problem and the dual problem got feasible, but the iterationhas exceeded the maxIteration and stopped.",
    :pdINF => "At least one of the primal problem and the dual problem is expected to be infeasible.",
    :pFEAS_dINF => "The primal problem has become feasible but the dual problem is expected to be infeasible.",
    :pINF_dFEAS => "The dual problem has become feasible but the primal problem is expected to be infeasible.",
    :pUNBD => "The primal problem is expected to be unbounded.",
    :dUNBD => "The dual problem is expected to be unbounded.")

function MOI.get(optimizer::Optimizer, ::MOI.RawStatusString)
	return RAW_STATUS[optimizer.phasevalue]
end
function MOI.get(optimizer::Optimizer, ::MOI.SolveTime)
	return optimizer.solve_time
end

function MOI.is_empty(optimizer::Optimizer)
    return iszero(optimizer.objconstant) &&
        optimizer.objsign == 1 &&
        isempty(optimizer.blockdims) &&
        isempty(optimizer.varmap) &&
        isempty(optimizer.b) &&
        optimizer.elemdata == []
end
function MOI.empty!(optimizer::Optimizer{T}) where T
    optimizer.objconstant = zero(Cdouble)
    optimizer.objsign = 1
    empty!(optimizer.blockdims)
    empty!(optimizer.varmap)
    empty!(optimizer.b)
    optimizer.X = PrimalSolution{T}(Matrix{T}[])
    optimizer.Z = VarDualSolution{T}(Matrix{T}[])
    optimizer.y = T[]
    optimizer.phasevalue = :noINFO
    optimizer.tempfile = mktempdir()
    optimizer.elemdata = []
    optimizer.primalobj = zero(T)
    optimizer.dualobj = zero(T)
end

function MOI.supports(
    optimizer::Optimizer,
    ::Union{MOI.ObjectiveSense,
            MOI.ObjectiveFunction{<:Union{MOI.SingleVariable,
                                          MOI.ScalarAffineFunction{T}}}}) where T
    return true
end

function MOI.supports_constraint(
    ::Optimizer, ::Type{MOI.VectorOfVariables}, ::Type{MOI.Reals})
    return false
end
const SupportedSets = Union{MOI.Nonnegatives, MOI.PositiveSemidefiniteConeTriangle}
function MOI.supports_constraint(
    ::Optimizer, ::Type{MOI.VectorOfVariables},
    ::Type{<:SupportedSets})
    return true
end
function MOI.supports_constraint(
    ::Optimizer, ::Type{MOI.ScalarAffineFunction{T}},
    ::Type{MOI.EqualTo{T}}) where T
    return true
end

function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike; kws...)
    return MOIU.automatic_copy_to(dest, src; kws...)
end
MOIU.supports_allocate_load(::Optimizer, copy_names::Bool) = !copy_names

function MOIU.allocate(optimizer::Optimizer, ::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    # To be sure that it is done before load(optimizer, ::ObjectiveFunction, ...), we do it in allocate
    optimizer.objsign = sense == MOI.MIN_SENSE ? -1 : 1
end
function MOIU.allocate(::Optimizer, ::MOI.ObjectiveFunction, ::Union{MOI.SingleVariable, MOI.ScalarAffineFunction}) end

function MOIU.load(::Optimizer, ::MOI.ObjectiveSense, ::MOI.OptimizationSense) end
# Loads objective coefficient α * vi
function load_objective_term!(optimizer::Optimizer{T}, α, vi::MOI.VariableIndex) where {T}
    blk, i, j = varmap(optimizer, vi)
    coef = optimizer.objsign * α
    if i != j
        coef /= 2
    end
    # in SDP format, it is max and in MPB Conic format it is min
    inputElement(optimizer, 0, blk, i, j, convert(T, coef))
end
function MOIU.load(optimizer::Optimizer, ::MOI.ObjectiveFunction, f::MOI.ScalarAffineFunction)
    obj = MOIU.canonical(f)
    optimizer.objconstant = f.constant
    for t in obj.terms
        if !iszero(t.coefficient)
            load_objective_term!(optimizer, t.coefficient, t.variable_index)
        end
    end
end
function MOIU.load(optimizer::Optimizer{T}, ::MOI.ObjectiveFunction, f::MOI.SingleVariable) where T
    load_objective_term!(optimizer, one(T), f.variable)
end

function new_block(optimizer::Optimizer, set::MOI.Nonnegatives)
    push!(optimizer.blockdims, -MOI.dimension(set))
    blk = length(optimizer.blockdims)
    for i in 1:MOI.dimension(set)
        push!(optimizer.varmap, (blk, i, i))
    end
end

function new_block(optimizer::Optimizer, set::MOI.PositiveSemidefiniteConeTriangle)
    push!(optimizer.blockdims, set.side_dimension)
    blk = length(optimizer.blockdims)
    for i in 1:set.side_dimension
        for j in 1:i
            push!(optimizer.varmap, (blk, i, j))
        end
    end
end

function MOIU.allocate_constrained_variables(optimizer::Optimizer,
                                             set::SupportedSets)
    offset = length(optimizer.varmap)
    new_block(optimizer, set)
    ci = MOI.ConstraintIndex{MOI.VectorOfVariables, typeof(set)}(offset + 1)
    return [MOI.VariableIndex(i) for i in offset .+ (1:MOI.dimension(set))], ci
end

function MOIU.load_constrained_variables(
    optimizer::Optimizer, vis::Vector{MOI.VariableIndex},
    ci::MOI.ConstraintIndex{MOI.VectorOfVariables},
    set::SupportedSets)
end

function MOIU.allocate_variables(model::Optimizer, nvars)
end

function MOIU.load_variables(optimizer::Optimizer{T}, nvars) where T
    @assert nvars == length(optimizer.varmap)
    dummy = isempty(optimizer.b)
    if dummy
        optimizer.b = [one(T)]
        optimizer.blockdims = [optimizer.blockdims; -1]
    end
    if dummy
        inputElement(optimizer, 1, length(optimizer.blockdims), 1, 1, one(T))
    end
end

function MOIU.allocate_constraint(optimizer::Optimizer{T},
                                  func::MOI.ScalarAffineFunction,
                                  set::MOI.EqualTo) where T
    push!(optimizer.b, MOI.constant(set))
    return AFFEQ{T}(length(optimizer.b))
end

function MOIU.load_constraint(m::Optimizer{T}, ci::AFFEQ,
                              f::MOI.ScalarAffineFunction, s::MOI.EqualTo) where T
    if !iszero(MOI.constant(f))
        throw(MOI.ScalarFunctionConstantNotZero{
            T, MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}(
                MOI.constant(f)))
    end
    f = MOIU.canonical(f) # sum terms with same variables and same outputindex
    for t in f.terms
        if !iszero(t.coefficient)
            blk, i, j = varmap(m, t.variable_index)
            coef = t.coefficient
            if i != j
                coef /= 2
            end
            inputElement(m, ci.value, blk, i, j, convert(T, coef))
        end
    end
end

function MOI.optimize!(m::Optimizer)
	start_time = time()
    # SDPA.initializeUpperTriangle(m.problem, false)
    redundant_F = initializeSolve(m)
    # SDPA.solve(m)
    inputname = "input.dat-s"
    outputname = "output.dat"
    full_input_path = joinpath(m.tempfile, inputname)
    full_output_path = joinpath(m.tempfile, outputname)
    sdpa_gmp_binary_solve!(m, full_input_path, full_output_path, redundant_entries = redundant_F)
    m.solve_time = time() - start_time
end

function MOI.get(m::Optimizer, ::MOI.TerminationStatus)
    status = m.phasevalue
    if status == :noINFO
        return MOI.OPTIMIZE_NOT_CALLED
    elseif status == :pFEAS
        return MOI.SLOW_PROGRESS
    elseif status == :dFEAS
        return MOI.SLOW_PROGRESS
    elseif status == :pdFEAS
        return MOI.OPTIMAL
    elseif status == :pdINF
        return MOI.INFEASIBLE_OR_UNBOUNDED
    elseif status == :pFEAS_dINF
        return MOI.DUAL_INFEASIBLE
    elseif status == :pINF_dFEAS
        return MOI.INFEASIBLE
    elseif status == :pdOPT
        return MOI.OPTIMAL
    elseif status == :pUNBD
        return MOI.DUAL_INFEASIBLE
    elseif status == :dUNBD
        return MOI.INFEASIBLE
    end
end

function MOI.get(m::Optimizer, ::MOI.PrimalStatus)
    status = m.phasevalue
    if status == :noINFO
        return MOI.UNKNOWN_RESULT_STATUS
    elseif status == :pFEAS
        return MOI.FEASIBLE_POINT
    elseif status == :dFEAS
        return MOI.UNKNOWN_RESULT_STATUS
    elseif status == :pdFEAS
        return MOI.FEASIBLE_POINT
    elseif status == :pdINF
        return MOI.UNKNOWN_RESULT_STATUS
    elseif status == :pFEAS_dINF
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif status == :pINF_dFEAS
        return MOI.INFEASIBLE_POINT
    elseif status == :pdOPT
        return MOI.FEASIBLE_POINT
    elseif status == :pUNBD
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif status == :dUNBD
        return MOI.INFEASIBLE_POINT
    end
end

function MOI.get(m::Optimizer, ::MOI.DualStatus)
    status = m.phasevalue
    if status == :noINFO
        return MOI.UNKNOWN_RESULT_STATUS
    elseif status == :pFEAS
        return MOI.UNKNOWN_RESULT_STATUS
    elseif status == :dFEAS
        return MOI.FEASIBLE_POINT
    elseif status == :pdFEAS
        return MOI.FEASIBLE_POINT
    elseif status == :pdINF
        return MOI.UNKNOWN_RESULT_STATUS
    elseif status == :pFEAS_dINF
        return MOI.INFEASIBLE_POINT
    elseif status == :pINF_dFEAS
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif status == :pdOPT
        return MOI.FEASIBLE_POINT
    elseif status == :pUNBD
        return MOI.INFEASIBLE_POINT
    elseif status == :dUNBD
        return MOI.INFEASIBILITY_CERTIFICATE
    end
end

MOI.get(m::Optimizer, ::MOI.ResultCount) = 1
function MOI.get(m::Optimizer, ::MOI.ObjectiveValue)
    return m.objsign * m.primalobj + m.objconstant
end
function MOI.get(m::Optimizer, ::MOI.DualObjectiveValue)
    return m.objsign * m.dualobj + m.objconstant
end
struct PrimalSolutionMatrix <: MOI.AbstractModelAttribute end
MOI.is_set_by_optimize(::PrimalSolutionMatrix) = true
MOI.get(optimizer::Optimizer, ::PrimalSolutionMatrix) = optimizer.X

struct DualSolutionVector <: MOI.AbstractModelAttribute end
MOI.is_set_by_optimize(::DualSolutionVector) = true
function MOI.get(optimizer::Optimizer, ::DualSolutionVector)
    return optimizer.y
end

struct DualSlackMatrix <: MOI.AbstractModelAttribute end
MOI.is_set_by_optimize(::DualSlackMatrix) = true
MOI.get(optimizer::Optimizer, ::DualSlackMatrix) = optimizer.Z

function block(optimizer::Optimizer, ci::MOI.ConstraintIndex{MOI.VectorOfVariables})
    return optimizer.varmap[ci.value][1]
end
function dimension(optimizer::Optimizer, ci::MOI.ConstraintIndex{MOI.VectorOfVariables})
    blockdim = optimizer.blockdims[block(optimizer, ci)]
    if blockdim < 0
        return -blockdim
    else
        return MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(blockdim))
    end
end
function vectorize_block(M, blk::Integer, s::Type{MOI.Nonnegatives})
    return diag(block(M, blk))
end
function vectorize_block(M::AbstractMatrix{T}, blk::Integer, s::Type{MOI.PositiveSemidefiniteConeTriangle}) where T
    B = block(M, blk)
    d = LinearAlgebra.checksquare(B)
    n = MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(d))
    v = Vector{T}(undef, n)
    k = 0
    for j in 1:d
        for i in 1:j
            k += 1
            v[k] = B[i, j]
        end
    end
    @assert k == n
    return v
end

function MOI.get(optimizer::Optimizer, ::MOI.VariablePrimal, vi::MOI.VariableIndex)
    blk, i, j = varmap(optimizer, vi)
    return block(MOI.get(optimizer, PrimalSolutionMatrix()), blk)[i, j]
end

function MOI.get(optimizer::Optimizer, ::MOI.ConstraintPrimal,
                 ci::MOI.ConstraintIndex{MOI.VectorOfVariables, S}) where S<:SupportedSets
    return vectorize_block(MOI.get(optimizer, PrimalSolutionMatrix()), block(optimizer, ci), S)
end
function MOI.get(m::Optimizer, ::MOI.ConstraintPrimal, ci::AFFEQ)
    return m.b[ci.value]
end

function MOI.get(optimizer::Optimizer, ::MOI.ConstraintDual,
                 ci::MOI.ConstraintIndex{MOI.VectorOfVariables, S}) where S<:SupportedSets
    return vectorize_block(MOI.get(optimizer, DualSlackMatrix()), block(optimizer, ci), S)
end
function MOI.get(optimizer::Optimizer, ::MOI.ConstraintDual, ci::AFFEQ)
    return -MOI.get(optimizer, DualSolutionVector())[ci.value]
end


text:
using Turing

function f()
  t = TArray(Int, 1);
  t[1] = 0;
  while true
    produce(t[1])
    t[1] = 1 + t[1]
  end
end

t = Task(f)

consume(t); consume(t)
a = copy(t);
consume(a); consume(a)

Base.@assert consume(t) == 2
Base.@assert consume(a) == 4

# Base.@assert TArray(Float64,  5)[1] != 0
Base.@assert tzeros(Float64, 5)[1]==0
Base.@assert tzeros(4)[1]==0


text:
using JLD2
using Statistics

mutable struct AutostepStepSize <: LoggerKeyData
    step_sizes::Array{Any}
    log_interval::Int

    function AutostepStepSize(logger_init_info)
        num_logged_steps = fld(logger_init_info[LoggerInitKey.TOTAL_STEPS], logger_init_info[LoggerInitKey.INTERVAL])
        println("The AutostepStepSize logger might be pretty inefficient and should probably only be used for Auto debug. If you are planning on running a sweep, it might be worth not logging step sizes.")
        new([], logger_init_info[LoggerInitKey.INTERVAL])
    end
end

function lg_step!(self::AutostepStepSize, env, agent, s, a, s_next, r, is_terminal, cur_step_in_episode, cur_step_total)
    if rem(cur_step_total, self.log_interval) == 0
        ind = fld(cur_step_total, self.log_interval)
        demon_lu = update(agent.demon_learner)
        opt = demon_lu.opt
        @assert demon_lu.opt isa Auto "Demon Optimizer is not Auto! Please set it to Auto or remove LoggerKey.AUTOSTEP_STEPSIZE from logger keys"

        step_size = []
        for k in keys(opt.α)
            α = opt.α[k]
            push!(step_size, copy(α))
        end
        push!(self.step_sizes, step_size)
    end
end

function lg_episode_end!(self::AutostepStepSize, cur_step_in_episode, cur_step_total)
end

function save_log(self::AutostepStepSize, save_dict::Dict)
    save_dict[:autostep_stepsize] = self.step_sizes
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
"""
# Prologue

- Books and algorithms
- Enter Fibonacci
    - An exponential algorithm
    - A polynomial algorithm
    - More careful analysis
- Big-O notation
"""
module Chapter0

export
    fib1,
    fib2

function fib1(n)
    if n == 0
        return 0
    end
    if n == 1
        return 1
    end
    return fib1(n - 1) + fib1(n - 2)
end

function fib2(n)
    if n == 0
        return 0
    end

    f = fill(NaN, n)
    f[0] = 0
    f[1] = 1
    for i ∈ 2:n
        f[i] = f[i - 1] + f[i - 2]
    end
    return f[n]
end

end # module


text:
module TestFailure
# This is for testing failure cases of a packages tests based on the process exiting abnormally
# See ../test/runtests.jl for details
end # module


text:
# construct the extensive formulation

using JuMP, Gurobi;
using NPZ;
# using NLopt;

function find_topology(ks, fs, gs)

    ### find network topology given degree sequence, steady states, and network dynamics
    # ks: list of degree sequence; fs, gs: lists of function values for self dynamics and interaction between adjacent nodes at steady states
    ###

    # initialize model
    model = Model(optimizer_with_attributes(() -> Gurobi.Optimizer(), "OutputFlag" => 1));

    # define variables
    n_node = length(ks);
    @variable(model, Adj[i in 1:n_node, j in 1:n_node], Bin);
    @variable(model, aux_var[i in 1:n_node]);

    # constraints
    # symmetry
    @constraint(model, symmetry[i in 1:n_node, j in 1:n_node], Adj[i,j] == Adj[j,i]);
    # degree sequence
    @constraint(model, degree[i in 1:n_node], sum(Adj[i,j] for j in 1:n_node) == ks[i]);

    # objective function
    @expression(model, cost[i in 1:n_node], fs[i] + sum(Adj[i,j]*gs[i,j] for j in 1:n_node) );
    # linearize the absolute term in the obj
    @constraint(model, LB1[i in 1:n_node], aux_var[i] >= cost[i]);
    @constraint(model, LB2[i in 1:n_node], aux_var[i] >= -cost[i]);

    @objective(model, Min, sum(aux_var[i] for i in 1:n_node) );

    # @NLobjective(model, Min, sum(cost[i]*cost[i] for i in 1:n_node) );

    # solve the model
    t0 = time();
    optimize!(model);
    t1 = time();

    if termination_status(m) != MOI.OPTIMAL
        @warn("===== Final solution not optimal =====")
        return
    end

    # return solution time and solution
    t_sol = t1 - t0;
    println("===== Solution time for ", n, "-node network is", "t_sol");

    obj = objective_value(model);
    println("===== Optimal objective is ", obj);

    Adj = value.(model[:Adj]);

    return obj, Adj, t_sol;
end

# test
function main()
    # load ks, gs, fs
    cd("C:/code/network_reconstruction/data");
    gs = npzread("gs.npy");
    fs = npzread("fs.npy");
    ks = npzread("ks.npy");
    obj, Adj, t_sol = find_topology(ks, fs, gs);
end

main()


text:
module Quantum

using LinearAlgebra
import Base: ==, length

export Swap, X, Y, Z, H, S, Sdag, T, Tdag, Rx, Ry, Rz, R, CNOT, CZ, CR, CRk
export QuantumRegister

mutable struct QuantumRegister
    qubit_product::Vector{Complex{Float64}}
    QuantumRegister(length::Int) = new(vcat([1], zeros(Complex{Float64}, 2^length - 1)))
end

==(register::QuantumRegister, vector::Vector) = register.qubit_product == vector
==(vector::Vector, register::QuantumRegister) = vector == register.qubit_product
==(a::QuantumRegister, b::QuantumRegister) = a.qubit_product == b.qubit_product

# Because `QuantumRegister` saves the tensor product of all its qubits,
# its length is 2^n (where n is the amount of qubits) and needs to be modified
length(register::QuantumRegister) = log2size(register.qubit_product)

# These functions are defined here because they give an `UndefVarError` bellow

single_qubit_gate(matrix::Matrix) =
    (register::QuantumRegister, at::Int) -> apply!(register, matrix, at)

controlled_gate(matrix::Matrix) = function(register::QuantumRegister, control::Int, target::Int)
    _control = control
    _target = target

    if(control == length(register))
        _control -= 1
        Swap(register, control, _control)
        if(_control == target)
            _target = control
        end
    end

    if(_control + 1 != _target)
        _target = _control + 1
        Swap(register, target, _target)
    end

    apply!(register, controlled_matrix(matrix), _control)

    Swap(register, _target, target)

    if(_control == _target)
        Swap(register, _control, control)
    end

    register
end

# This is where the real gates start

Swap(register::QuantumRegister, from::Int, to::Int) =
    apply!(register, swap_matrix(abs(to - from) + 1), min(from, to))

X = single_qubit_gate([0 1 ; 1 0])

Y = single_qubit_gate([0 -im ; im 0])

Z = single_qubit_gate([1 0 ; 0 -1])

H = single_qubit_gate([1/√2 1/√2 ; 1/√2 -1/√2])

S = single_qubit_gate([1 0 ; 0 im])

Sdag = single_qubit_gate([1 0 ; 0 -im])

T = single_qubit_gate([1 0 ; 0 exp(im * pi / 4)])

Tdag = single_qubit_gate([1 0 ; 0 exp(-im * pi / 4)])

Rx(register::QuantumRegister, at::Int, θ::Float64) =
    single_qubit_gate([cos(θ) -im * sin(θ) ; -im * sin(θ) cos(θ)])(register, at)

Ry(register::QuantumRegister, at::Int, θ::Float64) =
    single_qubit_gate([cos(θ / 2) -im * -sin(θ / 2) ; -im * sin(θ) cos(θ / 2)])(register, at)

Rz(register::QuantumRegister, at::Int, θ::Float64) =
    single_qubit_gate([0 1 ; 0 exp(im * θ)])(register, at)

# R is often used to denote Rz
R = Rz

CNOT = controlled_gate([0 1 ; 1 0])

CZ = controlled_gate([1 0 ; 0 -1])

CR(register::QuantumRegister, control::Int, target::Int, θ::Float64) =
    controlled_gate([0 1 ; 0 exp(im * θ)])(register, control, target)

CRk(register::QuantumRegister, control::Int, target::Int, k::Float64) =
    CR(register, control, target, 2π / 2 ^ k)

# Helper functions for the gates

controlled_matrix(matrix::Matrix) = kron([1 0 ; 0 0], I₂) + kron([0 0 ; 0 1], matrix)

function apply!(register::QuantumRegister, matrix::Matrix, from::Int)
    register.qubit_product = pad_matrix(matrix, length(register), from) * register.qubit_product

    register
end

pad_matrix(matrix::Matrix, tosize::Int, from::Int) =
    # kron together identity matrices with `matrix` in index `from`
    foldl(kron, vcat(fill(I₂, from - 1), [matrix], fill(I₂, tosize - from - log2size(matrix) + 1)))

function swap_matrix(span::Int)
    matrix = I + zeros(2^span, 2^span)

    for i in 1:2^(span - 1)
        if iseven(i)
            swap_rows!(matrix, i, i + 2^(span - 1) - 1)
        end
    end

    matrix
end

function swap_rows!(matrix::Matrix, x::Int, y::Int)
    for i in 1:size(matrix, 2)
        temp = matrix[x, i]
        matrix[x, i] = matrix[y, i]
        matrix[y, i] = temp
    end

    matrix
end

log2size(array::Array) = Int(log2(size(array, 1)))

I₂ = [1 0 ; 0 1]

end # module


text:
if !isdefined(:ExoplanetsSysSim) using ExoplanetsSysSim end

# Code for generating clustered planetary systems (once stable can move to src/planetary_system.jl)

function calc_hill_sphere(a::Float64, mu::Float64)
  a*(mu/3)^(1//3)
end

function calc_mutual_hill_radii{StarT<:StarAbstract}(ps::PlanetarySystem{StarT},pl1::Int64, pl2::Int64)
  mu = (ps.planet[pl1].mass+ps.planet[pl2].mass)/ps.star.mass
  a = 0.5*(ps.orbit[pl1].a+ps.orbit[pl2].a)
  calc_hill_sphere(a,mu)
end

function test_stability_circular(P::Vector{Float64},mass::Vector{Float64},star_mass::Float64, sim_param::SimParam)
   @assert length(P) == length(mass)
   const min_num_mutual_hill_radii = get_real(sim_param,"num_mutual_hill_radii")
   #println("# test_stability_circ: mass= ",mass," star_mass=",star_mass)
   #println("# test_stability_circ: mu= ",mass/star_mass, "; P= ",P)
   #print(  "# (aratio,Delta_H)= ")
   found_instability = false
   order = sortperm(P)
   a2 = semimajor_axis(P[order[1]],star_mass)
   for pl in 1:(length(P)-1)
       a1 = a2   # semimajor_axis(P[order[pl]],star_mass)
       a2 = semimajor_axis(P[order[pl+1]],star_mass)
       a = 0.5*(a1+a2)
       mu = (mass[order[pl]]+mass[order[pl+1]])/star_mass
       mutual_hill_radius = calc_hill_sphere(a,mu)
       #print("(",a2/a1,", ",(a2-a1)/mutual_hill_radius,"), ")
       if ! ( a2-a1  >= min_num_mutual_hill_radii*mutual_hill_radius )
          found_instability = true
          break
       end
   end # loop over neighboring planet pairs within cluster
   #println("# test_stability_circ = ",!found_instability)
   return !found_instability
end

function test_stability(P::Vector{Float64},mass::Vector{Float64},star_mass::Float64, sim_param::SimParam; ecc::Vector{Float64} = zeros(length(P)))
   @assert length(P) == length(mass) == length(ecc)
   const min_num_mutual_hill_radii = get_real(sim_param,"num_mutual_hill_radii")
   found_instability = false
   order = sortperm(P)
   a2 = semimajor_axis(P[order[1]],star_mass)
   for pl in 1:(length(P)-1)
       a1 = a2   # semimajor_axis(P[order[pl]],star_mass)
       a2 = semimajor_axis(P[order[pl+1]],star_mass)
       a = 0.5*(a1+a2)
       mu = (mass[order[pl]]+mass[order[pl+1]])/star_mass
       mutual_hill_radius = calc_hill_sphere(a,mu)
       e1 = ecc[order[pl]]
       e2 = ecc[order[pl+1]]
       if ! ( a2*(1-e2)-a1*(1+e1)  >= min_num_mutual_hill_radii*mutual_hill_radius )
          found_instability = true
          break
       end
   end # loop over neighboring planet pairs within cluster
   return !found_instability
end

function calc_if_near_resonance(P::Vector{Float64})
   @assert issorted(P)   # TODO: OPT: Could remove once know it is safe
   const resonance_width = 0.05   # TODO: FEATURE Make a model parameter?  Would need sim_param
   const resonance_width_factor = 1+resonance_width
   const period_ratios_to_check = [ 2.0, 1.5, 4/3, 5/4 ]
   result = falses(length(P))
   if length(P) >= 2
      for i in 1:(length(P)-1)
          for period_ratio in period_ratios_to_check
              if P[i]*period_ratio <= P[i+1] <= P[i]*period_ratio*resonance_width_factor
                 result[i] = true
                 result[i+1] = true
                 break
              end # if near resonance
          end # period_ratio
      end # planets
   end # at least two planets
   return result
end

function generate_planet_periods_sizes_masses_eccs_in_cluster( star::StarAbstract, sim_param::SimParam; n::Int64 = 1 )  # TODO: IMPORTANT: Make this function work and test before using for science
   @assert n>=1
   const generate_planet_mass_from_radius = get_function(sim_param,"generate_planet_mass_from_radius")
   const generate_sizes = get_function(sim_param,"generate_sizes")
   const generate_e_omega = get_function(sim_param,"generate_e_omega")
   const sigma_ecc::Float64 = haskey(sim_param,"sigma_hk") ? get_real(sim_param,"sigma_hk") : 0.0

   if n==1
      R = generate_sizes(star,sim_param)[1]
      mass = map(r->generate_planet_mass_from_radius(r,sim_param),R)
      P = [1.0] # generate_periods_power_law(star,sim_param)
      eccentricity::Float64,  omega_e::Float64 = generate_e_omega(sigma_ecc)::Tuple{Float64,Float64}
      ecc, omega = [eccentricity], [omega_e]
      return P, R, mass, ecc, omega
   end

   # If reach here, then at least 2 planets in cluster

   mean_R = generate_sizes(star,sim_param)[1]
   const sigma_log_radius_in_cluster = get_real(sim_param,"sigma_log_radius_in_cluster")
   #println("# mean_R = ",mean_R," sigma_log_radius_in_cluster= ",sigma_log_radius_in_cluster)

   const min_radius::Float64 = get_real(sim_param,"min_radius")
   const max_radius::Float64 = get_real(sim_param,"max_radius")
   #Rdist = LogNormal(log(mean_R),sigma_log_radius_in_cluster)
   Rdist = Truncated(LogNormal(log(mean_R),sigma_log_radius_in_cluster), min_radius, max_radius) # if we want clustered planet sizes
   R = rand(Rdist,n)
   #R = generate_sizes(star,sim_param, num_pl=n) # if want non-clustered planet sizes instead

   #println("# Rp = ", R)
   mass = map(r->generate_planet_mass_from_radius(r,sim_param),R)
   #println("# mass = ", mass)

   local ecc, omega
   ecc::Array{Float64,1} = Array{Float64}(n)
   omega::Array{Float64,1} = Array{Float64}(n)
   for i in 1:n
        (ecc[i],  omega[i]) = generate_e_omega(sigma_ecc)::Tuple{Float64,Float64}
   end

   const sigma_logperiod_per_pl_in_cluster = get_real(sim_param,"sigma_logperiod_per_pl_in_cluster")
   const min_period = get_real(sim_param,"min_period")
   const max_period = get_real(sim_param,"max_period")
   log_mean_P = 0.0 # log(generate_periods_power_law(star,sim_param))
   # Note: Currently, drawing all periods within a cluster at once and either keeping or rejecting the whole cluster
   #       Should we instead draw periods one at a time?
   max_period_ratio = max_period/min_period
   #Pdist = LogNormal(log_mean_P,sigma_logperiod_per_pl_in_cluster*n)
   Pdist = Truncated(LogNormal(log_mean_P,sigma_logperiod_per_pl_in_cluster*n), 1/sqrt(max_period_ratio), sqrt(max_period_ratio)) #Truncated unscaled period distribution to ensure that the cluster can fit in the period range [min_period, max_period] after scaling by a period scale
   local P
   found_good_periods = false
   attempts = 0
   max_attempts = 100  # Note: Should this be a parameter?
   while !found_good_periods && attempts<max_attempts
      attempts += 1
      P = rand(Pdist,n)
      if test_stability(P, mass, star.mass, sim_param; ecc=ecc)
        found_good_periods = true
     end
   end # while trying to draw periods

   #println("attempts for cluster: ", attempts)

   if !found_good_periods
      #println("# Warning: Did not find a good set of periods, sizes and masses for one cluster.")
      return fill(NaN,n), R, mass, ecc, omega  # Return NaNs for periods to indicate failed
   end
   return P, R, mass, ecc, omega    # Note can also return earlier if only one planet in cluster or if fail to generate a good set of values
end

function generate_num_planets_in_cluster_poisson(s::Star, sim_param::SimParam)
  const lambda::Float64 = exp(get_real(sim_param,"log_rate_planets_per_cluster"))
  const max_planets_in_cluster::Int64 = get_int(sim_param,"max_planets_in_cluster")
  ExoplanetsSysSim.generate_num_planets_poisson(lambda,max_planets_in_cluster,min_planets=1)
end

function generate_num_clusters_poisson(s::Star, sim_param::SimParam)
  const lambda::Float64 = exp(get_real(sim_param,"log_rate_clusters"))
  const max_clusters_in_sys::Int64 = get_int(sim_param,"max_clusters_in_sys")
  ExoplanetsSysSim.generate_num_planets_poisson(lambda,max_clusters_in_sys)
end

# This version generates clustered planetary systems, adaptation of python code from Matthias He
function generate_planetary_system_clustered(star::StarAbstract, sim_param::SimParam; verbose::Bool = false)  # TODO: Make this function work and test before using for science
  # load functions to use for drawing parameters
  const generate_num_clusters = get_function(sim_param,"generate_num_clusters")
  const generate_num_planets_in_cluster = get_function(sim_param,"generate_num_planets_in_cluster")
  const power_law_P = get_real(sim_param,"power_law_P")
  const min_period = get_real(sim_param,"min_period")
  const max_period = get_real(sim_param,"max_period")

  # To check if we have any giant stars in our catalog
  if star.radius > 5
      log_g = log10((6.67408e-8*star.mass*1.989e33)/((star.radius*6.955e10)^2))
      println("R_star: ", star.radius, "; log(g): ", log_g)
  end

  # Generate a set of periods, planet radii, and planet masses.
  attempt_system = 0
  max_attempts_system = 1 #Note: currently this should not matter; each system is always attempted just once
  local num_pl, Plist, Rlist, masslist, ecclist, omegalist
  valid_system = false
  while !valid_system && attempt_system < max_attempts_system

     # First, generate number of clusters (to attempt) and planets (to attempt) in each cluster
     num_clusters = generate_num_clusters(star,sim_param)::Int64
     #num_clusters = 1 # If we want to test singly-clustered model
     num_pl_in_cluster = map(x->generate_num_planets_in_cluster(star, sim_param)::Int64, 1:num_clusters)
     #num_pl_in_cluster = ones(Int64, num_clusters) ##### If we want a non-clustered model, setting each cluster to have 1 planet is equivalent to having no clusters
     num_pl = sum(num_pl_in_cluster)

     #println("num_clusters: ", num_clusters, " ; num_pl_in_clusters", num_pl_in_cluster)

     if( num_pl==0 )
       return PlanetarySystem(star)
     end
     Plist::Array{Float64,1} = Array{Float64}(num_pl)
     Rlist::Array{Float64,1} = Array{Float64}(num_pl)
     masslist::Array{Float64,1} = Array{Float64}(num_pl)
     ecclist::Array{Float64,1} = Array{Float64}(num_pl)
     omegalist::Array{Float64,1} = Array{Float64}(num_pl)
     @assert num_pl_in_cluster[1] >= 1
     pl_start = 1
     pl_stop = 0
     for c in 1:num_clusters
         pl_stop  += num_pl_in_cluster[c]
         Plist_tmp::Array{Float64,1}, Rlist[pl_start:pl_stop], masslist[pl_start:pl_stop], ecclist[pl_start:pl_stop], omegalist[pl_start:pl_stop] = generate_planet_periods_sizes_masses_eccs_in_cluster(star, sim_param, n = num_pl_in_cluster[c])
         valid_cluster = !any(isnan.(Plist_tmp)) #should this be looking for nans in Plist or Plist_tmp? Was Plist but I think it should be Plist_tmp!
         valid_period_scale = false
         attempt_period_scale = 0
         max_attempts_period_scale = 100
         while !valid_period_scale && attempt_period_scale<max_attempts_period_scale && valid_cluster
             attempt_period_scale += 1

             #period_scale::Array{Float64,1} = ExoplanetsSysSim.generate_periods_power_law(star,sim_param)
             period_scale::Array{Float64,1} = power_law_P!=-1.0 ? ExoplanetsSysSim.draw_power_law(power_law_P, min_period/min(Plist_tmp...), max_period/max(Plist_tmp...), 1) : exp(log(min_period/min(Plist_tmp...)).+rand()*log((max_period/max(Plist_tmp...))/(min_period/min(Plist_tmp...))))
             #Note: this ensures that the minimum and maximum periods will be in the range [min_period, max_period]
             #Warning: not sure about the behaviour when min_period/min(Plist_tmp...) > max_period/max(Plist_tmp...) (i.e. when the cluster cannot fit in the given range)?
             #TODO OPT: could draw period_scale more efficiently by computing the allowed regions in [min_period, max_period] given the previous cluster draws

             Plist[pl_start:pl_stop] = Plist_tmp .* period_scale

             #Note: Checking whether the periods are in the range [min_period, max_period] should be redundant now:
             periods_in_range = true
             for p in Plist_tmp .* period_scale # To check if the periods are still in the range [min_period, max_period] after scaling by period_scale
                if (p < min_period) || (p > max_period)
                    periods_in_range = false
                    println("Pscale = ", period_scale, " , P = ", p, " ; Pmax/Pmin = ", max(Plist_tmp...)/min(Plist_tmp...))
                end
             end

             if test_stability(Plist[1:pl_stop], masslist[1:pl_stop], star.mass, sim_param; ecc=ecclist[1:pl_stop]) && periods_in_range
                valid_period_scale = true
             end
         end  # while !valid_period_scale...

         #if attempt_period_scale > 1
             #println("attempts_period_scale: ", attempt_period_scale)
         #end

         if !valid_period_scale
            Plist[pl_start:pl_stop] = NaN
         end
         pl_start += num_pl_in_cluster[c]
     end # for c in 1:num_clusters
     isnanPlist::Array{Bool,1} = isnan.(Plist::Array{Float64,1})
     if any(isnanPlist)  # If any loop failed to generate valid planets, it should set a NaN in the period list
     keep::Array{Bool,1} = .!(isnanPlist)  # Currently, keeping clusters that could be fit, rather than throwing out entirely and starting from scratch.  Is this a good idea?  Matthias tried the other approach in his python code.
        num_pl = sum(keep)
        Plist = Plist[keep]
        Rlist = Rlist[keep]
        masslist = masslist[keep]
        ecclist = ecclist[keep]
        omegalist = omegalist[keep]
     end
     valid_system = false #this was set to true but I think it should be false?

#println("P: ", Plist)
#println("R: ", Rlist)
#println("M: ", masslist)
#println("ecc: ", ecclist)
#println("omega: ", omegalist)

     #Note: This would be for drawing each cluster separately and then accepting or rejecting the whole lot. By testing for stability before adding each cluster, this last test should be unnecessary.
     if length(Plist) > 0
         if test_stability(Plist, masslist, star.mass, sim_param; ecc=ecclist)
            valid_system = true
         else
            println("Warning: re-attempting system because it fails stability test even though its clusters each pass the test.")
            #Note: this should never happen because we check for stability before adding each cluster, and unstable additions are set to NaN and then discarded
         end
     else
         valid_system = true #this else statement is to allow for systems with no planets to pass
     end

     attempt_system += 1
  end # while !valid_system...

  if attempt_system > 1
    println("attempt_system: ", attempt_system)
  end

  # To print out periods, radii, and masses (for troubleshooting):
  #=
  i_sort = sortperm(Plist)
  Plist_sorted = sort(Plist)
  if length(Plist) > 1
      ratio_list = Plist_sorted[2:end]./Plist_sorted[1:end-1]
      if minimum(ratio_list) < 1.1
         println("P: ", Plist_sorted)
         println(Rlist[i_sort])
         println(masslist[i_sort])
         println(ecclist[i_sort])
         println(omegalist[i_sort])
      end
  end
  =#

  # Now assign orbits with given periods, sizes, and masses.
###const generate_e_omega = get_function(sim_param,"generate_e_omega")
###const sigma_ecc::Float64 = haskey(sim_param,"sigma_hk") ? get_real(sim_param,"sigma_hk") : 0.0
  const sigma_incl = deg2rad(get_real(sim_param,"sigma_incl"))
  const sigma_incl_near_mmr = deg2rad(get_real(sim_param,"sigma_incl_near_mmr"))
  const max_incl_sys = get_real(sim_param,"max_incl_sys")
    pl = Array{Planet}(num_pl)
    orbit = Array{Orbit}(num_pl)
    incl_sys = acos(cos(max_incl_sys*pi/180)*rand()) #acos(rand()) for isotropic distribution of system inclinations; acos(cos(X*pi/180)*rand()) gives angles from X (deg) to 90 (deg)
    idx = sortperm(Plist)       # TODO OPT: Check to see if sorting is significant time sink.  If so, could reduce redundant sortperm
    is_near_resonance = calc_if_near_resonance(Plist[idx])
    for i in 1:num_pl
      # if verbose   println("i=",i," idx=",idx," Plist=",Plist[idx] );     end
      #=
      if haskey(sim_param,"sigma_hk_one") && haskey(sim_param,"sigma_hk_multi")
         sigma_ecc = num_pl == 1 ? get_real(sim_param,"sigma_hk_one") : get_real(sim_param,"sigma_hk_multi")
      end
      =#
###(ecc,  omega) = generate_e_omega(sigma_ecc)::Tuple{Float64,Float64}  # WARNING: Not testing for stability after eccentricites drawn.

      sigma_incl_use = is_near_resonance[i] ? sigma_incl_near_mmr : sigma_incl
      incl_mut = sigma_incl_use*sqrt(randn()^2+randn()^2) # rand(Distributions.Rayleigh(sigma_incl))
      asc_node = 2pi*rand()
      mean_anom = 2pi*rand()
      incl =  incl_mut!=zero(incl_mut) ? acos( cos(incl_sys)*cos(incl_mut) + sin(incl_sys)*sin(incl_mut)*cos(asc_node) ) : incl_sys
      orbit[i] = Orbit(Plist[idx[i]], ecclist[idx[i]], incl, omegalist[idx[i]], asc_node, mean_anom)
      pl[i] = Planet( Rlist[idx[i]],  masslist[idx[i]] )
    end # for i in 1:num_pl

  return PlanetarySystem(star,pl,orbit)
end

## test_generate_planetary_system
function test_generate_planetary_system_clustered() # TODO: Update to test to not reply on generate_kepler_physical_catalog
  sim_param = setup_sim_param_model()
  cat_phys = generate_kepler_physical_catalog(sim_param)
end


text:
####
#### Local minimization with NLopt. Loaded on demand when `NLopt` is.
####

import NLopt

import MultistartOptimization: local_minimization

export NLoptLocalMethod

const NLopt_ret_success = Set([:SUCCESS, :STOPVAL_REACHED, :FTOL_REACHED, :XTOL_REACHED,
                               :MAXEVAL_REACHED, :MAXTIME_REACHED])

Base.@kwdef struct NLoptLocalMethod{S}
    algorithm::NLopt.Algorithm
    xtol_abs::Float64 = 1e-8
    xtol_rel::Float64 = 1e-8
    maxeval::Int = 0
    maxtime::Float64 = 0.0
    "Return values which are considered as “success”."
    ret_success::S = NLopt_ret_success
end

"""
$(SIGNATURES)

A wrapper for algorithms supported by `NLopt`. Used to construct the corresponding
optimization problem. All return values in `ret_success` are considered valid (`ret` is also
kept in the result), all negative return values are considered invalid.

See the NLopt documentation for the options. Defaults are changed slightly.
"""
function NLoptLocalMethod(algorithm::NLopt.Algorithm; options...)
    NLoptLocalMethod(; algorithm = algorithm, options...)
end

"""
$(SIGNATURES)

Solve `minimization_problem` using `local_method`, starting from `x`. Return a
`LocationValue`.
"""
function local_minimization(local_method::NLoptLocalMethod,
                            minimization_problem::MinimizationProblem, x)
    @unpack algorithm, xtol_abs, xtol_rel, maxeval, maxtime, ret_success = local_method
    @unpack objective, lower_bounds, upper_bounds = minimization_problem
    opt = NLopt.Opt(algorithm, length(x))
    opt.lower_bounds = lower_bounds
    opt.upper_bounds = upper_bounds

    # if a method objective(x,grad) exists, use it. otherwise assume objective is not differentiable (which was the previous behavior)
    opt.min_objective = applicable(objective,x,x) ? objective : nlopt_nondifferentiable_wrapper(objective)
    opt.xtol_abs = xtol_abs
    opt.xtol_rel = xtol_rel
    opt.maxeval = maxeval
    opt.maxtime = maxtime
    optf, optx, ret = NLopt.optimize(opt, x)
    ret ∈ ret_success ? (value = optf, location = optx, ret = ret) : nothing
end

function nlopt_nondifferentiable_wrapper(fn)
    function f̃(x,grad)              # wrapper for NLopt
        @argcheck isempty(grad)     # ensure no derivatives are asked for
        return fn(x)
    end
    return f̃
end


text:
using MacroTools
using UUIDs

function format_tag(name ::String, attributes::Dict{String, String}, inner::String = ""; opened = false, closed = false)
    attrs_string = join(
        ["$k=\"$v\"" for (k, v) in attributes],
        " "
    )
    tag = "<$name $attrs_string"
    if closed
        tag *= "/>"
    elseif opened
        tag *= ">"
    else
        tag *= ">$inner</$name>"
    end
end

function interpolate_string(s::String; kwargs...)
    result = s
    for (k, v) in kwargs
        result = replace(result, "{%$(k)%}" =>  v)
    end
    return result
end

function validate_index(name::AbstractString, index::AbstractString, checks::Vector)
    missings = filter(checks) do check
        !occursin(check[2], index)
    end
    if !isempty(missings)
        error(
            string(
                "Missing item", (length(missings)>1 ? "s" : ""), " ",
                join(getindex.(missings, 1), ", "),
                " in ", name
            )
        )
    end
end

macro var_str(s)
    return Symbol(s)
end

function parse_props(s)
    function make_prop(part)
        m = match(r"^(?<id>[A-Za-z]+[\w\-\:\.]*)\.(?<prop>[A-Za-z]+[\w\-\:\.]*)$", strip(part))
        if isnothing(m)
            error("expected <id>.<property>[,<id>.<property>...] in $(part)")
        end

        return (Symbol(m[:id]), Symbol(m[:prop]))
    end

    props_parts = split(s, ",", keepempty = false)

    return map(props_parts) do part
        return make_prop(part)
    end
end

function generate_hash()
    return strip(string(UUIDs.uuid4()), '-')
end

"""
    @callid_str"

Macro for crating Dash CallbackId.
Parse string in form "[{State1[, ...]}] Input1[, ...] => Output1[, ...]"

#Examples
```julia
    id1 = callid"{inputDiv.children} input.value => output1.value, output2.value"
```
"""
macro callid_str(s)
    rex = r"(\{(?<state>.*)\})?(?<input>.*)=>(?<output>.*)"ms
    m = match(rex, s)
    if isnothing(m)
        error("expected {state} input => output")
    end
    input = parse_props(strip(m[:input]))
    output = parse_props(strip(m[:output]))
    state = isnothing(m[:state]) ? Vector{IdProp}() : parse_props(strip(m[:state]))
    return CallbackId(state, input, output)
end


text:
# Source: https://github.com/FluxML/model-zoo/blob/master/vision/mnist/conv.jl
#
# Classifies MNIST digits with a convolutional network.
# Writes out saved model to the file "mnist_conv.bson".
# Demonstrates basic model construction, training, saving,
# conditional early-exit, and learning rate scheduling.
#
# This model, while simple, should hit around 99% test
# accuracy after training for approximately 20 epochs.

using Flux, Flux.Data.MNIST, Statistics
using Flux: onehotbatch, onecold, crossentropy, throttle
using Base.Iterators: repeated, partition
using Printf, BSON

# Load labels and images from Flux.Data.MNIST
@info("Loading data set")
train_labels = MNIST.labels()
train_imgs = MNIST.images()

# Bundle images together with labels and group into minibatchess
function make_minibatch(X, Y, idxs)
    X_batch = Array{Float32}(undef, size(X[1])..., 1, length(idxs))
    for i in 1:length(idxs)
        X_batch[:, :, :, i] = Float32.(X[idxs[i]])
    end
    Y_batch = onehotbatch(Y[idxs], 0:9)
    return (X_batch, Y_batch)
end
batch_size = 128
mb_idxs = partition(1:length(train_imgs), batch_size)
train_set = [make_minibatch(train_imgs, train_labels, i) for i in mb_idxs]

# Prepare test set as one giant minibatch:
test_imgs = MNIST.images(:test)
test_labels = MNIST.labels(:test)
test_set = make_minibatch(test_imgs, test_labels, 1:length(test_imgs))

# Define our model.  We will use a simple convolutional architecture with
# three iterations of Conv -> ReLU -> MaxPool, followed by a final Dense
# layer that feeds into a softmax probability output.
@info("Constructing model...")
model = Chain(
    # First convolution, operating upon a 28x28 image
    Conv((3, 3), 1=>16, pad=(1,1), relu),
    MaxPool((2,2)),

    # Second convolution, operating upon a 14x14 image
    Conv((3, 3), 16=>32, pad=(1,1), relu),
    MaxPool((2,2)),

    # Third convolution, operating upon a 7x7 image
    Conv((3, 3), 32=>32, pad=(1,1), relu),
    MaxPool((2,2)),

    # Reshape 3d tensor into a 2d one, at this point it should be (3, 3, 32, N)
    # which is where we get the 288 in the `Dense` layer below:
    x -> reshape(x, :, size(x, 4)),
    Dense(288, 10),

    # Finally, softmax to get nice probabilities
    softmax,
)

# choose device for testing
# device = gpu
device = cpu

train_set = device.(train_set)
test_set = device.(test_set)
model = device(model)

# Make sure our model is nicely precompiled before starting our training loop
model(train_set[1][1])

# `loss()` calculates the crossentropy loss between our prediction `y_hat`
# (calculated from `model(x)`) and the ground truth `y`.  We augment the data
# a bit, adding gaussian random noise to our image to make it more robust.
function loss(x, y)
    # We augment `x` a little bit here, adding in random noise
    x_aug = x .+ 0.1f0*device(randn(eltype(x), size(x)))

    y_hat = model(x_aug)
    return crossentropy(y_hat, y)
end
accuracy(x, y) = mean(onecold(model(x)) .== onecold(y))

# Train our model with the given training set using the ADAM optimizer and
# printing out performance against the test set as we go.
opt = ADAM(0.001)

@info("Beginning training loop...")
best_acc = 0.0
last_improvement = 0
for epoch_idx in 1:10
    # global best_acc, last_improvement
    # Train for a single epoch
    @show epoch_idx
    @time Flux.train!(loss, params(model), train_set, opt)
end


text:
"""
    Cell(;has_mine=false, flagged=false, revealed=false, mine_count=0, idx=(0,0))

Cell represents the state of a cell in Minesweeper
- `has_mine`: true if cell contains a mine
- `flagged`: true if cell has been flagged as having a mine
- `revealed`: true if the cell has been revealed
- `mine_count`: number of mines in adjecent cells
- `idx`: Cartesian indices of cell
"""
mutable struct Cell
    has_mine::Bool
    flagged::Bool
    revealed::Bool
    mine_count::Int
    idx::CartesianIndex
end

function Cell(;has_mine=false, flagged=false, revealed=false, mine_count=0, idx=(0,0))
    coords = CartesianIndex(idx)
    return Cell(has_mine, flagged, revealed, mine_count, coords)
end

"""
    Game(;dims=(12,12), n_mines=40, mines_flagged=0, mine_detonated=false, trials=0)

Minesweeper game object

- `cells`: an array of cells
- `dims`: a Tuple indicating dimensions of cell array
- `n_mines`: number of mines in the game
- `mines_flagged`: number of mines flagged
- `mine_detonated`: indicates whether a mine has been detonated
- `score`: score for game, which includes hits, misses, false alarms and correct rejections
- `trials`: the number of trials or moves
"""
mutable struct Game{T}
    cells::Array{Cell,2}
    dims::Tuple{Int,Int}
    n_mines::Int
    mines_flagged::Int
    mine_detonated::Bool
    score::T
    trials::Int
end

function Game(;dims=(12,12), n_mines=40, mines_flagged=0,
    mine_detonated=false, trials=0)
    score = (hits=0.0,false_alarms=0.0,misses=0.0,correct_rejections=0.0)
    cells = initilize_cells(dims)
    add_mines!(cells, n_mines)
    mine_count!(cells)
    return Game(cells, dims, n_mines, mines_flagged, mine_detonated,
        score, trials)
end

function initilize_cells(dims)
    return [Cell(idx=(r,c)) for r in 1:dims[1], c in 1:dims[2]]
end

function add_mines!(cells, n_mines)
    mines = sample(cells, n_mines, replace=false)
    for m in mines
        m.has_mine = true
    end
    return nothing
end

"""
    mine_count!(cells)

For each cell, omputes the number of neighboring cells containing a mine.
"""
function mine_count!(cells)
    for cell in cells
        neighbors = get_neighbors(cells, cell.idx)
        cell.mine_count = count(x->x.has_mine, neighbors)
    end
    return nothing
end

get_neighbors(game::Game, idx) = get_neighbors(game.cells, idx.I...)
get_neighbors(game::Game, cell::Cell) = get_neighbors(game.cells, cell.idx.I...)
get_neighbors(game::Game, x, y) = get_neighbors(game.cells, x, y)
get_neighbors(cells, idx) = get_neighbors(cells, idx.I...)

function get_neighbors(cells, xr, yr)
    v = [-1,0,1]
    X = xr .+ v; Y = yr .+ v
    neighbors = Vector{Cell}()
    for y in Y, x in X
        (x == xr) && (y == yr) ? (continue) : nothing
        in_bounds(cells, x, y) ? push!(neighbors, cells[x,y]) : nothing
    end
    return neighbors
end

get_neighbor_indices(game::Game, idx) = get_neighbor_indices(game.cells, idx.I...)
get_neighbor_indices(game::Game, x, y) = get_neighbor_indices(game.cells, x, y)
get_neighbor_indices(cells, idx) = get_neighbor_indices(cells, idx.I...)

function get_neighbor_indices(cells, xr, yr)
    v = [-1,0,1]
    X = xr .+ v; Y = yr .+ v
    neighbors = Vector{CartesianIndex}()
    for y in Y, x in X
        (x == xr) && (y == yr) ? (continue) : nothing
        in_bounds(cells, x, y) ? push!(neighbors, CartesianIndex(x,y)) : nothing
    end
    return neighbors
end

function in_bounds(cells, x, y)
    if (x < 1) || (y < 1)
        return false
    elseif (x > size(cells, 1)) || (y > size(cells, 2))
        return false
    else
        return true
    end
end

reveal_zeros!(game::Game, x, y) = reveal_zeros!(game.cells, CartesianIndex(x, y))
reveal_zeros!(game::Game, idx) = reveal_zeros!(game.cells, idx)

function reveal_zeros!(cells, idx)
    cells[idx].mine_count ≠ 0 ? (return nothing) : nothing
    indices = get_neighbor_indices(cells, idx)
    for i in indices
        c = cells[i]
        if !c.has_mine && !c.revealed
            c.revealed = true
            c.mine_count == 0 ? reveal_zeros!(cells, i) : nothing
        end
    end
    return nothing
end

reveal(game::Game) = reveal(game.cells::Array{Cell,2})

"""
    reveal(cells::Array{Cell,2})

Reveals game state in REPL
"""
function reveal(cells::Array{Cell,2})
    s = ""
    for x in 1:size(cells,1)
        for y in 1:size(cells,2)
            c = cells[x,y]
            if c.has_mine
                s *= string(" ", "💣", " ")
            else
                s *= string(" ", c.mine_count, " ")
            end
        end
        s *= "\r\n"
    end
    println(s)
end

Base.show(io::IO, game::Game) = Base.show(io::IO, game.cells)

function Base.show(io::IO, cells::Array{Cell,2})
    s = ""
    for x in 1:size(cells,1)
        for y in 1:size(cells,2)
            c = cells[x,y]
            if c.flagged
                s *= string(" ", "🚩", " ")
                continue
            end
            if c.revealed
                if c.has_mine
                    s *= string(" ", "💣", " ")
                    continue
                end
                if c.mine_count == 0
                    s *= string(" ", "∘", " ")
                else
                    s *= string(" ", c.mine_count, " ")
                end
            else
                s *= string(" ", "■", " ")
            end
        end
        s *= "\r\n"
    end
    println(s)
end

"""
    game_over(game)

Terminates game if mine is detonated or all cells have been revealed or flagged.
"""
function game_over(game)
    if game.mine_detonated
        return true
    end
    return all(x->x.revealed || x.flagged, game.cells)
end

"""
    compute_score!(game)

Computes hits, false alarms, misses and correct rejections
"""
function compute_score!(game)
    cells = game.cells
    n_mines = game.n_mines
    no_mines = prod(game.dims) - n_mines
    hits = count(x->x.flagged && x.has_mine, cells)
    misses = n_mines - hits
    false_alarms = count(x->x.flagged && !x.has_mine, cells)
    cr = no_mines - false_alarms
    game.score = (hits=hits,false_alarms=false_alarms,misses=misses,correct_rejections=cr)
    return nothing
end

"""
    select_cell!(game, x, y)

Select cell given a game and row and column indices.
"""
select_cell!(game, x, y) = select_cell!(game, CartesianIndex(x, y))

"""
    select_cell!(game, x, y)

Select cell given a game and cell object.
"""
select_cell!(game, cell::Cell) = select_cell!(game, cell.idx)

"""
    select_cell!(game, idx)

Select cell given a game and Cartesian index.
"""
function select_cell!(game, idx)
    cell = game.cells[idx]
    cell.revealed = true
    reveal_zeros!(game, idx)
    game.mine_detonated = cell.has_mine
    return nothing
end

function flag_cell!(cell, game, gui::Nothing)
    cell.flagged = true
    game.mines_flagged += 1
    return nothing
end

function setup()
    println("Select Difficulty")
    waiting = true
    levels = [20,50,70]
    i = 0
    while waiting
        d = input("1: easy 2: medium, 3: difficult")
        i = parse(Int, d)
        if i ∉ 1:3
            println("Please select 1, 2, or 3")
        else
            waiting = false
        end
    end
    return Game(;dims=(15,15), n_mines=levels[i])
end

function play()
    game = setup()
    play(game)
end

"""
    play(game)

Play Minesweeper in REPL

Run play() to use with default easy, medium and difficult settings
Pass game object to play custom game
Enter exit to quit
Enter row_index column_index to select cell (e.g. 1 2)
"""
function play(game)
    playing = true
    coords = fill(0, 2)
    println(game)
    while playing
        str = input("select coordinates")
        str == "exit" ? (return) : nothing
        try
            coords = parse.(Int, split(str, " "))
        catch
            error("Please give x and y coordinates e.g. 3 2")
        end
        select_cell!(game, coords...)
        if game.mine_detonated
            reveal(game)
            println("Game Over")
            println(" ")
            d = input("Press y to play again")
            d == "y" ? (game = setup()) : (playing = false)
        else
            println(game)
        end
    end
end

function input(prompt)
    println(prompt)
    return readline()
end

"""
    update!(game, gui::Nothing)

`update!` is used to display game state during simulations.
`update!` can be used with REPL or Gtk GUI if imported via import_gui()
"""
function update!(game, gui::Nothing)
    println(game)
end

"""
imports Gtk related code for the GUI
"""
function import_gui()
    path = @__DIR__
    include(path*"/Gtk_Gui.jl")
end


text:
"""
    Build a map of node IDs to elements
"""
function buildNodeToElemMap!(mesh::Mesh)
    startLog(main_perf_log, "buildNodeToElemMap!()")

    node_to_elem_map = mesh.node_to_elem_map

    for elem in mesh.elements
        for node in elem.nodes
            if !(node.id in keys(node_to_elem_map))
                node_to_elem_map[node.id] = Array{Element}(0)
            end

            push!(node_to_elem_map[node.id], elem)
        end
    end

    stopLog(main_perf_log, "buildNodeToElemMap!()")
end

"""
    Build lists of local elements and nodes
"""
function buildNodeAndElemLists!(mesh::Mesh)
    startLog(main_perf_log, "buildNodeAndElemLists!()")

    proc_id = MPI.Comm_rank(MPI.COMM_WORLD)

    for elem in mesh.elements
        if elem.processor_id == proc_id
            push!(mesh.local_elements, elem)
        end
    end

    for node in mesh.nodes
        if node.processor_id == proc_id
            push!(mesh.local_nodes, node)
        end
    end

    stopLog(main_perf_log, "buildNodeAndElemLists!()")
end

"""
    Called after adding all of the nodes and elements to the mesh.
    Partitions the mesh and assigns processor IDs
"""
function initialize!(mesh::Mesh)
    startLog(main_perf_log, "initialize!()")

    buildNodeToElemMap!(mesh)

    if MPI.Comm_size(MPI.COMM_WORLD) <= 2
        partition!(mesh, SimplePartitioner)
    else
        partition!(mesh, MetisPartitioner)
    end

    buildNodeAndElemLists!(mesh)

    stopLog(main_perf_log, "initialize!()")
end


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

using StatisticalRethinking
using StatsPlots

ProjDir = @__DIR__

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

if success(rc)

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

end