text:
import StaticArrays: SVector
import EffectiveWaves: AbstractAzimuthalSymmetry

using HCubature, Interpolations, Statistics

function material_scattering_coefficients_discrete(discrete_coefs, ω::T, source::AbstractSource{T,Acoustic{T,Dim}}, material::Material{Dim,Sphere{T,Dim}};
        basis_order::Int = 1,
        basis_field_order::Int = 2,
        rtol::T = 1e-2,
        maxevals::Int = Int(5e4)
    ) where {T,Dim}

    lm2n = lm_to_spherical_harmonic_index
    rθφ2xyz = radial_to_cartesian_coordinates

    v = regular_basis_function(source.medium,  ω)

    numdensity = number_density(material.species)

    function kernel(rθφ)
        x = rθφ2xyz(rθφ)

        vs = conj.(v(basis_order + basis_field_order, x))
        fs = discrete_coefs(x)

        (rθφ[1]^2 * sin(rθφ[2]) * numdensity) .* [
            sum(
                gaunt_coefficient(l,m,dl,dm,l1,m-dm) * vs[lm2n(l1,m-dm)] * fs[lm2n(dl,dm)]
            for dl = 0:basis_order for dm = -dl:dl for l1 in max(abs(m-dm),abs(dl-l)):(dl+l))
        for l = 0:basis_field_order for m = -l:l]
    end

    (v,err) = hcubature(kernel, SVector(0.0,0.0,-π), SVector(R-outer_radius(s1),π,π)
        ; rtol=rtol, maxevals = maxevals
    )

    println("Integration error of $err")

    return v

end

function discrete_system_residue(discrete_coefs, ω::T, source::AbstractSource{T,Acoustic{T,Dim}}, material::Material{Dim,Sphere{T,Dim}}, ::AbstractAzimuthalSymmetry{Dim};
        basis_order::Int = 1,
        mesh_points::Int = 11,
        rtol::T = 1e-2,
        maxevals::Int = Int(1e5),
        pair_corr = hole_correction_pair_correlation
    ) where {T,Dim}

    if length(material.species) > 1
        @warn "discrete_system has only been implemented for 1 species for now. Will use only first specie."
    end

    s1 = material.species[1]
    scale_number_density = one(T) - one(T) / material.numberofparticles
    bar_numdensity = scale_number_density * number_density(s1)

    R = outer_radius(material.shape)

    gs = regular_spherical_coefficients(source)(basis_order,origin(material.shape),ω);
    v = regular_basis_function(source.medium,  ω)

    Uout = outgoing_translation_matrix(ω, source.medium, material;
        basis_order = basis_order, tol = rtol
    )

    t_matrices = get_t_matrices(source.medium, material.species, ω, basis_order)
    t_diags = diag.(t_matrices)

    rθφ2xyz = radial_to_cartesian_coordinates

    len = basisorder_to_basislength(Acoustic{T,Dim}, basis_order)

    function incident_coefficients(rθ1::AbstractVector{T})
        lm2n = lm_to_spherical_harmonic_index
        vs = v(2basis_order, rθφ2xyz([rθ1; 0.0]))

        coefs = [
            t_diags[1][lm2n(l,m)] *
            sum(
                gaunt_coefficient(dl,0,l,m,l1,-m) * vs[lm2n(l1,-m)] * gs[lm2n(dl,0)]
            for l1 in abs(m):(2basis_order), dl in 0:basis_order)
        for l = 0:basis_order for m = -l:l]

        return coefs
    end

    ls, ms = spherical_harmonics_indices(basis_order)

    function kernel_function(rθ1::AbstractVector{T})
        x1 = rθφ2xyz(SVector(rθ1[1],rθ1[2],zero(T)))

        fun = function (rθφ::SVector{3,T})
            x2 = rθφ2xyz(rθφ)
            if pair_corr(x1,s1,x2,s1) ≈ zero(T)
                return zeros(Complex{T}, len)
            end
            # U = outgoing_translation_matrix(medium, basis_order, ω, x1 - x2)
            U = Uout(x1 - x2)

            U = U .* (bar_numdensity * pair_corr(x1,s1,x2,s1) * sin(rθφ[2]) * rθφ[1]^2)
            return [
                    t_diags[1][n] * sum(U[nd,n] * discrete_coefs(x2)[nd] for nd in 1:len)
            for n in 1:len]
        end

        return fun
    end

    rs = LinRange(0.0, R - outer_radius(s1), mesh_points);
    θs = LinRange(0.0, π, mesh_points);

    test_ker = kernel_function(SVector(mean(rs),θs[1]))
    (v2,err) = hcubature(test_ker, SVector(0.0,0.0,-π), SVector(R-outer_radius(s1),π,π)
        ; rtol=rtol, maxevals = maxevals
    );

    println("The estimated max coefficient of the integrated kernel is:")
    println("I = ", maximum(abs.(v2)) )
    println("with an estimated error of: ")
    println("E = ", err)

    relative_residues = [
        begin
            ker = kernel_function([r,θ])
            K = hcubature(ker, SVector(0.0,0.0,-π), SVector(R-outer_radius(s1),π,π)
                ; rtol=rtol, maxevals = maxevals
            )[1];
            inc = incident_coefficients([r,θ]);
            norm(- discrete_coefs(rθφ2xyz([r,θ,0.0])) + inc + K) / norm(inc)
        end
    for r in rs, θ in θs]

    return (mean(relative_residues),maximum(relative_residues))
end


text:
using Flux

function get_optimizer(parsed::Dict, prefix="")
    opt_key = prefix == "" ? "opt" : join([prefix, "opt"], "_")
    opt_string = parsed[opt_key]
    get_optimizer(opt_string, parsed, prefix)
end

function get_optimizer(opt_string, parsed::Dict, prefix)

    if opt_string == "Auto+Descent"
        auto = _init_optimizer(Auto, parsed, "$(prefix)_auto")
        sgd = get_optimizer("Descent", parsed, prefix)
        (auto, sgd)
    elseif contains(opt_string, '+')
        opt_1_str, opt_2_str = split(opt_string, '+')
        (get_optimizer(opt_1_str, parsed, prefix), get_optimizer(opt_2_str, parsed, prefix))
    else
        opt_type = if opt_string == "Auto"
            Auto
        else
            getproperty(Flux, Symbol(opt_string))
        end
        _init_optimizer(opt_type, parsed, prefix)
    end
end

_init_optimizer(opt, ::Dict) =
    throw("$(string(opt)) optimizer initialization not found.")

function _init_optimizer(opt_type::Union{Type{Descent}, Type{ADAGrad}, Type{ADADelta}}, parsed::Dict, prefix = "")
    eta_str = prefix == "" ? "eta" : join([prefix, "eta"], "_")
    try
        η = parsed[eta_str]
        opt_type(η)
    catch
        throw("$(opt_type) needs: $(eta_str) (float).")
    end
end

function _init_optimizer(opt_type::Union{Type{RMSProp},
                                         Type{Momentum},
                                         Type{Nesterov}},
                         parsed::Dict, prefix="")

    eta_str = prefix == "" ? "eta" : join([prefix, "eta"], "_")
    rho_str = prefix == "" ? "rho" : join([prefix, "rho"], "_")
    try
        η = parsed[eta_str]
        ρ = parsed[rho_str]
        opt_type(η, ρ)
    catch
        throw("$(opt_type) needs: $(eta_str) (float), and $(rho_str) (float).")
    end
end

function _init_optimizer(opt_type::Union{Type{ADAM},
                                         Type{RADAM},
                                         Type{NADAM},
                                         Type{AdaMax},
                                         Type{AMSGrad}},
                         parsed::Dict, prefix="")

    eta_str = prefix == "" ? "eta" : join([prefix, "eta"], "_")
    beta_str = prefix == "" ? "beta" : join([prefix, "beta"], "_")
    try
        η = parsed[eta_str]
        β = if beta_str ∈ keys(parsed)
            parsed[beta_str]
        else
            (parsed[join([beta_str, "m"], "_")],
             parsed[join([beta_str, "v"], "_")])
        end
        opt_type(η, β)
    catch
        throw("$(opt_type) needs: $(eta_str) (float), and $(beta_str) ((float, float)), or (beta_m, beta_v)).")
    end
end

function _init_optimizer(opt_type::Union{Type{Auto}}, parsed::Dict, prefix="")
    α_str = prefix == "" ? "eta" : join([prefix, "eta"], "_")
    α_init_str = prefix == "" ?  "alpha_init" : join([prefix, "alpha_init"], "_")
    try
        α = parsed[α_str]
        α_init = parsed[α_init_str]
        opt_type(α, α_init)
    catch
        throw("$(opt_type) needs: $(α_str) (float), and $(α_init_str) (float).")
    end
end

function get_linear_learner(parsed::Dict,
                            feature_size::Int,
                            num_actions::Int,
                            num_demons::Int,
                            num_tasks::Int,
                            prefix="",
                            feature_projector = nothing)

    learner_key = prefix == "" ? "learner" : join([prefix, "learner"], "_")
    opt = get_optimizer(parsed, prefix)
    lu = get_learning_update(parsed, opt, prefix)
    w_init_key = prefix == "" ? "w_init" : join([prefix, "w_init"], "_")
    w_init = w_init_key in keys(parsed) ? parsed[w_init_key] : 0

    learner_str = parsed[learner_key]
    demon_learner = if learner_str ∈ ["Q", "QLearner", "q"]
        LinearQLearner(lu,
                       feature_size,
                       num_actions,
                       num_demons,
                       w_init)
    elseif learner_str ∈ ["SR", "SRLearner", "sr"]
        SRLearner(lu,
                  feature_size,
                  num_demons,
                  num_actions,
                  num_tasks,
                  feature_projector,
                  w_init)
    elseif learner_str ∈ ["GPI", "gpi"]
        GPI(lu,
            feature_size,
            num_demons,
            num_actions,
            num_tasks,
            feature_projector,
            w_init)
    elseif learner_str ∈ ["LSTD", "LSTDLearner"]
        eta_str = prefix == "" ? "eta" : join([prefix, "eta"], "_")
        LSTDLearner(lu,
                    parsed[eta_str],
                    feature_size,
                    num_actions,
                    num_demons)
    elseif learner_str ∈ ["NoLearner", "nolearner"]
        NoLearner(collect(1:num_actions), num_demons)
    else
        throw(ArgumentError("Not a valid demon learner"))
    end
end

function get_linear_learner(parsed::Dict,
                            feature_size,
                            num_actions,
                            demons,
                            prefix="",
                            feature_projector=nothing)

    # Don't want to change call now. num_tasks only exists for SRHorde demons.
    num_tasks = if !hasproperty(demons, :num_tasks)
        0
    else
        demons.num_tasks
    end
    return get_linear_learner(parsed::Dict,
        feature_size,
        num_actions,
        length(demons),
        num_tasks,
        prefix,
        feature_projector)
end

function get_learning_update(parsed::Dict, opt, prefix="")
    lu_key = prefix == "" ? "update" : join([prefix, "update"], "_")
    lu = getproperty(Curiosity, Symbol(parsed[lu_key]))
    _init_learning_update(lu, opt, parsed, prefix)
end

_init_learning_update(lu_type, args...) =
    throw("$(string(lu_type)) does not have an init function.")

function _init_learning_update(lu_type::Union{Type{TabularRoundRobin}}, args...)
    lu_type()
end

function _init_learning_update(lu_type::Union{Type{TB}},
                               opt,
                               parsed::Dict,
                               prefix)
    λ_str = prefix == "" ? "lambda" : join([prefix, "lambda"], "_")
    try
        λ = parsed[λ_str]
        lu_type(lambda=λ, opt=opt)
    catch
        throw("$(lu_type) needs: $(λ_str) (float).")
    end
end

function _init_learning_update(lu_type::Union{Type{InterestTB}},
                                opt,
                                parsed::Dict,
                                prefix)
    λ_str = prefix == "" ? "lambda" : join([prefix, "lambda"], "_")
    interest_set_str = prefix == "" ? "interest_set" : join([prefix, "interest_set"], "_")
    norm_interest_str = prefix == "" ? "interest_set" : join([prefix, "normalize_interest"], "_")

    try
        lambda = parsed[λ_str]
        interest_set = parsed[interest_set_str]
        norm_interest = parsed[norm_interest_str]

        hdpi = get_hdpi(parsed, interest_set)
        lu_type(lambda=lambda, opt=opt, hdpi=hdpi, normalize_interest=norm_interest)

    catch e
        throw(e)
        # throw("$(lu_type) needs: $(interest_set_str) (string)")
    end
end

function _init_learning_update(lu_type::Union{Type{TB}, Type{PriorESARSA}, Type{PriorTB}},
                                opt,
                                parsed::Dict,
                                prefix)
    λ_str = prefix == "" ? "lambda" : join([prefix, "lambda"], "_")
    try
        λ = parsed[λ_str]

        if haskey(parsed, "clip_threshold")
            clip_threshold = parsed["clip_threshold"]
        else
            clip_threshold = Inf64
            println("Setting clip threshold to inf")
        end

        lu_type(lambda=λ, opt=opt, clip_threshold=clip_threshold)
    catch
        throw("$(lu_type) needs: $(λ_str) (float).")
    end
end

function _init_learning_update(lu_type::Union{Type{ETB},Type{EmphESARSA},Type{StateETB}},
                                opt,
                                parsed::Dict,
                                prefix)
    λ_str = prefix == "" ? "lambda" : join([prefix, "lambda"], "_")
    interest_set_str = prefix == "" ? "interest_set" : join([prefix, "interest_set"], "_")
    norm_interest_str = prefix == "" ? "interest_set" : join([prefix, "normalize_interest"], "_")

    # try
    lambda = parsed[λ_str]
    interest_set = parsed[interest_set_str]
    norm_interest = parsed[norm_interest_str]

    if haskey(parsed, "emphasis_clip_threshold")
        clip_threshold = parsed["emphasis_clip_threshold"]
    else
        clip_threshold = Inf64
        println("Setting clip threshold to inf")
    end

    hdpi = get_hdpi(parsed, interest_set)
    lu_type(lambda=lambda, opt=opt, hdpi=hdpi, clip_threshold=clip_threshold, normalize_interest=norm_interest)

    # catch e
    #     throw(e)
    # throw("$(lu_type) needs: $(interest_set_str) (string)")
    # end
end

function _init_learning_update(lu_type::Union{Type{ESARSA}, Type{SARSA}},
                                opt,
                                parsed::Dict,
                                prefix)
    λ_str = prefix == "" ? "lambda" : join([prefix, "lambda"], "_")
    trace_str = prefix == "" ? "trace" : join([prefix, "trace"], "_")
    trace = try
        getproperty(Curiosity, Symbol(parsed[trace_str]))
    catch
        throw("$(lu_type) needs trace type")
    end

    try
        λ = parsed[λ_str]
        lu_type(lambda=λ, opt=opt, trace=trace())
    catch
        throw("$(lu_type) needs: $(λ_str) (float).")
    end
end

function get_exploration_strategy(parsed, action_set)
    if parsed["exploration_strategy"] == "epsilon_greedy"
        EpsilonGreedy(parsed["exploration_param"])
    elseif parsed["exploration_strategy"] == "epsilon_greedy_decay"
        ϵGreedyDecay(
            parsed["ϵ_range"],
            parsed["decay_period"],
            parsed["warmup_steps"],
            action_set
        )
    else
        throw(ArgumentError("Not a Valid Exploration Strategy"))
    end
end

Base.@kwdef struct NoLearner <: Learner
    action_set::Array{Int,1}
    num_demons::Int
end

Curiosity.update!(learner::NoLearner, args...) = nothing

Base.get(π::NoLearner; state_t, action_t, kwargs...) =
    get_action_probs(π, state_t, nothing)[action_t]

function Curiosity.get_action_probs(l::NoLearner, features, state)
    return ones(l.num_demons,length(action_set)) ./ length(action_set)
end
(l::NoLearner)(ϕ,a) = zeros(l.num_demons)
(l::NoLearner)(ϕ) = [zeros(l.num_demons) for a in l.action_set]
function zero_eligibility_traces!(l::NoLearner) nothing end

function get_hdpi(parsed, interest_set)
    if (interest_set == "ttmaze")
        if (parsed["exploring_starts"] == false)
            @load "./src/data/dpi/ttmaze_exploring_starts_false.jld2" hdpi
            return hdpi
        end
    elseif (interest_set == "oned_tmaze")
        if (parsed["exploring_starts"] == "beg")
            @load "./src/data/dpi/oned_tmaze_beg.jld2" hdpi
            return hdpi
        end
    elseif (interest_set == "2dOpenWorld_center")
        @load "./src/data/dpi/2dOpenWorld_center.jld2" hdpi
        return hdpi
    end
    throw("no valid hdpi for string $(interest_set)")
end

ext: jl

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
function runexperiment(trialtype::String, args...; kwargs...)
	runexperiment(ExperimentType(trialtype), args...; kwargs...)
end

function runexperiment(::ExperimentType{:new}, name::AbstractString,
					   mp::ModelParams, delgraph::Array{Int64, 2},
					   tp::TrialParams, input;
					   numprocs = 8, execloc="./")
	writemparams(mp, name)
	writetparams(tp, name)
	writeinput(input, name)
	writedelgraph(delgraph, name)
	run(`mpirun -np $(numprocs) $(execloc)runtrial-mpi 0 $(name) $(name)`)
	spikes = loadspikes(name)
	inputtimes, inputids = loadinputtimes(name)
	synapses = loadsynapses(name)
	ExperimentOutput(name, name, tp, input, spikes, inputtimes, inputids, synapses)
end

function runexperiment(::ExperimentType{:new}, name::AbstractString,
					   m::Model,
					   tp::TrialParams, input; weights=:none,
					   numprocs = 8, execloc="./")
    println("dispatched here.")
	writemparams(m.params, name)
	writedelgraph(m.delgraph, name)
	writesyngraph(m.syngraph, name)
	writeneurons(m.neurons, name)
	writetparams(tp, name)
    if weights == :none
	    writeinput(input, name)
    elseif typeof(weights) <: Number
        println("and, again, here")
        writeinput(input, name; weight=weights)
    else
        writeinput(input, weights, name)
    end
	run(`mpirun -np $(numprocs) $(execloc)runtrial-mpi 0 $(name) $(name)`)
	spikes = loadspikes(name)
	inputtimes, inputids = loadinputtimes(name)
	synapses = loadsynapses(name)
	ExperimentOutput(name, name, tp, input, spikes, inputtimes, inputids, synapses)
end

function runexperiment(::ExperimentType{:resume}, modelname::AbstractString,
					   trialname::AbstractString, tp::TrialParams, input;
					   numprocs = 8, execloc="./")
	writetparams(tp, trialname)
	writeinput(input, trialname)
	run(`mpirun -np $(numprocs) $(execloc)runtrial-mpi 1 $(modelname) $(trialname)`)
	spikes = loadspikes(trialname)
	inputtimes, inputids = loadinputtimes(trialname)
	synapses = loadsynapses(trialname)
	ExperimentOutput(modelname, trialname, tp, input, spikes, inputtimes, inputids, synapses)
end

function stageexperiment(name::AbstractString, m::Model, tp::TrialParams, input;
						 dest="./")
	writemparams(m.params, dest*name)
	writedelgraph(m.delgraph, dest*name)
	writesyngraph(m.syngraph, dest*name)
	writeneurons(m.neurons, dest*name)
	writetparams(tp, dest*name)
	writeinput(input, dest*name)
end

function pgsearch(name::AbstractString;
				  basesize=3, threshold=18.0, duration=0.100,
				  maxgroups=50000, mingrouplen=3,
				  numprocs=8, execloc="/home/dahlbom/research/delnet/")
	run(`mpirun -np $(numprocs) $(execloc)pgsearch-mpi $(name) $(basesize) $(threshold) $(duration) $(maxgroups)`)
	spikes = loadspikes(name*"pg")
	groups = []
	t1 = 0.0
	t2 = t1 + duration
	idx1 = 1
	idx2 = 0
	for k ∈ 1:length(spikes)
		if spikes[k].t > t2
			idx2 = k-1
			if idx2-idx1+1 >= mingrouplen
				push!(groups, [Spike(s.n,s.t-spikes[idx1].t)
							   for s ∈ spikes[idx1:idx2]])
			end
			t1 += duration
			t2 += duration
			idx1 = idx2
		end
	end
	return groups
end

struct DenseGraph{T}
	contacts::Array{Int64,1}
	vals::Array{T,1}
	startidcs::Array{Int64,1}
	counts::Array{Int64,1}
end

function densegraph(g; offset=1)
	n = size(g)[1]
	numcontacts = zeros(Int64,n)
	startidcs = ones(Int64,n)

	# Set up memory
	for c ∈ 1:n
		for r ∈ 1:n
			if g[r,c] != 0
				numcontacts[r] += 1
			end
		end
	end
	for i ∈ 2:n
		startidcs[i] = startidcs[i-1] + numcontacts[i-1]
	end

	# Copy over contents into dense structure
	numgraphentries = sum(numcontacts)
	densegraph = zeros(Int64, numgraphentries)
	graphvals = zeros(typeof(g[1,1]), numgraphentries)
	counts = zeros(Int64,n)
	for r ∈ 1:n
		for c ∈ 1:n
			if g[r,c] != 0
				densegraph[startidcs[r]+counts[r]] = c
				graphvals[startidcs[r]+counts[r]] = g[r,c]
				counts[r] += 1
			end
		end
	end
	DenseGraph(densegraph, graphvals, startidcs, numcontacts)
end

function densifysynapses(ss, n; offset=1)
	numentries = length(ss)
	counts = zeros(Int64, n)
	startidcs = ones(Int64, n)
	contacts = zeros(Int64, numentries)
	vals = zeros(Float64, numentries)
	runningcounts = zeros(Int64, n)
	for s ∈ ss
		counts[s.source+1] += 1
	end
	for i ∈ 2:n
		startidcs[i] = startidcs[i-1] + counts[i-1]
	end

	for s ∈ ss
		contacts[startidcs[s.source+1] + runningcounts[s.source+1]] = s.dest+offset
		vals[startidcs[s.source+1] + runningcounts[s.source+1]] = s.strength
		runningcounts[s.source+1] += 1
	end
	DenseGraph(contacts, vals, startidcs, counts)
end

max_stars_repo_name: analogouscircuit/DelayNetTools

max_issues_repo_name: analogouscircuit/DelayNetTools

max_forks_repo_name: analogouscircuit/DelayNetTools


text:
"F, the residual level Fct"
function residual( z::Array{T,1}, alpha::Array{T,1}, c::SLiSeFilters.Config ):T
    X = calcX(z, c);
    W = calcW(z, c);
    Z = calcZ(z, c);
    Y = calcY(z, c);
    theta = calctheta(z, c);

    XZa = (X - Z) * alpha
    WYa = (W - Y) * conj(alpha)

    return (2*real(alpha' * (XZa + WYa -2*theta)) .+ real(caseh_h( c )))[1]
end

"The gradient of F"
function Fgradient( z::Array{T,1}, alpha::Array{T,1}, c::SLiSeFilters.Config )
    X = calcX(z, c);
    Z = calcZ(z, c);
    theta = calctheta(z, c);
    Wc = calcW(conj(z), c);
    Yc = calcY(conj(z), c);

    da = 4*(alpha' * (X - Z) + transpose(alpha) * (Wc - Yc) - theta' );

    dYc = calcDeltaYc( z, c );
    dWc = calcDeltaWc( z, c );
    dX = calcDeltaX( z, c );
    dZ = calcDeltaZ( z, c );
    dTheta = calcDeltaTheta( z, c );

    dz = 4*(alpha' * (dX - dZ) + transpose(alpha) * (dWc - dYc) - dTheta') .* transpose(alpha); # elementwise multiplication

    # we conjugate as we want column vectors
    return dz',da'
end

# the approximation to the hessian required for levMar
function hessianApprox( z::Array{T,1}, a::Array{T,1}, c::SLiSeFilters.Config )

    ddX = calcddX( z, c );
    ddZ = calcddZ( z, c );

    X = calcX( z, c);
    Z = calcZ( z, c);

    dX = calcDeltaX( z, c );
    dZ = calcDeltaZ( z, c );

    H2 = 2*[  (diagm(a)' * ( ddX - ddZ) * diagm(a)) (diagm(a)' * (dX - dZ)')  ;
              ( ( dX -  dZ ) * diagm(a) )          (X - Z)                  ];

    dz, da = Fgradient( z, a, c );

    H = H2;
    d = (0.5.*[dz;da]);
    d = (d[:])

    return H, d;
end


text:
"""

## Theoretical performance of a BHA

Function to compute the theoretical performance of a Bottom Hole Assembly (BHA)
in terms of building or dropping rate.

### Function
```julia
theoretical_performance!(data::Dict{Symbol, Any})

or

tp!(data::Dict{Symbol, Any})

```
### Arguments
```julia
* `data`      : Dictionary with input (and updated) data.
```
### Results
```julia
* `params`    : Object with all computational results.
```
"""
function tp!(data::Dict{Symbol, Any})
  if !(:wobrange in keys(data))
    println("Please specify a WOB range, e.g. :wobrange = 5:10:25.")
    data[:terminated] = true
    return
  end
  if !(:inclrange in keys(data))
    println("Please specify an inclination range, e.g. inclrange = 20:5:50.")
    data[:terminated] = true
    return
  end

  if :wob in keys(data)
    if typeof(data[:wobrange]) == StepRange{Int64, Int64}
      data[:wob] = float(data[:wobrange].start)
    else
      data[:wob] = float(data[:wobrange][1])
    end
  else
    if typeof(data[:wobrange]) == StepRange{Int64, Int64}
      push!(data, :wob =>  float(data[:wobrange].start))
    else
      push!(data, :wob =>  float(data[:wobrange][1]))
    end
  end

  if :inclination in keys(data[:welltrajectory])
    if typeof(data[:inclrange]) == StepRange{Int64, Int64}
      data[:welltrajectory][:inclination] = float(data[:inclrange].start)
    else
      data[:welltrajectory][:inclination] = float(data[:inclrange][1])
    end
  else
    if typeof(data[:inclrange]) == StepRange{Int64, Int64}
      push!(data[:welltrajectory], :inclination => float(data[:inclrange].start))
    else
      push!(data[:welltrajectory], :inclination => float(data[:inclrange][1]))
    end
  end
  push!(data[:welltrajectory], :dalpha => 0.0)

  checktrajectoryinput!(data)
  if :terminated in keys(data)
    return
  end
	println("Trajectory input data checks completed!")

  checkinput!(data)
  if :terminated in keys(data)
    return
  end
	println("BHA input data checks completed!\n")

  wobrange = data[:wobrange]
  inclrange = data[:inclrange]
  heading = data[:welltrajectory][:heading]
  holediameter = data[:welltrajectory][:diameter]

  if length(wobrange) > 0 && length(inclrange) > 0 && !isdefined(data, :welltrajectory)
    indxtab = Array{Tuple}(undef, length(wobrange)*length(inclrange))
    for (i, wob) in enumerate(wobrange)
      for (j, incl) in enumerate(inclrange)
        #wobstring = string(Int(wob/klbf))
        #inclstring = string(Int(incl/°))
        wobstring = string(Int(wob))
        inclstring = string(Int(incl))
        cd(data[:pdir]) do
          if isdir("tp0/$wobstring/$inclstring")
            rm("tp0/$wobstring/$inclstring", recursive=true)
          end
          (!isdir("tp0")) && mkdir("tp0")
          (!isdir("tp0/$wobstring")) && mkdir("tp0/$wobstring")
          (!isdir("tp0/$wobstring/$inclstring")) && mkdir("tp0/$wobstring/$inclstring")
          if isdir("tp/$wobstring/$inclstring")
            rm("tp/$wobstring/$inclstring", recursive=true)
          end
          (!isdir("tp")) && mkdir("tp")
          (!isdir("tp/$wobstring")) && mkdir("tp/$wobstring")
          (!isdir("tp/$wobstring/$inclstring")) && mkdir("tp/$wobstring/$inclstring")
          indxtab[(i - 1) * length(inclrange) + j] =
            (data, wob, incl, holediameter, heading)
        end
      end
    end
    pmap(tprun, indxtab)
    println()
  end
end


text:
# This script will generate a plot similar to random.png

using DaylightPlots, Plots

pyplot()

t = DateTime(2018,01,01):Minute(1):DateTime(2018,01,02)

x = randn(length(t))

daylight(t,x,loc=Point(42.35,-71.05),tz=tz"America/New_York",color=:black,leg=false,grid=false)

savefig("random.png")


text:
"""
    combine_indices(x, y)

Returns the combined axes of `x` and `y` for broadcasting operations.
"""
combine_indices(x::AbstractArray, y::AbstractArray) = combine_indices(axes(x), axes(y))
function combine_indices(x::Tuple, y::Tuple)
    return (combine_index(first(x), first(y)), combine_indices(tail(x), tail(y))...)
end
combine_indices(x::Tuple{Any}, y::Tuple{}) = (first(x),)
combine_indices(x::Tuple{}, y::Tuple{Any}) = (first(y),)
combine_indices(x::Tuple{}, y::Tuple{}) = ()

"""
    combine_index(x, y)

Returns the combination of `x` and `y`, creating a new index. New subtypes of
`AbstractIndex` should implement a `combine_index` method.
"""
function combine_index(x::Index, y::Index)
    return Index{combine_names(x, y)}(combine_keys(x, y), combine_values(x, y))
end
function combine_index(x::SimpleIndex, y::SimpleIndex)
    return SimpleIndex{combine_names(x, y)}(combine_values(x, y))
end
function combine_index(x::AbstractIndex, y::AbstractIndex)
    error("`combine_index` must be defined for new subtypes of AbstractIndex.")
end

"""
    combine_values(x, y)

Returns the combination of the values of `x` and `y`, creating a new index. New
subtypes of `AbstractIndex` may implement a unique `combine_values` method if
needed. Default behavior is to use the return of `promote_rule(x, y)` for the
type of the combined values.
"""
function combine_values(x::AbstractIndex, y::AbstractIndex)
    return combine_values(promote_values_rule(x, y), values(x), values(y))
end
combine_values(::Type{T}, x, y) where {T<:AbstractUnitRange} = T(x)

"""
    combine_keys(x::AbstractIndex, y::AbstractIndex)

Returns the combination of the keys of `x` and `y`, creating a new index. New
subtypes of `AbstractIndex` may implement a unique `combine_keys` method if
needed. Default behavior is to use the return of `promote_rule(x, y)` for the
type of the combined keys.
"""
function combine_keys(x::AbstractIndex, y::AbstractIndex)
    return combine_keys(promote_keys_rule(x, y),keys(x), keys(y))
end

combine_keys(::Union{}, x, y) = combine_keys(typeof(x), x, y)
combine_keys(::Type{T}, x, y) where {T<:Union{OneTo,OneToRange}} = T(length(x))
combine_keys(::Type{T}, x, y) where {T<:AbstractUnitRange} = T(first(x), last(x))
function combine_keys(::Type{T}, x, y) where {T<:Union{StepRange,AbstractStepRange}}
    return T(first(x), step(x), last(x))
end
function combine_keys(::Type{T}, x, y) where {T<:Union{LinRange,AbstractLinRange}}
    return T(first(x), last(x), length(x))
end
function combine_keys(::Type{T}, x, y) where {T<:Union{StepRangeLen,AbstractStepRangeLen}}
    return T(first(x), step(x), length(x), x.offset)
end
combine_keys(::Type{T}, x, y) where {T<:AbstractVector} = copy(x)

"""
    combine_names(a, b)

Returns the combined name of `a` and `b`. The standard rules are:

* nothing + nothing = name
* nothing + name = name
* name + nothing = name
* name1 + name2 = name1

## Examples
```jldoctest
julia> a, b, n, = Index{:a}(1:10), Index{:b}(1:10), Index(1:10);

julia> combine_names(a, b)
:a

julia> combine_names(a, n)
:a

julia> combine_names(b, n)
:b

julia> combine_names(n, n)

julia> combine_names(b, a)
:b
```
"""
combine_names(a::Union{Symbol,Nothing}, b::AbstractIndex) = combine_names(a, dimnames(b))
combine_names(a::AbstractIndex, b::Union{Symbol,Nothing}) = combine_names(dimnames(a), b)
combine_names(a::AbstractIndex, b::AbstractIndex) = combine_names(dimnames(a), dimnames(b))

combine_names(a::Symbol, b::Symbol) = a
combine_names(::Nothing, b::Symbol) = b
combine_names(a::Symbol, ::Nothing) = a
combine_names(::Nothing, ::Nothing) = nothing


text:
##
##  function detect_polyA_in_a_string(
##    fq_seq::String,
##    minimum_polyA_length::Int64,
##    maximum_non_A_symbols::Int64;
##    maximum_search_fragment_length::Int64=50,
##    debug=false
##    )::Bool
##

minimum_polyA_length=10
min_l = Int64(minimum_polyA_length / 2)
re = Regex("(A+[GTC])?(A{$min_l,})([GTC]A+)?")

seq = "AGAGAATTACAAATCAGAAGGGGAAGATTGACTGTTTAATAAAATGTGCTGGGAATCAATGGCAAATATACAATAAAATAAAATTATATTTTTATTCATATTATACCCCAAATAAATTTCAGATATACATAAAACCCAAAGTAAAATATT"
@test !detect_polyA_in_a_string(seq,10,1,re,debug=false)

seq = "GAATGTATGGTAGGAATGTATTCTCTTGTAGGAATGTAAATCTGTATTAAAAGGGGGTCCAAGCCAGGCCCCCAGGTCTTCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACAAAAAAAA"
@test detect_polyA_in_a_string(seq,10,1,re,debug=false)

seq = "AGAGAATTACAAATCAGAAGGGGAAGATTGACTGTTTAATAAAATGTGCTGGGAATCAATGGCAAATATACAATAAAATAAAATTATATTTTTATTCATATTATACCCCAAATAAATTTCAGATATACATAAAACCCAAAGTAAAATATT"
@test detect_polyA_in_a_string(seq,10,3,re,debug=false)

seq = "AGAGAATTACAAATCAGAAGGGGAAGATTGACTGTTTAATAAAATGTGCTGGGAATCAATGGCAAATATACAATAAAATAAAATTATATTTTTATTCATATTATACCCCAAATAAATTTCAGATATACATAAAACCCAAAGTAAAATATT"
@test !detect_polyA_in_a_string(seq,10,1,re,debug=false)

seq = "TCAAAAAAATAATATGGTAATAATAATAAAAGCAGTGCCCAGAGAACAAGGCTTGTTGGCTGTTCCACCCCAGGGGGCCCCTTGCACAGGCGGTGCCATCTCTGCCTCCCAAAGCTCTAAGAGCCACTGTCCCCCATCCCAAGAGA"
@test !detect_polyA_in_a_string(seq,10,1,re,debug=false)

seq = "TCAAAAAAATAATATGGTAATAATAATAAAAGCAGTGCCCAGAGAACAAGGCTTGTTGGCTGTTCCACCCCAGGGGGCCCCTTGCACAGGCGGTGCCATCTCTGCCTCCCAAAGCTCTAAGAGCCACTGTCCCCCATCCCAAGAGA"
@test !detect_polyA_in_a_string(seq,10,1,re,debug=false)

seq = "GGGGCGCGGCGGCGCGGGGGGGGAAAAAAAAACCCCCCCCCCCCCCCGGGGGGCCCC"
@test !detect_polyA_in_a_string(seq,10,1,re,debug=false)

seq = "GGGGCGCGGCGGCGCGGGGGGGGAAAATAAAAACCCCCCCCCCCCCCCGGGGGGCCCC"
@test detect_polyA_in_a_string(seq,10,1,re,debug=false)

seq = "AAGAAAAAAACAAACAACCCCATCAAAAAGTGGGCAAAGGATATAAATAGACAC"
@test detect_polyA_in_a_string(seq,10,1,re,debug=false)

seq = "TCAGGAGATCGAGACCATCCTGGCTAACATGGTGAAACCCCGTCTCTACTAAAAATAACAAAAAAAAA"
@test detect_polyA_in_a_string(seq,10,1,re,debug=false)

seq = "TGTGTGTGTGTGTGGTAAAAAAAAAAGT"
@test detect_polyA_in_a_string(seq,10,1,re,maximum_search_fragment_length=20)

seq = "TGTGTGTGTGTGTGGTAAAAAAAAAA"
@test detect_polyA_in_a_string(seq,10,1,re,maximum_search_fragment_length=20)

seq = "TGTGTGTGTGTGTGGTAAAAAAAAATAGGGG"
@test detect_polyA_in_a_string(seq,10,1,re,maximum_search_fragment_length=20)

seq = "TGTGTGTGTGTGTGGATAAAAAAAAA"
@test detect_polyA_in_a_string(seq,10,1,re,maximum_search_fragment_length=20)

seq = "TGTGTGTGTGTGTGGAATAAAAAAAAA"
@test detect_polyA_in_a_string(seq,10,1,re,maximum_search_fragment_length=20)

seq = "TGTGTGTGTGTGTGGAAAAAAAAATAA"
@test detect_polyA_in_a_string(seq,10,1,re,maximum_search_fragment_length=20)

seq = "TGTGTGTGTGTGTGGAATAAAAAAAAATAA"
@test detect_polyA_in_a_string(seq,10,1,re,maximum_search_fragment_length=20)

seq = "AAAAAAAAAAAAAAAAAAAATGGAACGCAGGGCAGGAACTCGTATTTGGGGGGAGATGGG"
@test detect_polyA_in_a_string(seq,10,1,re,debug=false)

seq = "TGTGTGTGTGTGTGGAAAAAAAAAT"
@test !detect_polyA_in_a_string(seq,10,1,re,maximum_search_fragment_length=20)

seq = "TGTGTGTGTGTGTGGTAAAAAAAAA"
@test !detect_polyA_in_a_string(seq,10,1,re,maximum_search_fragment_length=20)

seq = "TGTGTGTGTGTGTGGAAAAAAAAA"
@test !detect_polyA_in_a_string(seq,10,1,re,maximum_search_fragment_length=20)

seq = "AAAAAAAAAAAAAAATGTGTGTGTGTGTGGAAAAAAAAA"
@test !detect_polyA_in_a_string(seq,10,1,re,maximum_search_fragment_length=20)

seq = "AAAAAAAAAAAAAAA"
@test detect_polyA_in_a_string(seq,10,1,re,maximum_search_fragment_length=20)

seq = "AAAAAAAAA"
@test !detect_polyA_in_a_string(seq,10,1,re,maximum_search_fragment_length=20)

max_stars_repo_name: laurabiggins/polyAterminus

max_issues_repo_name: laurabiggins/polyAterminus

max_forks_repo_name: laurabiggins/polyAterminus


text:
module PathCompare
    export Distance

    # http://en.wikipedia.org/wiki/Dynamic_time_warping
    function DTWDistance{T}(s::Array{T}, t::Array{T})
        n, m = size(s, 1), size(t, 1)
        DTW = zeros(n+1,m+1)
        DTW[:,1], DTW[1,:] = Inf, Inf
        DTW[1,1] = 0

        for i=2:n+1
            for j=2:m+1
                cost = sum((s[i-1,:] - t[j-1,:]).^2)

                DTW[i,j] = cost + min(DTW[i-1, j  ], #insertion
                                      DTW[i  , j-1], #deletion
                                      DTW[i-1, j-1]) #match
            end
        end

        return DTW[n+1,m+1]
    end

    # note. there is also a windowed version for better performance!  See the wikipedia article
    function DTWDistance{T}(s::Array{T}, t::Array{T}, w::Int64)
        n, m = size(s, 1), size(t, 1)
        w = max(w, abs(n-m))
        DTW = zeros(n+1,m+1)
        #println("sizeof DTW: ", size(DTW))
        DTW[:,:] = Inf
        DTW[1,1] = 0

        for i=2:n+1
            for j=max(2, i-w):min(m+1, i+w)
                #println("i=", i, "/", n, " j=", j, "/", m, " size(s)=", size(s), " size(t)=", size(t))
                cost = sum((s[i-1,:] - t[j-1,:]).^2)

                DTW[i,j] = cost + min(DTW[i-1, j  ], #insertion
                                      DTW[i  , j-1], #deletion
                                      DTW[i-1, j-1]) #match
                #println("cost is ", cost, " i=", i, "/", n, " j=", j, "/", m, " DTW=", DTW[i,j])
            end
        end

        return DTW[n+1,m+1]
    end

    # Point distance metrics

    function HaversineDistance{T}(lat1::T, lon1::T, lat2::T, lon2::T)
        REarth = 6371.
        lat, lon = (lat1-lat2)*π/180., (lon1-lon2)*π/180.
        lat1, lat2 = lat1*π/180., lat2*π/180.
        a = sin(lat/2)*sin(lat/2)+cos(lat1)*cos(lat2)*sin(lon/2)*sin(lon/2)
        d = 2*atan2(sqrt(a), sqrt(1-a))
        d = REarth*d

        return d
    end

    function meanHaversideDistance{T}(lat1::T, lon1::T, lat2::T, lon2::T)
        return mean(HaversideDistance{T}(lat1, lon1, lat2, lon2))
    end
end


text:
# LSL.jl: Julia interface for Lab Streaming Layer
# Copyright (C) 2019 Samuel Powell

# wrap.jl: automatically wrap C API using Clang.jl

# Note: in current wrap the following manual changes are made
# 1. Remove duplicate signature
# 2. Change signature of push/pull const *char to Cchar from Cstring
# 3. Guard Clong functions with @static if !Sys.iswindows() && Sys.WORD_SIZE == 64

@info("NOTE MANUAL CHANGES in deps/wrap.jl")

using Clang

const LIBLSL_INCLUDE = joinpath(@__DIR__, "..", "deps", "usr", "include") |> normpath
const LIBLSL_HEADERS = [joinpath(LIBLSL_INCLUDE, header) for header in ["lsl_constants.h", "lsl_c.h"]]

wc = init(; headers = LIBLSL_HEADERS,
            output_file = joinpath(@__DIR__, "liblsl_api.jl"),
            common_file = joinpath(@__DIR__, "liblsl_common.jl"),
            clang_includes = vcat(LIBLSL_INCLUDE, CLANG_INCLUDE),
            clang_args = ["-I", joinpath(LIBLSL_INCLUDE, "..")],
            header_wrapped = (root, current)->root == current,
            header_library = x->"liblsl",
            clang_diagnostics = true,
            )

run(wc)


text:
using Cryptography
using FactCheck

NumTestRepetitions = 100


text:
using Latexify
using Test

arr = [1 2; 3 4]

@test latexify(arr) ==
raw"\begin{equation}
\left[
\begin{array}{cc}
1 & 2 \\
3 & 4 \\
\end{array}
\right]
\end{equation}
"

arr = [1,2,:(x/y),4]

@test latexify(arr) ==
raw"\begin{equation}
\left[
\begin{array}{c}
1 \\
2 \\
\frac{x}{y} \\
4 \\
\end{array}
\right]
\end{equation}
"

@test latexify(arr; transpose=true) ==
raw"\begin{equation}
\left[
\begin{array}{cccc}
1 & 2 & \frac{x}{y} & 4 \\
\end{array}
\right]
\end{equation}
"

@test latexify((1.0, 2), (3, 4)) ==
raw"\begin{equation}
\left[
\begin{array}{cc}
1.0 & 3 \\
2 & 4 \\
\end{array}
\right]
\end{equation}
"

@test latexify(((1.0, 2), (3, 4))) ==
raw"\begin{equation}
\left[
\begin{array}{cc}
1.0 & 3 \\
2 & 4 \\
\end{array}
\right]
\end{equation}
"

@test latexify(:(x = [1 2] * [1, 2] * [1 2; 3 4])) ==
raw"$x = \left[
\begin{array}{cc}
1 & 2 \\
\end{array}
\right] \cdot \left[
\begin{array}{c}
1 \\
2 \\
\end{array}
\right] \cdot \left[
\begin{array}{cc}
1 & 2 \\
3 & 4 \\
\end{array}
\right]$"

@test latexify(:(x = $arr)) ==
raw"$x = \left[
\begin{array}{c}
1 \\
2 \\
\frac{x}{y} \\
4 \\
\end{array}
\right]$"


text:
##################################################
# Hessian of `f(::AbstractArray)::AbstractArray` #
##################################################

# hessian #
#---------#

"""
    ReverseDiff.hessian(f, input::AbstractArray, cfg::HessianConfig = HessianConfig(input))

Given `f(input::AbstractArray{<:Real})::Real`, return `f`s Hessian w.r.t. to the given
`input`.

Note that `cfg` can be preallocated and reused for subsequent calls.

If possible, it is highly recommended to use `ReverseDiff.HessianTape` to prerecord `f`.
Otherwise, this method will have to re-record `f`'s execution trace for every subsequent
call.
"""
function hessian(f, input::AbstractArray, cfg::HessianConfig = HessianConfig(input))
    ∇f = x -> gradient(f, x, cfg.gradient_config)
    return jacobian(∇f, input, cfg.jacobian_config)
end

# hessian! #
#----------#

"""
    ReverseDiff.hessian!(result::AbstractArray, f, input::AbstractArray, cfg::HessianConfig = HessianConfig(input))

    ReverseDiff.hessian!(result::DiffResult, f, input::AbstractArray, cfg::HessianConfig = HessianConfig(result, input))

Returns `result`. This method is exactly like `ReverseDiff.hessian(f, input, cfg)`, except
it stores the resulting Hessian in `result` rather than allocating new memory.

If `result` is a `DiffResults.DiffResult`, the primal value `f(input)` and the gradient
`∇f(input)` will be stored in it along with the Hessian `H(f)(input)`.
"""
function hessian!(result, f, input::AbstractArray, cfg::HessianConfig = HessianConfig(input))
    ∇f = x -> gradient(f, x, cfg.gradient_config)
    jacobian!(result, ∇f, input, cfg.jacobian_config)
    return result
end

function hessian!(result::DiffResult, f, input::AbstractArray,
                  cfg::HessianConfig = HessianConfig(result, input))
    ∇f! = (y, x) -> begin
        gradient_result = DiffResult(zero(eltype(y)), y)
        gradient!(gradient_result, f, x, cfg.gradient_config)
        result = DiffResults.value!(result, value(DiffResults.value(gradient_result)))
        return y
    end
    jacobian!(DiffResults.hessian(result), ∇f!,
              DiffResults.gradient(result), input,
              cfg.jacobian_config)
    return result
end

##########################
# Executing HessianTapes #
##########################

"""
    ReverseDiff.hessian!(tape::Union{HessianTape,CompiledHessian}, input)

Assuming `tape` represents a function of the form `f(::AbstractArray{<:Real})::Real`,
return the Hessian `H(f)(input)`.
"""
function hessian!(tape::Union{HessianTape,CompiledHessian}, input::AbstractArray)
    result = construct_result(output_hook(tape), input_hook(tape))
    hessian!(result, tape, input)
    return result
end

"""
    ReverseDiff.hessian!(result::AbstractArray, tape::Union{HessianTape,CompiledHessian}, input)

    ReverseDiff.hessian!(result::DiffResult, tape::Union{HessianTape,CompiledHessian}, input)

Returns `result`. This method is exactly like `ReverseDiff.hessian!(tape, input)`, except
it stores the resulting Hessian in `result` rather than allocating new memory.

If `result` is a `DiffResults.DiffResult`, the primal value `f(input)` and the gradient
`∇f(input)` will be stored in it along with the Hessian `H(f)(input)`.
"""
function hessian!(result::AbstractArray, tape::Union{HessianTape,CompiledHessian}, input::AbstractArray)
    seeded_forward_pass!(tape, input)
    seeded_reverse_pass!(result, tape)
    return result
end

function hessian!(result::DiffResult, tape::Union{HessianTape,CompiledHessian}, input::AbstractArray)
    seeded_forward_pass!(tape, input)
    seeded_reverse_pass!(DiffResult(DiffResults.gradient(result), DiffResults.hessian(result)), tape)
    result = DiffResults.value!(result, func_hook(tape)(input))
    return result
end

######################
# Hessian API Errors #
######################

const HESS_MULTI_ARG_ERR_MSG = "Taking the Hessian of a function with multiple arguments is not yet supported"

hessian(f, xs::Tuple, ::HessianConfig) = error(HESS_MULTI_ARG_ERR_MSG)
hessian(f, xs::Tuple) = error(HESS_MULTI_ARG_ERR_MSG)
hessian!(outs::Tuple, f, xs::Tuple, ::HessianConfig) = error(HESS_MULTI_ARG_ERR_MSG)
hessian!(outs::Tuple, f, xs::Tuple) = error(HESS_MULTI_ARG_ERR_MSG)


text:
export plot_avg_longterm_temps, plot_glacier_dataset

function plot_glacier_dataset(gdirs_climate, PDE_refs)
    if plots
        # gdirs_climate = (dates, gdirs, longterm_temps, annual_temps)
        gdirs = gdirs_climate[2]
        theme = Theme(fontsize=18, font="TeX Gyre Heros")
        set_theme!(theme)
        # Ice surface velocities
        figV = Makie.Figure(resolution=(1200, 1800))
        n = 5
        m = floor(Int, length(gdirs_climate[2])/n)
        axsV = [Axis(figV[i, j]) for i in 1:n, j in 1:m]
        hidedecorations!.(axsV)
        tightlimits!.(axsV)
        Vs, hm_Vs = [],[]
        for (i, ax) in enumerate(axsV)
            ax.aspect = DataAspect()
            glacier_gd = xr.open_dataset(gdirs[i].get_filepath("gridded_data"))
            nx = glacier_gd.y.size # glacier extent
            ny = glacier_gd.x.size # really weird, but this is inversed
            Δx = abs(gdirs[i].grid.dx)
            Δy = abs(gdirs[i].grid.dy)
            rgi_id = gdirs[i].rgi_id
            name = gdirs[i].name
            ax.title = "$name - $rgi_id"
            Vx = PDE_refs["V̄x_refs"][i]
            Vy = PDE_refs["V̄y_refs"][i]
            V = (abs.(Vx) .+ abs.(Vy))./2
            push!(Vs,V)
            push!(hm_Vs, Makie.heatmap!(ax, (1:nx).*Δx, (1:ny).*Δy, fillZeros(V), colormap=:speed))
        end
        minV = minimum(minimum.(Vs))
        maxV = maximum(maximum.(Vs))
        Colorbar(figV[2:3,m+1], limits=(minV,maxV), colormap=:speed, label="Ice surface velocity")
        Label(figV[0, :], text = "Glacier dataset", textsize = 30)
        display(figV)
        Makie.save(joinpath(root_plots, "glaciers_V.pdf"), figV, pt_per_unit = 1)

        # Ice thickness
        figH = Makie.Figure(resolution=(1200, 1800))
        axsH = [Axis(figH[i, j]) for i in 1:n, j in 1:m]
        hidedecorations!.(axsH)
        tightlimits!.(axsH)
        Hs, hm_Hs = [],[]
        for (i, ax) in enumerate(axsH)
            ax.aspect = DataAspect()
            glacier_gd = xr.open_dataset(gdirs[i].get_filepath("gridded_data"))
            nx = glacier_gd.y.size # glacier extent
            ny = glacier_gd.x.size # really weird, but this is inversed
            Δx = abs(gdirs[i].grid.dx)
            Δy = abs(gdirs[i].grid.dy)
            rgi_id = gdirs[i].rgi_id
            name = gdirs[i].name
            ax.title = "$name - $rgi_id"
            H = glacier_gd.consensus_ice_thickness.data # initial ice thickness conditions for forward model
            push!(Hs, fillNaN(H))
            push!(hm_Hs, Makie.heatmap!(ax, H, colormap=:ice))
        end
        minH = minimum(minimum.(Hs))
        maxH = maximum(maximum.(Hs))
        Colorbar(figH[2:3,m+1], limits=(minH,maxH), colormap=:ice, label="Ice thickness (m)")
        Label(figH[0, :], text = "Glacier dataset", textsize = 30)
        display(figH)
        Makie.save(joinpath(root_plots, "glaciers_H.pdf"), figH, pt_per_unit = 1)
    end
end

function plot_avg_longterm_temps(climate, gdirs)
    if plots
        mean_longterm_temps, labels = [],[]
        for (climate_glacier, gdir) in zip(climate, gdirs)
            push!(mean_longterm_temps, climate_glacier["longterm_temps"].temp.data)
            push!(labels, gdir.rgi_id)
        end
        display(Plots.plot(mean_longterm_temps, label=permutedims(labels),
                            xlabel="Years", ylabel="Mean longterm air temperature (°C)", legend=:topright;
                            palette=palette(:tab10,15)))
    end
end

# TODO: correctly retrieve the glacier coordinates to plot them in `imshow` as an extent
function plot_monthly_map(climate, variable, year)
    climate = climate[variable].where(climate.time.dt.year == year, drop=true).groupby("time.month")
    fig_clim, ax_clim = pplt.subplots([1:6, 7:12], axheight=2)
    fig_clim.format(
        abc=true, abcloc="ul", suptitle= ("$year - monthly $variable")
    )
    for mon in 1:12
        if variable == "temp"
            m_var = ax_clim[mon].imshow(climate.mean()[mon], cmap="Thermal",
                    vmin=minimum(climate.mean().data), vmax=maximum(climate.mean().data)) # set common min max temp
        else
            m_var = ax_clim[mon].imshow(climate.sum()[mon], cmap="DryWet",
                    vmin=climate.sum().min().data, vmax=climate.sum().max().data) # set common min max precipitation
        end
        ax_clim[mon].set_title(Dates.monthname(mon))
        if(mon == 12)
            if variable == "temp"
                fig_clim.colorbar(m_var, label="Air temperature (°C)")
            else
                fig_clim.colorbar(m_var, label="Accumulated $variable (mm)")
            end
        end
    end
end


text:
export ActualModel

MLJBase.@mlj_model mutable struct ActualModel <: PDFEstimator
    dist::Distribution = Normal()
end

function MLJBase.fit(model::ActualModel, verbosity, X)
    n_vars = size(X, 2)

    fitresult = n_vars
    cache = nothing
    report = nothing
    fitresult, cache, report
end

function MLJBase.transform(model::ActualModel, fitresult, Xnew)
    xnew = collect(MLJBase.matrix(Xnew)')
    n_vars = fitresult

    ynew = model.dist isa UnivariateDistribution ? pdf.(model.dist, vec(xnew)) : pdf(model.dist, xnew)
    ynew = reshape(ynew, size(ynew, 1), 1)
    ynew = DataFrame(ynew, :auto)
end

function MLJBase.fitted_params(model::ActualModel, fitresult)
    n_vars = fitresult

    (
        n_vars=n_vars,
    )
end


text:
# Global set variables logicset, Υ hold
integer(x::AbstractString) = parse(Int, strip(x))
Base.range(x::AbstractString) = range(integer(match(r"^[0-9]+", x).match),
                               stop = integer(match(r"[0-9]+$", x).match))
# logicaloccursin(y::Symbol) = any( [ y ∈ keys(logicset[i]) for i in 1:length(logicset) ] )
# logicaloccursin(x::LogicalCombo, y::Symbol) = y ∈ keys(x)

"""
    logicalparse

Takes a command and parses it into logical calls that either assigning
additional feasible variable ranges or constrain the relationship between
variables.

```julia
logicalparse(command::String; logicset::LogicalCombo = LogicalCombo(), verbose=true)
logicalparse(command::String, logicset::LogicalCombo; ...)
logicalparse(commands::Array{String,1}, logicset::LogicalCombo; ...)
logicalparse(commands::Array{String,1}; ...)
```
### Arguments
* `verbose` : specifies to print to screen or not

### Operators
There are numerous operators available to be used in the logical parse command.

### Examples
```julia
julia> myset = logicalparse("a, b, c in 1:3")
a,b,c in 1:3             feasible outcomes 27 ✓          :3 3 3

julia> myset = logicalparse("a == b", myset)
a == b                   feasible outcomes 9 ✓           :1 1 2

julia> myset = logicalparse("a > c", myset)
a > c                    feasible outcomes 3 ✓           :3 3 1

julia> myset = logicalparse("c != 1", myset)
c != 1                   feasible outcomes 1 ✓✓          :3 3 2
```
"""
function logicalparse(
    commands::Array{String,1};
    logicset::LogicalCombo = LogicalCombo(),
    verbose=true)

  clearcounter()

  verbose && println("")

  any(occursin.(";", commands)) &&
    (commands = split.(commands, ";") |> Iterators.flatten |> collect .|> strip .|> string)

  for command in commands
      logicset = logicalparse(command, logicset=logicset, verbose=verbose)
  end
  logicset
end

function logicalparse(
    command::String;
    logicset::LogicalCombo = LogicalCombo(),
    verbose=true)

    clearcounter()

    # A vector of non-standard operators to ignore
    generators = "in"
    superset = "xor|iff|if|then|notthen"
    metaset = "XOR|IFF|IF|THEN|AND|OR"

    exclusions = join([generators, superset, metaset],"|")

    occursin(";", command) &&
      return logicalparse(string.(strip.(split(command, ";"))),
        logicset=logicset, verbose=verbose)

    verbose && print(command)
    # verbose && println("\n" * join(logicset.commands), "\n")

    (strip(command) == "" || command[1] == '#') && return logicset |>
      (x -> (((strip(command) != "") && push!(x.commands, command)); return x))

    occursin(r"∈|\bin\b", command) &&
      return definelogicalset(logicset, command) |>
        reportfeasible(command, verbose)  |>
        (x -> (push!(x.commands, command); return x))

    # Check for the existance of any symbols in logicset
    varcheck = eachmatch(r"[a-zA-Z][0-9a-zA-Z_.]*", command)

    # Replace command
    commandout = command |>
     x -> replace(x, r"(\b|\s)true(\b|\s)" =>" 1 = 1 ") |>
     x -> replace(x, r"(\b|\s)false(\b|\s)"=>" 1 = 0 ") |>
     x -> replace(x, r"^(IF|if|If)"=>"")

    if occursin(Regex("([><=|!+\\-^&]{1,4})"), commandout)
      return metaoperator(commandout, logicset, verbose = verbose) |>
        reportfeasible(command, verbose) |>
        (x -> (push!(x.commands, command); return x))
    end

    println("\t Warning! { $command } not interpreted!")
end

logicalparse(commands::Array{String,1}, logicset::LogicalCombo; I...) =
  logicalparse(commands, logicset=logicset, I...)

logicalparse(command::String, logicset::LogicalCombo; I...) =
  logicalparse(command, logicset=logicset, I...)


text:
# record.jl
# =========
#
# Representation of a record of a BCF file.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/GeneticVariation.jl/blob/master/LICENSE

mutable struct Record
    # data and filled range
    data::Vector{UInt8}
    filled::UnitRange{Int}
    # fields
    sharedlen::UInt32
    indivlen::UInt32
end

"""
    BCF.Record()
Create an unfilled BCF record.
"""
function Record()
    return Record(UInt8[], 1:0, 0, 0)
end

function isfilled(record::Record)
    return !isempty(record.filled)
end

function datarange(record::Record)
    return record.filled
end

function checkfilled(record::Record)
    if !isfilled(record)
        throw(ArgumentError("unfilled BCF record"))
    end
end

function Base.:(==)(record1::Record, record2::Record)
    if isfilled(record1) == isfilled(record2) == true
        r1 = datarange(record1)
        r2 = datarange(record2)
        return length(r1) == length(r2) && memcmp(pointer(record1.data, first(r1)), pointer(record2.data, first(r2)), length(r1)) == 0
    else
        return isfilled(record1) == isfilled(record2) == false
    end
end

function Base.copy(record::Record)
    return Record(record.data[record.filled], record.filled, record.sharedlen, record.indivlen)
end

function Record(base::Record;
                chrom=nothing, pos=nothing, qual=nothing,
                id=nothing, ref=nothing, alt=nothing,
                filter=nothing, info=nothing, genotype=nothing)
    checkfilled(base)
    data = base.data[1:24]

    if chrom != nothing
        store!(data, 0, convert(Int32, chrom - 1))
    end
    if pos != nothing
        store!(data, 4, convert(Int32, pos - 1))
    end
    if qual != nothing
        store!(data, 12, convert(Float32, qual))
    end

    offset = boffset = 24
    if id == nothing
        offset, boffset = copyvec!(data, offset, base.data, boffset)
    else
        offset = storestr!(data, offset, string(id))
        boffset = skipvec(base.data, boffset)
    end

    if ref == nothing
        offset, boffset = copyvec!(data, offset, base.data, boffset)
    else
        ref = string(ref)
        offset = storestr!(data, offset, ref)
        if sizeof(ref) != rlen(base)
            store!(data, 8, convert(Int32, sizeof(ref)))
        end
        boffset = skipvec(base.data, boffset)
    end

    n = n_allele(base)
    if alt == nothing
        for _ in 1:n-1
            offset, boffset = copyvec!(data, offset, base.data, boffset)
        end
    else
        if !isa(alt, Vector)
            alt = [alt]
        end
        for x in alt
            offset = storestr!(data, offset, string(x))
        end
        for _ in 1:n-1
            boffset = skipvec(base.data, boffset)
        end
        if length(alt) != n
            n_allele_info = load(UInt32, data, 16)[1]
            n_allele_info = (n_allele_info & 0x0000ffff) | (UInt32(length(alt) + 1) << 16)
            store!(data, 16, n_allele_info)
        end
    end

    if filter == nothing
        offset, boffset = copyvec!(data, offset, base.data, boffset)
    else
        if !isa(filter, Vector)
            filter = [filter]
        end
        offset = storevec!(data, offset, convert(Vector{Int8}, filter .- 1))
        boffset = skipvec(base.data, boffset)
    end

    n = n_info(base)
    if info == nothing
        for _ in 1:n
            # copy key and value(s)
            offset, boffset = copyvec!(data, offset, base.data, boffset)
            offset, boffset = copyvec!(data, offset, base.data, boffset)
        end
    else
        if !isa(info, AbstractDict)
            throw(ArgumentError("info must be an abstract dictionary"))
        end
        keyvec = Vector{Int8}(undef, 1)
        for (key, val) in info
            if !isa(key, Integer)
                throw(ArgumentError("info key must be an integer"))
            end
            keyvec[1] = key
            offset = storevec!(data, offset, keyvec)
            if isa(val, String)
                offset = storestr!(data, offset, val)
            elseif !isa(val, Vector)
                offset = storevec!(data, offset, [val])
            else
                offset = storevec!(data, offset, val)
            end
        end
        if length(info) != n
            if length(info) > typemax(UInt16)
                throw(ArgumentError("too many info fields"))
            end
            n_allele_info = load(UInt32, data, 16)[1]
            n_allele_info = (n_allele_info & 0xffff0000) | UInt32(length(info))
            store!(data, 16, n_allele_info)
        end
        for _ in 1:n
            # skip key and value(s)
            boffset = skipvec(base.data, boffset)
            boffset = skipvec(base.data, boffset)
        end
    end
    sharedlen = offset

    n = n_format(base)
    if genotype == nothing
        N = n_sample(base)
        for _ in 1:n
            # copy key
            offset, boffset = copyvec!(data, offset, base.data, boffset)
            # copy value(s)
            head, boffset = loadvechead(base.data, boffset)
            offset = storevechead!(data, offset, head)
            for _ in 1:N
                offset, boffset = copyvecbody!(data, offset, base.data, boffset, head)
            end
        end
    else
        error("modifying genotype is not yet supported")
    end

    return Record(data, 1:lastindex(data), sharedlen, offset - sharedlen)
end

function Base.show(io::IO, record::Record)
    print(io, summary(record), ':')
    if isfilled(record)
        println(io)
        println(io, "   chromosome: ", chrom(record))
        println(io, "     position: ", pos(record))
        println(io, "   identifier: ", id(record))
        println(io, "    reference: ", ref(record))
        println(io, "    alternate: ", join(alt(record), ' '))
        println(io, "      quality: ", qual(record))
        println(io, "       filter: ", join(filter(record), ' '))
          print(io, "  information: ", info(record))
    else
        print(io, " <not filled>")
    end
end

# Accessor functions
# ------------------

"""
    chrom(record::Record)::Int
Get the chromosome index of `record`.
"""
function chrom(record::Record)::Int
    checkfilled(record)
    return load(Int32, record.data, 0)[1] % Int + 1
end

"""
    pos(record::Record)::Int
Get the reference position of `record`.
Note that the position of the first base is 1 (i.e. 1-based coordinate).
"""
function pos(record::Record)::Int
    checkfilled(record)
    return load(Int32, record.data, 4)[1] % Int + 1
end

"""
    rlen(record::Record)::Int
Get the length of `record` projected onto the reference sequence.
"""
function rlen(record::Record)::Int
    checkfilled(record)
    return load(Int32, record.data, 8)[1] % Int
end

"""
    qual(record::Record)::Float32
Get the quality score of `record`.
Note that `0x7F800001` (signaling NaN) is interpreted as a missing value.
"""
function qual(record::Record)
    checkfilled(record)
    # 0x7F800001 is a missing value.
    return load(Float32, record.data, 12)[1]
end

function n_allele(rec::Record)
    checkfilled(rec)
    return (load(Int32, rec.data, 16)[1] >> 16) % Int
end

function n_info(rec::Record)
    checkfilled(rec)
    return (load(Int32, rec.data, 16)[1] & 0x0000FFFF) % Int
end

function n_format(rec::Record)
    checkfilled(rec)
    return (load(UInt32, rec.data, 20)[1] >> 24) % Int
end

function n_sample(rec::Record)
    checkfilled(rec)
    return (load(UInt32, rec.data, 20)[1] & 0x00FFFFFF) % Int
end

function id(rec::Record)
    checkfilled(rec)
    offset = 24
    return loadstr(rec.data, offset)[1]
end

"""
    ref(record::Record)::String
Get the reference bases of `record`.
"""
function ref(record::Record)::String
    checkfilled(record)
    # skip ID
    offset = 24
    len, offset = loadveclen(record.data, offset)
    # load REF
    return loadstr(record.data, offset + len)[1]
end

"""
    alt(record::Record)::Vector{String}
Get the alternate bases of `record`.
"""
function alt(record::Record)
    checkfilled(record)
    # skip ID and REF
    offset = 24
    len, offset = loadveclen(record.data, offset)
    len, offset = loadveclen(record.data, offset + len)
    # load ALTs
    N = n_allele(record) - 1
    alt = Vector{String}(undef, N)
    for n in 1:N
        str, offset = loadstr(record.data, offset + len)
        alt[n] = str
        len = sizeof(str)
    end
    return alt
end

"""
    filter(record::Record)::Vector{Int}
Get the filter indexes of `record`.
"""
function filter(record::Record)::Vector{Int}
    checkfilled(record)
    # skip ID, REF and ALTs
    offset = 24
    len = 0
    for _ in 1:n_allele(record)+1
        len, offset = loadveclen(record.data, offset + len)
    end
    # load FILTER
    return loadvec(record.data, offset + len)[1] .+ 1
end

"""
    info(record::Record, [simplify::Bool=true])::Vector{Tuple{Int,Any}}
Get the additional information of `record`.
When `simplify` is `true`, a vector with a single element is converted to the
element itself and an empty vector of the void type is converted to `nothing`.
"""
function info(record::Record; simplify::Bool=true)::Vector{Tuple{Int,Any}}
    checkfilled(record)
    # skip ID, REF, ALTs and FILTER
    offset::Int = 24
    len = 0
    for _ in 1:n_allele(record)+2
        len, offset = loadveclen(record.data, offset + len)
    end
    offset += len
    # load INFO
    ret = Vector{Tuple{Int,Any}}(undef, n_info(record))
    for i in 1:lastindex(ret)
        key, offset = loadvec(record.data, offset)
        @assert length(key) == 1
        val, offset = loadvec(record.data, offset)
        if simplify
            if isa(val, Vector) && length(val) == 1
                val = val[1]
            elseif isa(val, Vector{Nothing}) && isempty(val)
                val = nothing
            end
        end
        ret[i] = (key[1], val)
    end
    return ret
end

"""
    info(record::Record, key::Integer, [simplify::Bool=true])
Get the additional information of `record` with `key`.
When `simplify` is `true`, a vector with a single element is converted to the
element itself and an empty vector of the Nothing type is converted to `nothing`.
"""
function info(record::Record, key::Integer; simplify::Bool=true)
    checkfilled(record)
    # skip ID, REF, ALTs and FILTER
    offset::Int = 24
    len = 0
    for _ in 1:n_allele(record) + 2
        len, offset = loadveclen(record.data, offset + len)
    end
    offset += len
    # load INFO
    for _ in 1:n_info(record)
        k, offset = loadvec(record.data, offset)
        @assert length(k) == 1
        if k[1] == key
            val, offset = loadvec(record.data, offset)
            if simplify
                if isa(val, Vector) && length(val) == 1
                    val = val[1]
                elseif isa(val, Vector{Nothing}) && isempty(val)
                    val = nothing
                end
            end
            return val
        else
            offset = skipvec(record.data, offset)
        end
    end
    throw(KeyError(key))
end

"""
    genotype(record::Record)::Vector{Tuple{Int,Vector{Any}}}
Get the genotypes of `record`.
BCF genotypes are encoded by field, not by sample like VCF.
"""
function genotype(record::Record)::Vector{Tuple{Int,Vector{Any}}}
    checkfilled(record)
    offset::Int = record.sharedlen
    N = n_sample(record)
    ret = Tuple{Int,Vector{Any}}[]
    for j in 1:n_format(record)
        key, offset = loadvec(record.data, offset)
        @assert length(key) == 1
        head, offset = loadvechead(record.data, offset)
        vals = Vector{Any}[]
        for n in 1:N
            val, offset = loadvecbody(record.data, offset, head)
            push!(vals, val)
        end
        push!(ret, (key[1], vals))
    end
    return ret
end

function genotype(record::Record, index::Integer)
    return [(k, geno[index]) for (k, geno) in genotype(record)]
end

function genotype(record::Record, index::Integer, key::Integer)
    checkfilled(record)
    N = n_sample(record)
    offset::Int = record.sharedlen
    for j in 1:n_format(record)
        k, offset = loadvec(record.data, offset)
        @assert length(k) == 1
        head, offset = loadvechead(record.data, offset)
        for n in 1:N
            if k[1] == key && n == index
                return loadvecbody(record.data, offset, head)[1]
            else
                offset = skipvecbody(record.data, offset, head)
            end
        end
    end
    throw(KeyError(key))
end

function genotype(record::Record, indexes::AbstractVector{<:Integer}, key::Integer)
    checkfilled(record)
    N = n_sample(record)
    offset::Int = record.sharedlen
    for j in 1:n_format(record)
        k, offset = loadvec(record.data, offset)
        @assert length(k) == 1
        head, offset = loadvechead(record.data, offset)
        if k[1] == key
            vals = Vector{Any}[]
            for n in 1:N
                val, offset = loadvecbody(record.data, offset, head)
                push!(vals, val)
            end
            return vals[indexes]
        else
            for n in 1:N
                offset = skipvecbody(record.data, offset, head)
            end
        end
    end
    throw(KeyError(key))
end

function genotype(record::Record, ::Colon, key::Integer)
    return genotype(record, 1:n_sample(record), key)
end

function gt(x::Int8)
    allele = (x >> 1) - 1
    phased = x & 1 != 0
    return allele, phased
end

# Resolve type encoding.
function bcftype(b::UInt8)
    b &= 0x0f
    if b == 0x00
        # this is not defined in the BCF specs v2.2 but used in htslib
        return :void
    elseif b == 0x01
        return :int8
    elseif b == 0x02
        return :int16
    elseif b == 0x03
        return :int32
    elseif b == 0x05
        return :float32
    elseif b == 0x07
        return :character
    else
        return :unused
    end
end

function bcftypetag(typ::Symbol, len::Int)
    b = 0x00
    if typ == :void
        b |= 0x00
        @assert len == 0
    elseif typ == :int8
        b |= 0x01
    elseif typ == :int16
        b |= 0x02
    elseif typ == :int32
        b |= 0x03
    elseif typ == :float32
        b |= 0x05
    elseif typ == :character
        b |= 0x07
    else
        throw(ArgumentError("invalid type name: $(typ)"))
    end

    if len < 0
        throw(ArgumentError("length must be non-negative"))
    elseif len < 15
        b |= (len % UInt8) << 4
    else
        b |= (15 % UInt8)  << 4
    end
    return b
end

function bcftypesize(typ::Symbol)
    return typ == :int8      ? 1 :
           typ == :int16     ? 2 :
           typ == :int32     ? 4 :
           typ == :float32   ? 4 :
           typ == :character ? 1 :
           typ == :void      ? 0 :
           error("size unknown")
end

function load(::Type{T}, data::Vector{UInt8}, offset::Int) where T
    if offset + sizeof(T) > length(data)
        throw(BoundsError())
    end
    return unsafe_load(convert(Ptr{T}, pointer(data, offset + 1))), offset + sizeof(T)
end

# Load a string; this is basically same as `loadvec` but type-stable.
function loadstr(data::Vector{UInt8}, offset::Int)
    (typ, len), offset = loadvechead(data, offset)
    @assert typ == :character
    return String(data[offset+1:offset+len]), offset
end

function loadvec(data::Vector{UInt8}, offset::Int)
    head, offset = loadvechead(data, offset)
    return loadvecbody(data, offset, head)
end

function loadveclen(data::Vector{UInt8}, offset::Int)
    (_, len), offset = loadvechead(data, offset)
    return len, offset
end

function loadvechead(data::Vector{UInt8}, offset::Int)
    b = data[offset+1]
    typ = bcftype(b)
    len::Int = b >> 4
    if len < 15
        return (typ, len), offset + 1
    end
    lenvec, offset = loadvec(data, offset + 1)
    @assert length(lenvec) == 1
    return (typ, Int(lenvec[1])), offset
end

function loadvecbody(data::Vector{UInt8}, offset::Int, head::Tuple{Symbol,Int})
    t, len = head
    T = t == :int8 ? Int8 :
        t == :int16 ? Int16 :
        t == :int32 ? Int32 :
        t == :float32 ? Float32 :
        t == :character ? UInt8 : Nothing
    ret = Vector{T}(undef, len)
    for i in 1:len
        ret[i], offset = load(T, data, offset)
    end
    if t == :character
        return String(ret), offset
    else
        return ret, offset
    end
end

function skipvec(data::Vector{UInt8}, offset::Int)
    (t, len), offset = loadvechead(data, offset)
    return offset + bcftypesize(t) * len
end

function skipvecbody(data::Vector{UInt8}, offset::Int, head::Tuple{Symbol,Int})
    t, len = head
    return offset + bcftypesize(t) * len
end

function store!(data::Vector{UInt8}, offset::Int, val::T) where T
    if offset + sizeof(T) > length(data)
        throw(BoundsError())
    end
    unsafe_store!(convert(Ptr{T}, pointer(data, offset + 1)), val)
    return offset + sizeof(T)
end

function storestr!(data::Vector{UInt8}, offset::Int, str::String)
    if !isascii(str)
        throw(ArgumentError("string must be ASCII"))
    end
    len = sizeof(str)
    offset = storevechead!(data, offset, (:character, len))
    if length(data) < offset + len
        resize!(data, offset + len)
    end
    memcpy(pointer(data, offset + 1), pointer(str), len)
    return offset + len
end

function storevec!(data::Vector{UInt8}, offset::Int, vec::Vector{T}) where T
    t = T == Int8 ? :int8 :
        T == Int16 ? :int16 :
        T == Int32 ? :int32 :
        T == Float32 ? :float32 :
        error("unsupported type: $(T)")
    len = length(vec)
    offset = storevechead!(data, offset, (t, len))
    return storevecbody!(data, offset, vec)
end

function storevechead!(data::Vector{UInt8}, offset::Int, head::Tuple{Symbol,Int})
    t, len = head
    if length(data) < offset + 1
        resize!(data, offset + 1)
    end
    data[offset+1] = bcftypetag(t, len)
    if len < 15
        return offset + 1
    elseif len < typemax(Int8)
        return storevec!(data, offset + 1, Int8[len])
    elseif len < typemax(Int16)
        return storevec!(data, offset + 1, Int16[len])
    elseif len < typemax(Int32)
        return storevec!(data, offset + 1, Int32[len])
    else
        error("too long vector")
    end
end

function storevecbody!(data::Vector{UInt8}, offset::Int, vec::Vector{T}) where T
    len = length(vec)
    n = len * sizeof(T)
    if length(data) < offset + n
        resize!(data, offset + n)
    end
    memcpy(pointer(data, offset + 1), pointer(vec), n)
    return offset + n
end

function copyvec!(dst::Vector{UInt8}, doffset::Int, src::Vector{UInt8}, soffset::Int)
    head, soffset = loadvechead(src, soffset)
    doffset = storevechead!(dst, doffset, head)
    return copyvecbody!(dst, doffset, src, soffset, head)
end

function copyvecbody!(dst::Vector{UInt8}, doffset::Int, src::Vector{UInt8}, soffset::Int, head::Tuple{Symbol,Int})
    t, len = head
    n = bcftypesize(t) * len
    if length(dst) < doffset + n
        resize!(dst, doffset + n)
    end
    copyto!(dst, doffset + 1, src, soffset + 1, n)
    return doffset + n, soffset + n
end

function memcpy(dst::Ptr, src::Ptr, n::Int)
    ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t), dst, src, n)
    return nothing
end

function memcmp(p1::Ptr, p2::Ptr, n::Integer)
    return ccall(:memcmp, Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t), p1, p2, n)
end


text:
using JuMP

m = Model()

type Recipe
	name::AbstractString
	#nfacs::Int64
	time::Float64
	outs::Dict{AbstractString, Float64}
	ins::Dict{AbstractString, Float64}
	#Recipe(n, t) = new(n, 0, t, Dict{AbstractString, Float64}(), Dict{AbstractString, Float64}())
	Recipe(name, time, speed) = new(name, time/speed, Dict{AbstractString, Float64}(), Dict{AbstractString, Float64}())
end

function production(r::Recipe, nfacs)
	@printf "%s: %0.1f\n" r.name nfacs
	@printf "\tProduces\n"
	for (k,v) in r.outs
		@printf "\t\t%0.2f %s per min / %0.2f per sec\n" 60v * nfacs k v*nfacs
	end
	@printf "\tConsumes\n"
	for (k,v) in r.ins
		@printf "\t\t%0.2f %s per min / %0.2f per sec\n" 60v * nfacs k v*nfacs
	end
end

@enum Rnum PJ AOP HOC LOC PGSOL HOSOL LOSOL ROK

recips = Dict{Rnum, Recipe}()

macro IN(rnum, iname, iqty)
	:(recips[$rnum].ins[$iname] = $iqty / recips[$rnum].time)
end

macro OUT(rnum, iname, iqty)
	:(recips[$rnum].outs[$iname] = $iqty / recips[$rnum].time)
end

recips[PJ] = Recipe("Pumpjacks", 1.0, 1)
@IN PJ "Pump" 1
@OUT PJ "Crude Oil" 0.2

recips[AOP] = Recipe("Advanced Oil Processing", 5.0, 1.)
@IN AOP "Crude Oil" 200.
@IN AOP "Water" 100.
@OUT AOP "Heavy Oil" 10.
@OUT AOP "Light Oil" 45.
@OUT AOP "Petroleum Gas" 55.

recips[HOC] = Recipe("Heavy Oil Cracking", 3.0 , 1.25)
@IN HOC "Heavy Oil" 40.
@IN HOC "Water" 30.
@OUT HOC "Light Oil" 30.

recips[LOC] = Recipe("Light Oil Cracking", 5.0 , 1.25)
@IN LOC "Light Oil" 30.
@IN LOC "Water" 30.
@OUT LOC "Petroleum Gas" 20.

recips[PGSOL] = Recipe("PG Solid Fuel", 3., 1.25)
@IN PGSOL "Petroleum Gas" 20.
@OUT PGSOL "Solid Fuel" 1

recips[LOSOL] = Recipe("LO Solid Fuel", 3., 1.25)
@IN LOSOL "Light Oil" 10.
@OUT LOSOL "Solid Fuel" 1

recips[HOSOL] = Recipe("HO Solid Fuel", 3., 1.25)
@IN HOSOL "Heavy Oil" 20.0
@OUT HOSOL "Solid Fuel" 1

recips[ROK] = Recipe("Rocket Fuel Assembly", 30., 1.25)
@IN ROK "Solid Fuel" 10.
@OUT ROK "Rocket Fuel" 1

@variable(m, PJs >= 0) #, Int)
@variable(m, AOPs >= 0) #, Int)
@variable(m, HOCs >= 0) #, Int)
@variable(m, LOCs >= 0) #, Int)
@variable(m, PGSOLs >= 0) #, Int)
@variable(m, HOSOLs >= 0) #, Int)
@variable(m, LOSOLs >= 0) #, Int)
@variable(m, ROKs >= 0) #, Int)

macro FIN(Fn, F, R) :($Fn * recips[$F].ins[$R]) end
macro FOUT(Fn, F, R) :($Fn * recips[$F].outs[$R] )end

@constraint(m, AOPs == 10)

@constraint(m, @FOUT(PJs, PJ, "Crude Oil") <= 2400)
@constraint(m, @FIN(AOPs, AOP, "Crude Oil") == @FOUT(PJs, PJ, "Crude Oil"))
@constraint(m, @FIN(HOCs, HOC, "Heavy Oil") + @FIN(HOSOLs, HOSOL, "Heavy Oil") == @FOUT(AOPs, AOP, "Heavy Oil"))
@constraint(m, @FIN(LOCs, LOC, "Light Oil") + @FIN(LOSOLs, LOSOL, "Light Oil") == @FOUT(AOPs, AOP, "Light Oil") + @FOUT(HOCs, HOC, "Light Oil"))
@constraint(m, @FIN(PGSOLs, PGSOL, "Petroleum Gas") == @FOUT(AOPs, AOP, "Petroleum Gas") + @FOUT(LOCs, LOC, "Petroleum Gas"))
@constraint(m, @FIN(ROKs, ROK, "Solid Fuel") == @FOUT(PGSOLs, PGSOL, "Solid Fuel") + @FOUT(HOSOLs, HOSOL, "Solid Fuel") + @FOUT(LOSOLs, LOSOL, "Solid Fuel"))

@objective(m, Max, sum([PGSOLs, HOSOLs, LOSOLs]))
solve(m)

production(recips[PJ], getvalue(PJs))
production(recips[AOP], getvalue(AOPs))
production(recips[HOC], getvalue(HOCs))
production(recips[LOC], getvalue(LOCs))
production(recips[PGSOL], getvalue(PGSOLs))
production(recips[LOSOL], getvalue(LOSOLs))
production(recips[HOSOL], getvalue(HOSOLs))
production(recips[ROK], getvalue(ROKs))

max_stars_repo_name: lawless-m/factorio-utils

max_issues_repo_name: lawless-m/factorio-utils

max_forks_repo_name: lawless-m/factorio-utils


text:
using XRayQuadrics
using LinearAlgebra
using Colors
using Plots

n = [0;0;1]
c = [1;1;1]

s = Plane(c,n)
qs = Quadric(s)

c = Cylinder(3, [1;1;1], [1;0;0])
qc = Quadric(c)

p = Paraboloid(3, [1;1;1], [0;0;1])
qp = Quadric(p)

h = Hyperboloid(3,8, [2;2;2], [0;1;0])
qh = Quadric(h)

ray = Particle([0;0;10],[0;0;-1],10e3)

tp = in_out(ray, qp)
tc = in_out(ray, qc)
th = in_out(ray, qh)

Ec = eigvals(qc.Q[1:3,1:3])
Eh = eigvals(qh.Q[1:3,1:3])
Ep = eigvals(qp.Q[1:3,1:3])
vc = eigvecs(qc.Q[1:3,1:3])
vh = eigvecs(qh.Q[1:3,1:3])
vp = eigvecs(qp.Q[1:3,1:3])

classify(qc)

ax = [1;1;1]
ax = ax/norm(ax)
center = [0;0;0]
R = 1
b = 1
c = Cylinder(R, center, ax)
p = Paraboloid(b, center, ax)
h = Hyperboloid(R, b, center, ax)

p1 = Plane(c.a*0.01, ax)
p2 = Plane(c.a*6, ax)

qc = Quadric(c)
qp = Quadric(p)
qh = Quadric(h)

tqc = TruncatedQuadric(c, [p1, p2], true)
tqp = TruncatedQuadric(p, [p1, p2], true)
tqh = TruncatedQuadric(h, [p1, p2], true)

# f = Figure()
# scatter([0, 1],[0, 1],[0, 1], plot_color=:red)
# Makie.plot(tq)
# Makie.plot!(tq)
# s = GLMakie.plot(tqp)

plotlyjs(reuse=true)
splot = Plots.surface(tqp, color=:blue, plot_caps=true, size=(1920,1080))
display(splot)
Plots.surface!(tqh, color=:red)
Plots.surface!(tqc, color=:green)
Plots.plot!([center center+3*ax][1,:], [center center+3*ax][2,:], [center center+3*ax][3,:], color=:blue, linewidth=4)
# current()

nray = 1000
rays = Array{Particle, 1}(undef, nray)

chits = zeros(1,nray)
phits = zeros(1,nray)
hhits = zeros(1,nray)

crays = zeros(3,nray)
prays = zeros(3,nray)
hrays = zeros(3,nray)

cnorm = zeros(3,nray)
pnorm = zeros(3,nray)
hnorm = zeros(3,nray)

for i = 1:nray
    dir = rand(3,) - 0.5*ones(3,)
    dir = dir/norm(dir)
    rays[i] = Particle(center + 1*ax, dir, 10)

    (cin, cout) = in_out(rays[i], qc)
    (pin, pout) = in_out(rays[i], qp)
    (hin, hout) = in_out(rays[i], qh)
    chits[i] = max(cin, cout)
    phits[i] = max(pin, pout)
    hhits[i] = max(hin, hout)
    crays[:,i] = rays[i].r0 + chits[i]*rays[i].v
    prays[:,i] = rays[i].r0 + phits[i]*rays[i].v
    hrays[:,i] = rays[i].r0 + hhits[i]*rays[i].v
    cnorm[:,i] = -normal(qc, crays[:,i])
    pnorm[:,i] = normal(qp, prays[:,i])
    hnorm[:,i] = normal(qh, hrays[:,i])
end

# Plots.quiver(prays[1,:], prays[2,:], prays[3,:], quiver=(pnorm[1,:], pnorm[2,:], pnorm[3,:]))

Plots.scatter!([prays[1,:] prays[1,:]+pnorm[1,:]], [prays[2,:] prays[2,:]+pnorm[2,:]], [prays[3,:] prays[3,:]+pnorm[3,:]], color=:blue, markersize=1)

Plots.scatter!([hrays[1,:] hrays[1,:]+hnorm[1,:]], [hrays[2,:] hrays[2,:]+hnorm[2,:]], [hrays[3,:] hrays[3,:]+hnorm[3,:]], color=:red, markersize=1)

Plots.scatter!([crays[1,:] crays[1,:]+cnorm[1,:]], [crays[2,:] crays[2,:]+cnorm[2,:]], [crays[3,:] crays[3,:]+cnorm[3,:]], color=:green, markersize=1)

a = 9
Plots.xlims!((center[1]-a,center[1]+a))
Plots.ylims!((center[2]-a,center[2]+a))
Plots.zlims!((center[3]-a,center[3]+a))

display(splot)


text:
# Hello world example, similar to the Boost.Python hello world

using CxxWrap
using Base.Test
using Compat

# Wrap the functions defined in C++
wrap_modules(CxxWrap._l_parametric)

import ParametricTypes.TemplateType, ParametricTypes.NonTypeParam

p1 = TemplateType{ParametricTypes.P1, ParametricTypes.P2}()
p2 = TemplateType{ParametricTypes.P2, ParametricTypes.P1}()

println("Dumping object p1:")
dump(p1)

@test ParametricTypes.get_first(p1) == 1
@test ParametricTypes.get_second(p2) == 1
@test typeof(ParametricTypes.get_first(p1)) == Int32
@test typeof(ParametricTypes.get_second(p2)) == Int32

@test ParametricTypes.get_first(p2) == 10.
@test ParametricTypes.get_second(p1) == 10.
@test typeof(ParametricTypes.get_first(p2)) == Float64
@test typeof(ParametricTypes.get_second(p1)) == Float64

@test ParametricTypes.TemplateDefaultType{ParametricTypes.P1}() != nothing

nontype1 = ParametricTypes.NonTypeParam{Int32, Int32(1)}()
@test ParametricTypes.get_nontype(nontype1) == 1

nontype2 = ParametricTypes.NonTypeParam{UInt32, UInt32(2)}()
@test ParametricTypes.get_nontype(nontype2) == UInt32(2)

nontype3 = ParametricTypes.NonTypeParam{Int32, Int32(1)}(3)
@test ParametricTypes.get_nontype(nontype3) == 3

nontype4 = ParametricTypes.NonTypeParam{Int64, Int64(64)}()
@test ParametricTypes.get_nontype(nontype4) == Int64(64)

concr = ParametricTypes.ConcreteTemplate{Float64}()
@test isa(concr, ParametricTypes.AbstractTemplate{Float64})
@test isa(concr, ParametricTypes.AbstractTemplate)
@test isa(concr, ParametricTypes.ConcreteTemplate)
abst = ParametricTypes.to_base(concr)
@test isa(abst, ParametricTypes.AbstractTemplate{Float64})
@test isa(abst, ParametricTypes.AbstractTemplate)

f3 = ParametricTypes.Foo3{Int32, Bool, Float32}()
@test length(methods(ParametricTypes.foo3_method)) == 6
f2 = ParametricTypes.Foo2{Float64}()
@test length(methods(ParametricTypes.foo2_method)) == 2

@test length(methods(ParametricTypes.foo3_free_method)) == 6
ParametricTypes.foo3_free_method(f3)

@test supertype(ParametricTypes.Foo3{Float64,ParametricTypes.P1,Float32}) == ParametricTypes.AbstractTemplate{Float64}

vec1 = ParametricTypes.CppVector{Float64}()
vec2 = ParametricTypes.CppVector2{Float64, Float32}()
@test isa(vec1, AbstractVector{Float64})
@test isa(vec2, AbstractVector{Float64})


text:
using Localize

@i18n Base.Math begin

zh"""
    sin(x)

正弦函数
"""
sin(::Number)

# uncomment this when utils for gettext is done
# @lang zh_TW """
#     sin(x)
#
# 正弦函數
# """
# sin(::Number)
#
# fr"""
#     sin(x)
#
# fonction sinusoïdale
# """
# sin(::Number)
#
# ja"""
#     sin(x)
#
# 正弦関数
# """
# sin(::Number)

end

# This fails because of the IdDict that current doc system used
# using LinearAlgebra
#
# @i18n LinearAlgebra begin
#
# zh"""
#     sin(A::AbstractMatrix)
#
# 计算方阵 `A` 的正弦函数。
# """
# sin(::AbstractMatrix{<:Real})
#
# ja"""
#     sin(A::AbstractMatrix)
#
# 行列 `A` の正弦関数
# """
# sin(::AbstractMatrix{<:Real})
#
# end


text:
# ----------
# MHD Analysis Method, providing MHD related quantities function
# ----------

# Scale Decomposition FUnction
function ScaleDecomposition(B1::Array,B2::Array,B3::Array;kf=[1,5],Lx = 2π)
    k1,k2 = minimum(kf),maximum(kf);
    nx,ny,nz = size(B1);
    T    = Float32;
    grid = ThreeDGrid(nx, Lx, T = T);

    B1h = zeros(ComplexF32,(div(nx,2)+1,ny,nz));
    B2h = zeros(ComplexF32,(div(nx,2)+1,ny,nz));
    B3h = zeros(ComplexF32,(div(nx,2)+1,ny,nz));
    Bxhf = copy(B1h);
    Byhf = copy(B2h);
    Bzhf = copy(B3h);
    mul!(B1h, grid.rfftplan, B1);
    mul!(B2h, grid.rfftplan, B2);
    mul!(B3h, grid.rfftplan, B3);

    for i in 1:div(nz,2)+1, j in 1:ny, k in 1:nx
       x,y,z = grid.kr[i],grid.l[j],grid.m[k];
       rr    = sqrt(x^2+y^2+z^2);
       if (( rr >=  k1) && (rr <= k2))
           Bxhf[i,j,k] = B1h[i,j,k];
           Byhf[i,j,k] = B2h[i,j,k];
           Bzhf[i,j,k] = B3h[i,j,k];
       end
    end
    cB1,cB2,cB3 = zeros(T,size(B1)),zeros(T,size(B1)),zeros(T,size(B1));
    ldiv!(cB1, grid.rfftplan,Bxhf);
    ldiv!(cB2, grid.rfftplan,Byhf);
    ldiv!(cB3, grid.rfftplan,Bzhf);
    return cB1,cB2,cB3;
end

#Kenitic Helicity
function h_k(iv,jv,kv;L=2π)
# V ⋅ ( ∇ × V )
	dlx,dly,dlz  = (L/size(iv)[1]),(L/size(iv)[2]),(L/size(iv)[3]);
	dV = dlx*dly*dlz;
	cV1,cV2,cV3 = Curl(iv,jv,kv;Lx=L);
	h_k = @. (cV1*iv + cV2*jv + cV3*kv)*dV
	return h_k
end

#Megnetic Helicity
function h_m(ib,jb,kb)
# A ⋅ B
   Ax,Ay,Az = VectorPotential(ib,jb,kb);
   return Ax.*ib + Ay.*jb + Az.*kb;
end

#Vector Potential
function VectorPotential(B1,B2,B3;L=2π)
#=
    funtion of computing B = ∇ × A using the fourier method
     fft(∇×Vector) -> im * k × A
      | i j k  |
      | x y z  |  =  (y*A3 - z*A2) i - (x*A3 - z*A1) j + (x*A2 - y*A1) k
      |A1 A2 A3|

    Note: We are using the Coulomb gauge ∇ ⋅ A = 0
    Using the relations J = ∇ × (∇ × B) and ∇ ⋅ A = 0 in peroideric condition,
    we will arrive Jₖ = - k² Aₖ.

    Aₖ = (k × B)ᵢ/ k²
=#
    nx,ny,nz = size(B1);
    T    = Float32;
    grid = ThreeDGrid(nx, L, T = T);

    B1h = zeros(ComplexF32,(div(nx,2)+1,ny,nz));
    B2h = zeros(ComplexF32,(div(nx,2)+1,ny,nz));
    B3h = zeros(ComplexF32,(div(nx,2)+1,ny,nz));
    Axh = copy(B1h);
    Ayh = copy(B2h);
    Azh = copy(B3h);
    mul!(B1h, grid.rfftplan, B1);
    mul!(B2h, grid.rfftplan, B2);
    mul!(B3h, grid.rfftplan, B3);

    for i in 1:div(nz,2)+1, j in 1:ny, k in 1:nx
       x,y,z = grid.kr[i],grid.l[j],grid.m[k];
       k²    = x^2 + y^2 + z^2;
       Axh[i,j,k] = im*(y*B3h[i,j,k] - z*B2h[i,j,k])/k²;
       Ayh[i,j,k] = im*(z*B1h[i,j,k] - x*B3h[i,j,k])/k²;
       Azh[i,j,k] = im*(x*B2h[i,j,k] - y*B1h[i,j,k])/k²;
       Axh[i,j,k] = ifelse(k² == 0, 0,  Axh[i,j,k]);
       Ayh[i,j,k] = ifelse(k² == 0, 0,  Ayh[i,j,k]);
       Azh[i,j,k] = ifelse(k² == 0, 0,  Azh[i,j,k]);
    end

    A1,A2,A3 = zeros(T,size(B1)),zeros(T,size(B1)),zeros(T,size(B1));
    ldiv!(A1, grid.rfftplan, deepcopy(Axh));
    ldiv!(A2, grid.rfftplan, deepcopy(Ayh));
    ldiv!(A3, grid.rfftplan, deepcopy(Azh));
    return A1,A2,A3;

end

#Checking for anagular momentum using center point as a reference (r = 0 at the center)
function getL(iv,jv,kv;L=2π)
    nx,ny,nz = size(iv);
    T    = Float32;
    grid = ThreeDGrid(nx, L, T = T);
    Lᵢ,Lⱼ,Lₖ  = getL(iv,jv,kv,grid)
    return Lᵢ,Lⱼ,Lₖ
end

function getL(iv,jv,kv,grid)
    nx,ny,nz = size(iv);
    T    = Float32;
    Lᵢ,Lⱼ,Lₖ = zeros(nx,ny,nz),zeros(nx,ny,nz),zeros(nx,ny,nz);
    for k in 1:nz, j in 1:ny, i in 1:nx
       x,y,z = grid.x[i],grid.y[j],grid.z[k];
       r = sqrt(x^2 + y^2);
       # L = r × p => iv,jv,kv
       # |  i  j  k |
       # |  x  y  z | = (y*kv - z*jv) i - (x*kv - z*iv) j + (x*jv - y*iv) k;
       # | iv jv kv |
       Lᵢ[i,j,k] =   (y*kv[i,j,k] - z*jv[i,j,k]);
       Lⱼ[i,j,k] =  -(x*kv[i,j,k] - z*iv[i,j,k]);
       Lₖ[i,j,k] =   (x*jv[i,j,k] - y*iv[i,j,k]);
    end
    return Lᵢ,Lⱼ,Lₖ
end

max_stars_repo_name: MHDFlows/MHDFlows

max_issues_repo_name: MHDFlows/MHDFlows

max_forks_repo_name: MHDFlows/MHDFlows


text:
let
    # roundtrip test
    # TODO this could be done with Quickcheck if we had a way of
    # generating "interesting" HTML documents
    doc = open("$testdir/fixtures/example.html") do example
        read(example, String) |> parsehtml
    end
    io = IOBuffer()
    print(io, doc)
    seek(io, 0)
    newdoc = read(io, String) |> parsehtml
    @test newdoc == doc
end

tests = [
    "30",  # regression test for issue #30
    "multitext",  # regression test for multiple HTMLText in one HTMLElement
    "varied",  # relatively complex example
]
for test in tests
    let
        doc = open("$testdir/fixtures/$(test)_input.html") do example
            read(example, String) |> parsehtml
        end
        io = IOBuffer()
        print(io, doc.root, pretty=true)
        seek(io, 0)
        @test read(io, String) ==
            read(open("$testdir/fixtures/$(test)_output.html"), String)
    end
end


text:
export PhasePartition
# Thermodynamic states
export ThermodynamicState,
    PhaseDry,
    PhaseDry_given_pT,
    PhaseEquil,
    PhaseNonEquil,
    TemperatureSHumEquil,
    LiquidIcePotTempSHumEquil,
    LiquidIcePotTempSHumEquil_given_pressure,
    LiquidIcePotTempSHumNonEquil,
    LiquidIcePotTempSHumNonEquil_given_pressure

"""
    PhasePartition

Represents the mass fractions of the moist air mixture.

# Constructors

    PhasePartition(q_tot::Real[, q_liq::Real[, q_ice::Real]])
    PhasePartition(ts::ThermodynamicState)

See also [`PhasePartition_equil`](@ref)

# Fields

$(DocStringExtensions.FIELDS)
"""
struct PhasePartition{FT <: Real}
    "total specific humidity"
    tot::FT
    "liquid water specific humidity (default: `0`)"
    liq::FT
    "ice specific humidity (default: `0`)"
    ice::FT
end

PhasePartition(q_tot::FT, q_liq::FT) where {FT <: Real} =
    PhasePartition(q_tot, q_liq, zero(FT))
PhasePartition(q_tot::FT) where {FT <: Real} =
    PhasePartition(q_tot, zero(FT), zero(FT))

"""
    ThermodynamicState{FT}

A thermodynamic state, which can be initialized for
various thermodynamic formulations (via its sub-types).
All `ThermodynamicState`'s have access to functions to
compute all other thermodynamic properties.
"""
abstract type ThermodynamicState{FT} end

"""
    PhaseEquil{FT} <: ThermodynamicState

A thermodynamic state assuming thermodynamic equilibrium (therefore, saturation adjustment
may be needed).

# Constructors

    PhaseEquil(e_int, ρ, q_tot)

# Fields

$(DocStringExtensions.FIELDS)
"""
struct PhaseEquil{FT, PS} <: ThermodynamicState{FT}
    "parameter set (e.g., planet parameters)"
    param_set::PS
    "internal energy"
    e_int::FT
    "density of air (potentially moist)"
    ρ::FT
    "total specific humidity"
    q_tot::FT
    "temperature: computed via [`saturation_adjustment`](@ref)"
    T::FT
end
function PhaseEquil(
    e_int::FT,
    ρ::FT,
    q_tot::FT,
    maxiter::Int = 3,
    tol::FT = FT(1e-1),
    sat_adjust::Function = saturation_adjustment,
    param_set::PS = MTPS(),
) where {FT <: Real, PS}
    # TODO: Remove these safety nets, or at least add warnings
    # waiting on fix: github.com/vchuravy/GPUifyLoops.jl/issues/104
    q_tot_safe = clamp(q_tot, FT(0), FT(1))
    T = sat_adjust(e_int, ρ, q_tot_safe, maxiter, tol, param_set)
    return PhaseEquil{FT, PS}(param_set, e_int, ρ, q_tot_safe, T)
end
function PhaseEquil(
    e_int::FT,
    ρ::FT,
    q_tot::FT,
    sat_adjust::Function,
) where {FT <: Real, PS}
    return PhaseEquil(e_int, ρ, q_tot, 3, FT(1e-1), sat_adjust, MTPS())
end
function PhaseEquil(
    e_int::FT,
    ρ::FT,
    q_tot::FT,
    maxiter::Int,
    tol::FT,
    param_set::PS,
) where {FT <: Real, PS}
    return PhaseEquil(
        e_int,
        ρ,
        q_tot,
        maxiter,
        tol,
        saturation_adjustment,
        MTPS(),
    )
end

"""
    PhaseDry{FT} <: ThermodynamicState

A dry thermodynamic state (`q_tot = 0`).

# Constructors

    PhaseDry(e_int, ρ)

# Fields

$(DocStringExtensions.FIELDS)
"""
struct PhaseDry{FT, PS} <: ThermodynamicState{FT}
    "parameter set (e.g., planet parameters)"
    param_set::PS
    "internal energy"
    e_int::FT
    "density of dry air"
    ρ::FT
end
PhaseDry(e_int::FT, ρ::FT, param_set::PS = MTPS()) where {FT, PS} =
    PhaseDry{FT, PS}(param_set, e_int, ρ)

"""
    PhaseDry_given_pT(p, T)

Constructs a [`PhaseDry`](@ref) thermodynamic state from:

 - `p` pressure
 - `T` temperature
"""
function PhaseDry_given_pT(
    p::FT,
    T::FT,
    param_set::PS = MTPS(),
) where {FT <: Real, PS}
    e_int = internal_energy(T, param_set)
    ρ = air_density(T, p, param_set)
    return PhaseDry{FT, PS}(param_set, e_int, ρ)
end

"""
    LiquidIcePotTempSHumEquil(θ_liq_ice, ρ, q_tot)

Constructs a [`PhaseEquil`](@ref) thermodynamic state from:

 - `θ_liq_ice` liquid-ice potential temperature
 - `ρ` (moist-)air density
 - `q_tot` total specific humidity
 - `tol` tolerance for saturation adjustment
 - `maxiter` maximum iterations for saturation adjustment
"""
function LiquidIcePotTempSHumEquil(
    θ_liq_ice::FT,
    ρ::FT,
    q_tot::FT,
    maxiter::Int = 30,
    tol::FT = FT(1e-1),
    param_set::PS = MTPS(),
) where {FT <: Real, PS}
    T = saturation_adjustment_q_tot_θ_liq_ice(
        θ_liq_ice,
        ρ,
        q_tot,
        maxiter,
        tol,
        param_set,
    )
    q_pt = PhasePartition_equil(T, ρ, q_tot, param_set)
    e_int = internal_energy(T, q_pt, param_set)
    return PhaseEquil{FT, PS}(param_set, e_int, ρ, q_tot, T)
end
LiquidIcePotTempSHumEquil(
    θ_liq_ice::FT,
    ρ::FT,
    q_tot::FT,
    param_set::PS,
) where {FT <: Real, PS} =
    LiquidIcePotTempSHumEquil(θ_liq_ice, ρ, q_tot, 30, FT(1e-1), param_set)

"""
    LiquidIcePotTempSHumEquil_given_pressure(θ_liq_ice, p, q_tot)

Constructs a [`PhaseEquil`](@ref) thermodynamic state from:

 - `θ_liq_ice` liquid-ice potential temperature
 - `p` pressure
 - `q_tot` total specific humidity
 - `tol` tolerance for saturation adjustment
 - `maxiter` maximum iterations for saturation adjustment
"""
function LiquidIcePotTempSHumEquil_given_pressure(
    θ_liq_ice::FT,
    p::FT,
    q_tot::FT,
    maxiter::Int = 30,
    tol::FT = FT(1e-1),
    param_set::PS = MTPS(),
) where {FT <: Real, PS}
    T = saturation_adjustment_q_tot_θ_liq_ice_given_pressure(
        θ_liq_ice,
        p,
        q_tot,
        maxiter,
        tol,
        param_set,
    )
    ρ = air_density(T, p, PhasePartition(q_tot), param_set)
    q = PhasePartition_equil(T, ρ, q_tot, param_set)
    e_int = internal_energy(T, q, param_set)
    return PhaseEquil{FT, PS}(param_set, e_int, ρ, q.tot, T)
end
LiquidIcePotTempSHumEquil_given_pressure(
    θ_liq_ice::FT,
    p::FT,
    q_tot::FT,
    param_set::PS,
) where {FT <: Real, PS} = LiquidIcePotTempSHumEquil_given_pressure(
    θ_liq_ice,
    p,
    q_tot,
    30,
    FT(1e-1),
    param_set,
)

"""
    TemperatureSHumEquil(T, p, q_tot)

Constructs a [`PhaseEquil`](@ref) thermodynamic state from temperature.

 - `T` temperature
 - `p` pressure
 - `q_tot` total specific humidity
"""
function TemperatureSHumEquil(
    T::FT,
    p::FT,
    q_tot::FT,
    param_set::PS = MTPS(),
) where {FT <: Real, PS}
    ρ = air_density(T, p, PhasePartition(q_tot), param_set)
    q = PhasePartition_equil(T, ρ, q_tot, param_set)
    e_int = internal_energy(T, q, param_set)
    return PhaseEquil{FT, PS}(param_set, e_int, ρ, q_tot, T)
end

"""
   	 PhaseNonEquil{FT} <: ThermodynamicState

A thermodynamic state assuming thermodynamic non-equilibrium (therefore, temperature can
be computed directly).

# Constructors

    PhaseNonEquil(e_int, q::PhasePartition, ρ)

# Fields

$(DocStringExtensions.FIELDS)

"""
struct PhaseNonEquil{FT, PS} <: ThermodynamicState{FT}
    "parameter set (e.g., planet parameters)"
    param_set::PS
    "internal energy"
    e_int::FT
    "density of air (potentially moist)"
    ρ::FT
    "phase partition"
    q::PhasePartition{FT}
end
function PhaseNonEquil(
    e_int::FT,
    ρ::FT,
    q::PhasePartition{FT} = q_pt_0(FT),
    param_set::PS = MTPS(),
) where {FT, PS}
    return PhaseNonEquil{FT, PS}(param_set, e_int, ρ, q)
end

"""
    LiquidIcePotTempSHumNonEquil(θ_liq_ice, ρ, q_pt)

Constructs a [`PhaseNonEquil`](@ref) thermodynamic state from:

 - `θ_liq_ice` liquid-ice potential temperature
 - `ρ` (moist-)air density
 - `q_pt` phase partition
and, optionally
 - `tol` tolerance for non-linear equation solve
 - `maxiter` maximum iterations for non-linear equation solve
"""
function LiquidIcePotTempSHumNonEquil(
    θ_liq_ice::FT,
    ρ::FT,
    q_pt::PhasePartition{FT},
    maxiter::Int = 5,
    tol::FT = FT(1e-1),
    param_set::PS = MTPS(),
) where {FT <: Real, PS}
    T = air_temperature_from_liquid_ice_pottemp_non_linear(
        θ_liq_ice,
        ρ,
        maxiter,
        tol,
        q_pt,
        param_set,
    )
    e_int = internal_energy(T, q_pt, param_set)
    return PhaseNonEquil{FT, PS}(param_set, e_int, ρ, q_pt)
end

"""
    LiquidIcePotTempSHumNonEquil_given_pressure(θ_liq_ice, p, q_pt)

Constructs a [`PhaseNonEquil`](@ref) thermodynamic state from:

 - `θ_liq_ice` liquid-ice potential temperature
 - `p` pressure
 - `q_pt` phase partition
"""
function LiquidIcePotTempSHumNonEquil_given_pressure(
    θ_liq_ice::FT,
    p::FT,
    q_pt::PhasePartition{FT},
    param_set::PS = MTPS(),
) where {FT <: Real, PS}
    T = air_temperature_from_liquid_ice_pottemp_given_pressure(
        θ_liq_ice,
        p,
        q_pt,
        param_set,
    )
    ρ = air_density(T, p, q_pt, param_set)
    e_int = internal_energy(T, q_pt, param_set)
    return PhaseNonEquil{FT, PS}(param_set, e_int, ρ, q_pt)
end

"""
    fixed_lapse_rate_ref_state(z::FT,
                               T_surface::FT,
                               T_min::FT) where {FT<:AbstractFloat}

Fixed lapse rate hydrostatic reference state
"""
function fixed_lapse_rate_ref_state(
    z::FT,
    T_surface::FT,
    T_min::FT,
    param_set::PS = MTPS(),
) where {FT <: AbstractFloat, PS}
    Γ = FT(grav) / FT(cp_d)
    z_tropopause = (T_surface - T_min) / Γ
    H_min = FT(R_d) * T_min / FT(grav)
    T = max(T_surface - Γ * z, T_min)
    p = FT(MSLP) * (T / T_surface)^(FT(grav) / (FT(R_d) * Γ))
    T == T_min && (p = p * exp(-(z - z_tropopause) / FT(H_min)))
    ρ = p / (FT(R_d) * T)
    return T, p, ρ
end

"""
    tested_convergence_range(FT, n)

A range of input arguments to thermodynamic state constructors
 - `e_int` internal energy
 - `ρ` (moist-)air density
 - `q_tot` total specific humidity
 - `q_pt` phase partition
 - `T` air temperature
 - `θ_liq_ice` liquid-ice potential temperature
that are tested for convergence in saturation adjustment.

Note that the output vectors are of size ``n*n_RH``, and they
should span the input arguments to all of the constructors.
"""
function tested_convergence_range(
    n::Int,
    ::Type{FT},
    param_set::PS = MTPS(),
) where {FT, PS}
    n_RS1 = 10
    n_RS2 = 20
    n_RS = n_RS1 + n_RS2
    z_range = range(FT(0), stop = FT(2.5e4), length = n)
    relative_sat1 = range(FT(0), stop = FT(1), length = n_RS1)
    relative_sat2 = range(FT(1), stop = FT(1.02), length = n_RS2)
    relative_sat = [relative_sat1..., relative_sat2...]
    T_min = FT(150)
    T_surface = FT(350)

    args =
        fixed_lapse_rate_ref_state.(
            z_range,
            Ref(T_surface),
            Ref(T_min),
            Ref(param_set),
        )
    T, p, ρ = getindex.(args, 1), getindex.(args, 2), getindex.(args, 3)

    p = collect(Iterators.flatten([p for RS in 1:n_RS]))
    ρ = collect(Iterators.flatten([ρ for RS in 1:n_RS]))
    T = collect(Iterators.flatten([T for RS in 1:n_RS]))
    relative_sat = collect(Iterators.flatten([relative_sat for RS in 1:n]))

    # Additional variables
    q_sat = q_vap_saturation.(T, ρ, Ref(param_set))
    q_tot = min.(relative_sat .* q_sat, FT(1))
    q_pt = PhasePartition_equil.(T, ρ, q_tot, Ref(param_set))
    e_int = internal_energy.(T, q_pt, Ref(param_set))
    θ_liq_ice = liquid_ice_pottemp.(T, ρ, q_pt, Ref(param_set))
    return e_int, ρ, q_tot, q_pt, T, p, θ_liq_ice
end

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
# ---
# title: 1289. Minimum Falling Path Sum II
# id: problem1289
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/minimum-falling-path-sum-ii/description/>
# hidden: true
# ---
#
# Given a square grid of integers `arr`, a _falling path with non-zero shifts_
# is a choice of exactly one element from each row of `arr`, such that no two
# elements chosen in adjacent rows are in the same column.
#
# Return the minimum sum of a falling path with non-zero shifts.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [[1,2,3],[4,5,6],[7,8,9]]
#     Output: 13
#     Explanation:
#     The possible falling paths are:
#     [1,5,9], [1,5,7], [1,6,7], [1,6,8],
#     [2,4,8], [2,4,9], [2,6,7], [2,6,8],
#     [3,4,8], [3,4,9], [3,5,7], [3,5,9]
#     The falling path with the smallest sum is [1,5,7], so the answer is 13.
#
#
#
#
# **Constraints:**
#
#   * `1 <= arr.length == arr[i].length <= 200`
#   * `-99 <= arr[i][j] <= 99`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
import FinanceLib as Fl

import FinanceLib.FixedIncomes.MoneyMarkets as MM
import FinanceLib.FixedIncomes as FI

@testset "FinanceLib.FixedIncomes                                     " begin

  @testset "MoneyMarkets" begin
    @test MM.tBillR(150,98_000,100_000) ≈ 0.048
    @test MM.tBillD(0.048, 150, 100_000) == 2_000

    @test MM.holdingPerYield(98, 95, 5) == 0.020408163265306145
    @test MM.effAnnYield(150, 98, 95, 5) == 0.05038831660532006
    @test MM.moneyMktYield(150, 98, 95, 5) == 0.04897959183673475

    @test MM.twrr([4,6,5.775,6.72,5.508],[1,-0.5,0.225,-0.6]) ==
    0.06159232319186159
    @test MM.twrr(1.0, [100, 112, 142.64], [0, 20.])==0.21027878787878795

  end

end

include("Bonds/runtests.jl")


text:
export VariableBoundEntropy

struct VariableBoundEntropy end

function score(::Type{T}, objective::BetheFreeEnergy, ::VariableBoundEntropy, variable::RandomVariable, scheduler) where { T <: InfCountingReal }
    mapping = let d = degree(variable)
        (m) -> convert(T, (d - 1) * score(DifferentialEntropy(), m))
    end
    return getmarginal(variable, marginal_skip_strategy(objective)) |> schedule_on(scheduler) |> map(T, mapping)
end