text:
#=
--------------------------------------------------------------------------------------------------------------------------------
UC models: general interface
--------------------------------------------------------------------------------------------------------------------------------
=#

"""
    penalty_eigen(λ::Float64)

Return penalty value for a single eigenvalue λ.
"""
penalty_eigen(λ::Float64) = abs(λ) < 1 ? abs(λ)/(1-abs(λ)) : Inf;
penalty_eigen(λ::Complex{Float64}) = abs(λ) < 1 ? abs(λ)/(1-abs(λ)) : Inf;

"""
    fmin_uc_models(θ_unbound::FloatVector, lb::FloatVector, ub::FloatVector, transform_id::Array{Int64,1}, model_structure::Function, settings::UCSettings)

Return fmin for the UC model specified by model_structure(settings).

# Arguments
- `θ_unbound`: Model parameters (with unbounded support)
- `lb`: Lower bound for the parameters
- `ub`: Upper bound for the parameters
- `transform_id`: Type of transformation required for the parameters (0 = none, 1 = generalised log, 2 = generalised logit)
- `model_structure`: Function to setup the state-space structure
- `settings`: Settings for model_structure
- `tightness`: Controls the strength of the penalty (if any)
"""
function fmin_uc_models(θ_unbound::FloatVector, lb::FloatVector, ub::FloatVector, transform_id::Array{Int64,1}, model_structure::Function, uc_settings::UCSettings, tightness::Float64)

    # Compute parameters with bounded support
    θ = copy(θ_unbound);
    for i=1:length(θ)
        if transform_id[i] == 1
            θ[i] = get_bounded_log(θ_unbound[i], lb[i]);
        elseif transform_id[i] == 2
            θ[i] = get_bounded_logit(θ_unbound[i], lb[i], ub[i]);
        end
    end

    # Kalman status and settings
    status = KalmanStatus();
    model_instance, model_penalty = model_structure(θ, uc_settings);

    if ~isinf(model_penalty)
        settings = ImmutableKalmanSettings(model_instance...);

        # Compute loglikelihood for t = 1, ..., T
        for t=1:size(settings.Y,2)
            kfilter!(settings, status);
        end

        # Return fmin
        return -status.loglik + tightness*model_penalty;
    else
        return 1/eps();
    end
end

"""
    forecast(settings::KalmanSettings, h::Int64)

Compute the h-step ahead forecast for the data included in settings.

# Arguments
- `settings`: KalmanSettings struct
- `h`: Forecast horizon

    forecast(settings::KalmanSettings, X::FloatVector, h::Int64)

Compute the h-step ahead forecast for the data included in settings, starting from X.

# Arguments
- `settings`: KalmanSettings struct
- `X`: Last known value of the latent states
- `h`: Forecast horizon
"""
function forecast(settings::KalmanSettings, h::Int64)

    # Initialise Kalman status
    status = KalmanStatus();

    # Compute the period referring to the last observation of each series
    last_observations = zeros(settings.n) |> Array{Int64,1};
    for i=1:settings.n
        last_observations[i] = findall(.~ismissing.(settings.Y[i,:]))[end];
    end

    # Starting point for the forecast
    starting_point = minimum(last_observations);

    # Filter for t=1,...,T
    for t=1:settings.T
        kfilter!(settings, status);
    end

    # Initial forecast: series with a shorter history are forecasted until they match the others
    fc = zeros(settings.m, settings.T-starting_point+h);
    if starting_point < settings.T
        fc[:,1:settings.T-starting_point] = hcat(status.history_X_post[starting_point+1:settings.T]...);
    end

    # h-steps ahead forecast of the states from last observed point
    fc[:,settings.T-starting_point+1:end] = hcat(kforecast(settings, status.X_post, h)...);

    # Compute forecast for Y
    Y_fc = settings.B*fc;
    if settings.T-starting_point > 0
        Y_fc[last_observations.==settings.T, settings.T-starting_point] .= NaN;
    end

    # Return forecast for Y
    return Y_fc;
end

forecast(settings::KalmanSettings, X::FloatVector, h::Int64) = settings.B*hcat(kforecast(settings, X, h)...);

#=
--------------------------------------------------------------------------------------------------------------------------------
VARIMA model
--------------------------------------------------------------------------------------------------------------------------------
=#

"""
    varma_structure(θ::FloatVector, settings::VARIMASettings)

VARMA(p,q) representation similar to the form reported in Hamilton (1994) for ARIMA(p,q) models.

# Arguments
- `θ`: Model parameters
- `settings`: VARIMASettings struct
"""
function varma_structure(θ::FloatVector, settings::VARIMASettings)

    # Initialise
    ϑ = copy(θ);
    I_n = Matrix(I, settings.n, settings.n) |> FloatMatrix;
    UT_n = UpperTriangular(ones(settings.n, settings.n)) |> FloatMatrix;
    UT_n[I_n.==1] .= 0;

    # Observation equation
    B = [I_n reshape(ϑ[1:settings.nnq], settings.n, settings.nq) zeros(settings.n, settings.nr-settings.nq-settings.n)];
    R = Symmetric(I_n*1e-8);

    # Transition equation: coefficients
    C = [reshape(ϑ[settings.nnq+1:settings.nnq+settings.nnp], settings.n, settings.np) zeros(settings.n, settings.nr-settings.np);
         Matrix(I, settings.nr-settings.n, settings.nr-settings.n) zeros(settings.nr-settings.n, settings.n)];

    # VARMA(p,q) var-cov matrix
    V1 = Diagonal(ϑ[settings.nnq+settings.nnp+1:settings.nnq+settings.nnp+settings.n]) |> FloatMatrix;

    # Transition equation: variance
    V = Symmetric(cat(dims=[1,2], V1, zeros(settings.nr-settings.n, settings.nr-settings.n)));

    # Companion form for the moving average part
    companion_vma = [B[:,settings.n+1:settings.n+settings.nq];
                     Matrix(I, settings.nq-settings.n, settings.nq-settings.n) zeros(settings.nq-settings.n, settings.n)];

    # Compute penalty
    varma_penalty = 0.0;
    for λ = [eigvals(C); eigvals(companion_vma)]
        varma_penalty += penalty_eigen(λ);
    end

    # Return state-space structure
    return (settings.Y, B, R, C, V), varma_penalty;
end

"""
    varima(θ::FloatVector, settings::VARIMASettings)

Return KalmanSettings for a varima(d,p,q) model with parameters θ.

# Arguments
- `θ`: Model parameters
- `settings`: VARIMASettings struct

    varima(settings::VARIMASettings, args...; tightness::Float64=1.0)

Estimate varima(d,p,q) model.

# Arguments
- `settings`: VARIMASettings struct
- `args`: Arguments for Optim.optimize
- `tightness`: Controls the strength of the penalty for the non-causal / non-invertible case (default = 1)
"""
function varima(θ::FloatVector, settings::VARIMASettings)

    # Compute state-space parameters
    model_instance, _ = varma_structure(θ, settings);
    output = ImmutableKalmanSettings(model_instance...);

# TBD: update the warnings
#=
    # Warning 1: invertibility (in the past)
    eigval_ma = eigvals(companion_form(output.B[2:end]));
    if maximum(abs.(eigval_ma)) >= 1
        @warn("Invertibility (in the past) is not properly enforced! \n Re-estimate the model increasing the degree of differencing.");
    end

    # Warning 2: causality
    eigval_ar = eigvals(output.C);
    if maximum(abs.(eigval_ar)) >= 1
        @warn("Causality is not properly enforced! \n Re-estimate the model increasing the degree of differencing.");
    end

    # Warning 3: parameter redundancy
    intersection_ar_ma = intersect(eigval_ar, eigval_ma);
    if length(intersection_ar_ma) > 0
        @warn("Parameter redundancy! \n Check the AR and MA polynomials.");
    end
=#

    # Return output
    return output
end

function varima(settings::VARIMASettings, args...; tightness::Float64=1.0)

    # Starting point
    θ_starting = 1e-4*ones(settings.nnp+settings.nnq+settings.n);

    # Bounds
    lb = [-0.99*ones(settings.nnp+settings.nnq); 1e-6*ones(settings.n)];
    ub = [0.99*ones(settings.nnp+settings.nnq);  Inf*ones(settings.n)];
    transform_id = [2*ones(settings.nnp+settings.nnq); ones(settings.n)] |> Array{Int64,1};

    # Estimate the model
    res = Optim.optimize(θ_unbound->fmin_uc_models(θ_unbound, lb, ub, transform_id, varma_structure, settings, tightness), θ_starting, args...);

    # Apply bounds
    θ_minimizer = copy(res.minimizer);
    for i=1:length(θ_minimizer)
        if transform_id[i] == 1
            θ_minimizer[i] = get_bounded_log(θ_minimizer[i], lb[i]);
        elseif transform_id[i] == 2
            θ_minimizer[i] = get_bounded_logit(θ_minimizer[i], lb[i], ub[i]);
        end
    end

    # Return output
    return varima(θ_minimizer, settings);
end

"""
    arima(settings::VARIMASettings, args...; tightness::Float64=1.0)

Define an alias of the varima function for arima models.
"""
arima(settings::VARIMASettings, args...; tightness::Float64=1.0) = varima(settings, args..., tightness=tightness);

"""
    forecast(settings::KalmanSettings, h::Int64, varima_settings::VARIMASettings)

Compute the h-step ahead forecast for the data included in settings (in the varima_settings.Y_levels scale).

# Arguments
- `settings`: KalmanSettings struct
- `h`: Forecast horizon
- `varima_settings`: VARIMASettings struct
"""
function forecast(settings::KalmanSettings, h::Int64, varima_settings::VARIMASettings)

    # VARIMA
    if varima_settings.d > 0

        Y = zeros(varima_settings.d, varima_settings.n);

        # Loop over each series
        for i=1:varima_settings.n

            # Initialise Y_all
            Y_all = zeros(varima_settings.d, varima_settings.d);

            # Last observed point
            last_observation = findall(.~ismissing.(varima_settings.Y_levels[i,:]))[end];

            # The first row of Y_all is the data in levels
            Y_all[1,:] = varima_settings.Y_levels[i, last_observation-varima_settings.d+1:last_observation];

            # Differenced data, ex. (1-L)^d * Y_levels
            for j=1:varima_settings.d-1
                Y_all[1+j,:] = [NaN * ones(1,j) permutedims(diff(Y_all[j,:]))];
            end

            # Cut Y_all
            Y[:,i] = permutedims(Y_all[:,end]);
        end

        # Initial cumulated forecast
        fc_differenced = forecast(settings, h) .+ varima_settings.μ;
        fc = zeros(size(fc_differenced));

        # Loop over d to compute a prediction for the levels
        for i=1:varima_settings.n
            starting_point = findall(.~isnan.(fc_differenced[i,:]))[1];
            fc[i,starting_point:end] .= cumsum(fc_differenced[i,starting_point:end], dims=2);

            for j=varima_settings.d:-1:1
                fc[i,starting_point:end] .+= Y[j,i];
                if j != 1
                    fc[i,starting_point:end] = cumsum(fc[i,starting_point:end]);
                end
            end
        end

        # Insert NaNs
        fc[isnan.(fc_differenced)] .= NaN;

    # VARMA
    else
        # Compute forecast for varima_settings.Y (adjusted by its mean)
        fc = forecast(settings, h) .+ varima_settings.μ;
    end

    return fc;
end


text:
mutable struct Student
    id::Int                     # Unique ID for student
    total_credits::Int          # The total number of credit hours the student has earned
    gpa::Float64                # The student's GPA
    total_points::Float64       # The total number of points the student has earned

    attributes::Dict            # A dictionary that can store any kind of student attribute
    stopout::Bool               # Indicates whether the student has stopped out. (False if the student has, True if still enrolled)
    stopsem::Bool               # The term the student stopped out.
    termcredits::Int            # The number of credits the student has enrolled in for a given term.
    performance::Dict           # Stores the grades the student has made in each course.
    graduated::Bool             # Indicates wheter the student has graduated.
    gradsem::Int                # The term the student has graduated.
    termpassed::Array{Int}      # An array that represents the term in which the student passed each course.

    # Constructor
    function Student(id::Int, attributes::Dict)
        this = new()
        this.id = id
        this.termcredits = 0
        this.performance = Dict()
        this.gpa = 0.0
        this.total_credits = 0
        this.total_points = 0
        this.attributes = attributes

        return this
    end
end

# Helper Functions

# Returns an array of students
function simpleStudents(number)
    students = Student[]
    for i=1:number
        student = Student(i, Dict())
        push!(students, student)
    end
    return students
end

ext: jl

max_stars_repo_name: M-S-H/CASL

max_issues_repo_name: M-S-H/CASL

max_forks_repo_name: M-S-H/CASL


text:
using ElectronicStructure
using Artifacts
import JLD2

"""
    get_molecular_data()

Return a `Dict` containing a few examples of PySCF electronic structure
calculations. They are retrieved in `JLD2` format using the `Artifacts`
framework. The `JLD2` file is read and the result is returned.

The data returned was generated using the function `generate_pyscf_molecules`
in *compare_formats.jl*.
"""
function get_molecular_data()
    rootpath = artifact"molecular_data"
    filename = joinpath(rootpath, "molecular_data_examples.jld2")
    data = JLD2.load(filename)
    return data
end


text:
module BumpOnTailSimulation

using LaTeXStrings
using LinearAlgebra
using Particles
using Particles.BumpOnTail
using Plots

# simulation parameters

const dt = 1e-1             # timestep
const T = 50                # final time
const nₜ = Int(div(T, dt))  # nb. of timesteps
const nₕ = 16               # nb. of elements
const p  = 3                # spline degree

const nₚ = Int(5e4)         # nb. of particles

const vmin = -10
const vmax = +10

# bump-on tail instability parameters: κ, ε, a, v₀, σ
# in order: spartial perturbation wave number and amplitude of s.p., fast particle fraction, speed, and temperature
const params = (κ = 0.3, ε = 0.03, a = 0.1, v₀= 4.5, σ = 0.5, χ = 1.0)

const L = 2π / params.κ     # domain length
const h = L/nₕ              # element width

const IP = VPIntegratorParameters(dt, nₜ, nₜ+1, nₕ, nₚ)

function run()

    # B-spline Poisson solver
    poisson = PoissonSolverPBSplines(p, IP.nₕ, L)

    # initial data
    P = BumpOnTail.draw_accept_reject(nₚ, params)
    # P = BumpOnTaildraw_importance_sampling(nₚ, params)

    # initial potential
    solve!(poisson, P.x, P.w)

    # plot initial particle distribution
    F = plot_particles(P.x, P.v, P.w, 0, L, vmin, vmax)
    savefig("bump-on-tail-init.png")

    F = plot_distribution(P.x, P.v, P.w, 0, L, vmin, vmax; nx=48, nv=32)
    savefig("bump-on-tail-init-f.png")

    # create integrate cache
    IC = VPIntegratorCache(IP)

    # integrate all time steps
    # for t = 1:nₜ
        integrate_vp!(P, poisson, params, IP, IC)
    # end

    # plot diagnostics
    plot(0:dt:dt*nₜ, IC.W, yaxis = :log, legend = :none)
    savefig("bump-on-tail-W.png")

    plot(0:dt:dt*nₜ, IC.K, yaxis = :log, legend = :none)
    savefig("bump-on-tail-K.png")

    plot(0:dt:dt*nₜ, IC.W .+ IC.K, yaxis = :log, legend = :none)
    savefig("bump-on-tail-E.png")

    # create animation and save to file
    anim = @animate for t in 0:IP.nₜ
        plot_particles(IC.X[:,t+1], IC.V[:,t+1], P.w, 0, L, vmin, vmax)
    end
    gif(anim, "bump-on-tail.gif", fps=10)
end

end # module

using .BumpOnTailSimulation

BumpOnTailSimulation.run()


text:
__precompile__(false)
module FluidFlowSimulation
export run_simulation, advance_in_time!, parameters, readglobal, initialize!

using FFTW, InplaceRealFFT, FluidTensors, FluidFields, LinearAlgebra, GlobalFileHelper, Interpolations, DelimitedFiles

include("Globals.jl")

using .Globals

include("macros.jl")
include("simulationtypes.jl")
include("Forcing_methods/spectrum.jl")
include("utils.jl")
include("vectorfunctions.jl")
include("time_step_functions/time_step_functions.jl")
include("stats.jl")
include("spec.jl")
include("output.jl")
include("models.jl")

function run_simulation()
    par = readglobal()
    s = parameters(par)
    lengthtime = haskey(par,:timeDuration) ? parse(Float64,par[:timeDuration]) : Inf
    nt = haskey(par,:nt) ? parse(Int,par[:nt]) : (typemax(Int) - s.iteration[])
    run_simulation(s,nt,lengthtime)
end

@par function run_simulation(s::@par(AbstractSimulation),totalnsteps::Integer,lenghtime::Real)
    @info(s)
    @info("Simulation started.")

    initialize!(s)

    @assert totalnsteps + s.iteration[] >= s.dtoutput
    @assert totalnsteps + s.iteration[] >= s.dtstats
    @show finalstep = s.iteration[] + totalnsteps
    @show finaltime = s.time[]+lenghtime

    while (s.iteration[]<finalstep && s.time[]<finaltime)
        advance_in_time!(s)
    end

    calculate_rhs!(s)
    hasforcing(s) && s.forcing(s)

    writestats(s)

    if !(mod(s.iteration[],s.dtoutput) == 0)
        real!(s.u)
        hasdensity(s) && real!(s.densitystratification.ρ)
        haspassivescalar(s) && real!(s.passivescalar.φ)
        hasles(s) && real!(s.lesmodel.tau)

        writeoutput(s)
    end

    if typeof(s.forcing) <: RfForcing
        writedlm("R.$(s.iteration[])",s.forcing.R)
    end

    return s
end

function initialize!(s::AbstractSimulation)
    s.iteration[] == 0 && writeheader(s)
    myfourier!(s.u)
    #dealias!(s.u)

    if haspassivescalar(s)
        myfourier!(s.passivescalar.φ)
        dealias!(s.passivescalar.φ)
    end

    if hasdensity(s)
        myfourier!(s.densitystratification.ρ)
        dealias!(s.densitystratification.ρ)
    end

    calculate_rhs!(s)

    initialize!(s.forcing,s)
    initialize!(s.timestep,s.rhs,s)
    initialize!(s.passivescalar,s)
    initialize!(s.densitystratification,s)

    return nothing
end

end # module


text:
using Plots

function square(c,side)
  x = [ c[1]-side/2, c[1]+side/2, c[1]+side/2, c[1]-side/2, c[1]-side/2]
  y = [ c[2]-side/2, c[2]-side/2, c[2]+side/2, c[2]+side/2, c[2]-side/2]
  return x, y
end

# Generate points
p = [ 10*rand(2) for i in 1:100 ]

# Plot a square around point 1
plot(square(p[1],10),seriestype=[:shape],
            linewdith=0.5,fillalpha=0.2,
            color="green",label="")

# Images
images = [ -2, -1, 0, 1, 2 ]
for j in images, k in images
  if j != 0 || k != 0 # Do not plot (0,0) again
    xim = [ p[i][1] + j*10 for i in 1:100 ]
    yim = [ p[i][2] + k*10 for i in 1:100 ]
    scatter!(xim,yim,color="orange",label="")
  end
end

# Central plot
x = [ p[i][1] for i in 1:100 ]
y = [ p[i][2] for i in 1:100 ]
scatter!(x,y,xlabel="x",ylabel="y",label="",color="lightblue")

# Mark point one as red
scatter!([p[1][1]],[p[1][2]],color="red",label="",size=4)

plot!(xlim=[-15,25],ylim=[-15,25],
      size=(400,400),framestyle=:box)

hline!(10*images,color="black",label="") # horizontal lines
vline!(10*images,color="black",label="") # vertical lines

savefig("pbc2.pdf")


text:
context("SimpleRectangle") do
context("constructors") do
    s = SimpleRectangle(1,2,3,4)
    @fact typeof(s) --> SimpleRectangle{Int}
    @fact width(s) --> 3
    @fact height(s) --> 4
    @fact area(s) --> 12
    @fact xwidth(s) --> 4
    @fact yheight(s) --> 6
    @fact minimum(s) --> Point(1,2)
    @fact maximum(s) --> Point(4,6)
    @fact origin(s) --> Point(1,2)
    h = HyperRectangle(s)
    @fact typeof(h) --> HyperRectangle{2,Int}
    @fact h --> HyperRectangle(1,2,3,4)
    h = HyperRectangle{3,Int}(s)
    @fact h --> HyperRectangle(1,2,0,3,4,0)

    @fact SimpleRectangle(Vec(2,3)) --> SimpleRectangle(0,0,2,3)
end

context("setops") do
    r1 = SimpleRectangle(0,0,2,2)
    r2 = SimpleRectangle(2,2,2,2)
    r3 = SimpleRectangle(1,1,2,2)
    r4 = SimpleRectangle(0.5,0.5,3.,3.)
    @fact intersect(r1,r2) --> SimpleRectangle(2,2,0,0)
    @fact intersect(r1,r3) --> SimpleRectangle(1,1,1,1)
    @fact intersect(r1,r4) --> SimpleRectangle(0.5,0.5,1.5,1.5)
end

context("indexing") do
    r = SimpleRectangle(0,0,2,2)
    a = eye(4)
    @fact a[r] --> eye(2)
    a[r] = eye(2)*2
    @fact a --> [2 0 0 0; 0 2 0 0; 0 0 1 0; 0 0 0 1]
end

end


text:
@testset "RMSprop" begin
    run_optim_tests(RMSprop(),
                    skip = ("Trigonometric", "Powell", "Extended Powell"),
                    f_increase_exceptions = ("Hosaki",),
                    convergence_exceptions = (("Polynomial", 1), ("Polynomial", 2), ("Rosenbrock", 1), ("Rosenbrock", 2),
                                              ("Extended Rosenbrock", 1), ("Extended Rosenbrock", 2),
                                              ("Penalty Function I", 1),
                                              ("Penalty Function I", 2)),
                    iteration_exceptions = (("Rosenbrock", 10000),
                                            ("Extended Rosenbrock", 12000),
                                            ("Fletcher-Powell", 10000),
                                            ("Paraboloid Diagonal", 10000),
                                            ("Paraboloid Random Matrix", 10000),
                                            ("Penalty Function I", 10000),),
                    show_name = debug_printing)

    function f_gd_1(x)
      (x[1] - 5.0)^2
    end

    function g_gd_1(storage, x)
      storage[1] = 2.0 * (x[1] - 5.0)
    end

    initial_x = [0.0]

    d = OnceDifferentiable(f_gd_1, g_gd_1, initial_x)

    results = Optim.optimize(d, initial_x, RMSprop())
    @test_throws ErrorException Optim.x_trace(results)
    @test Optim.g_converged(results)
    @test norm(Optim.minimizer(results) - [5.0]) < 0.01
    @test summary(results) == "Root Mean Square propagation"

    eta = 0.9

    function f_gd_2(x)
      (1.0 / 2.0) * (x[1]^2 + eta * x[2]^2)
    end

    function g_gd_2(storage, x)
      storage[1] = x[1]
      storage[2] = eta * x[2]
    end

    d = OnceDifferentiable(f_gd_2, g_gd_2, [1.0, 1.0])

    results = Optim.optimize(d, [1.0, 1.0], RMSprop())
    @test_throws ErrorException Optim.x_trace(results)
    @test Optim.g_converged(results)
    @test norm(Optim.minimizer(results) - [0.0, 0.0]) < 0.01
end


text:
include("LAS.jl")

using .ListenAttendSpell

ListenAttendSpell.main(n_epochs=300, saved_results=true, log_weights=false)


text:
using IntervalContractors
using IntervalArithmetic
using Test

include("abs_rev.jl")
include("libieeep1788_rev.jl")
include("Non1788tests.jl")

@testset "power_rev for odd power" begin
    x = -∞..∞
    a = -8..27
    power_rev(a, x, 3)[2] == Interval(-2.0000000000000004, 3.0000000000000004)
end


text: # Nothing, at the moment


text:
using Documenter
using SequentialRaytrace

push!(LOAD_PATH,"../src/")

makedocs(sitename="SequentialRaytrace")


text:
module Multisets

import Base: show, length, getindex, collect, union, intersect, isempty, first
import Base: push!, setindex!, delete!, hash
import Base: (*), (+), (-), (|), (&)
import Base: issubset, Set, (==), (<), (<=), (>), (>=)

export Multiset, set_short_show, set_julia_show, set_braces_show

"""
A `Multiset` is an unordered collection of things with repetition permitted.
A new `Multiset` container is created with `Multiset{T}()` where `T` is the
type of the objects held in the multiset. If `T` is omitted, it defaults
to `Any`.

A `Multiset` can be created from a collection `list` (such as a `Vector` or
`Set`) with `Multiset(list)`. If an element is repeated in `list` it has
the appropriate multiplicity.
"""
struct Multiset{T}
  data::Dict{T,Int}
  function Multiset{T}() where T
    d = Dict{T,Int}()
    new(d)
  end
end
Multiset() = Multiset{Any}()

function Multiset(list::AbstractArray{T,d}) where {T,d}
  M = Multiset{T}()
  for x in list
    push!(M,x)
  end
  return M
end

function Multiset(A::Base.AbstractSet{T}) where T
  M = Multiset{T}()
  for x in A
    push!(M,x)
  end
  return M
end

"""
`clean!(M)` removes elements of multiplicy 0 from the underlying data
structure supporting `M`.
"""
function clean!(M::Multiset)
  for x in keys(M.data)
    if M[x]==0
      delete!(M.data,x)
    end
  end
  nothing
end

"""
For a `M[t]` where `M` is a `Multiset` returns the
multiplicity of `t` in `M`. A value of `0` means that
`t` is not a member of `M`.
"""
function getindex(M::Multiset{T}, x)::Int where T
  if haskey(M.data,x)
    return M.data[x]
  end
  return 0
end

"""
`push!(M,x,incr)` increases the multiplicity of `x` in `M`
by `incr` (which defaults to 1). `incr` can be negative, but
it is not possible to decrease the multiplicty below 0.
"""
function push!(M::Multiset{T}, x, incr::Int=1) where T
  if haskey(M.data,x)
    M.data[x] += incr
  else
    M.data[x] = incr
  end
  if M.data[x] < 0
    M.data[x]=0
  end
  return M
end

function setindex!(M::Multiset{T}, m::Int, x) where T
  mm = max(m,0)
  M.data[x] = mm
  return mm
end

function delete!(M::Multiset, x)
  if haskey(M.data,x)
    delete!(M.data,x)
  end
  return M
end

function length(M::Multiset)
  total = 0
  for v in values(M.data)
    total += v
  end
  return total
end

isempty(M::Multiset) = length(M)==0

function collect(M::Multiset{T}) where T
  n = length(M)
  result = Array{T,1}(undef,n)   #  Vector{T}(n)
  i = 0
  for (k,v) in M.data
    for _=1:v
      i += 1
      result[i] = k
    end
  end
  try
    sort!(result)
  catch
  end
  return result
end

function braces_string(M::Multiset{T}) where T
  elts = collect(M)
  n = length(elts)
  str = "{"
  for k=1:n
    str *= string(elts[k])
    if k<n
      str *= ","
    end
  end
  str *= "}"
  return str
end

short_string(M::Multiset{T})  where T= "Multiset{$T} with $(length(M)) elements"

function julia_string(M::Multiset{T}) where T
  elts = collect(M)
  n = length(elts)
  str = "Multiset($T["
  q = ""
  if T <: AbstractString
    q = "\""
  end
  for k=1:n
    str *= q*string(elts[k])*q
    if k<n
      str *= ","
    end
  end
  str *= "])"
  return str
end

# This variable controls printing:
# 0 -- {x,y,z}
# 1 -- Multiset{T} with n elements
# 2--  Multiset{T}(x,y,z)

const multi_show_braces = 0
const multi_show_short  = 1
const multi_show_julia  = 2
multi_show_flag = multi_show_braces

"""
Set show display mode for multisets, like this:

`Multiset{Int64} with 7 elements`

See also `set_braces_show` and `set_julia_show`.
"""
set_short_show() = (global multi_show_flag = multi_show_short; nothing)

"""
Set braces display mode for multisets, like this:

`{1,2,2,3,3,3,3}`

See also `set_short_show` and `set_julia_show`.
"""
set_braces_show() = (global multi_show_flag = multi_show_braces; nothing)

"""
Set Julia style display mode for multisets, like this:

`Multiset(Int64[1,2,2,3,3,3,3])`

See also `set_short_show` and `set_braces_show`.
"""
set_julia_show()  = (global multi_show_flag = multi_show_julia; nothing)

function show(io::IO, M::Multiset)
  if multi_show_flag == multi_show_short
    print(io, short_string(M))
  end
  if multi_show_flag == multi_show_braces
    print(io, braces_string(M))
  end
  if multi_show_flag == multi_show_julia
    print(io, julia_string(M))
  end
end

# private helper for union
function _union(A::Multiset{S}, B::Multiset{S}) where S
  M = Multiset{S}()
  for (x,v) in A.data
    M[x] = max(v,B[x])
  end
  for (y,w) in B.data
    M[y] = max(A[y],w)
  end
  return M
end

function type_convert(T, A::Multiset)
  M = Multiset{T}()
  for (x,v) in A.data
    xx = convert(T,x)
    M[xx] = v
  end
  return M
end

"""
`union(A,B)` for multisets creates a new multiset in which the
multiplicity of `x` is `max(A[x],B[x])`.
"""
function union(A::Multiset{S}, B::Multiset{T}) where {S,T}
  if S==T
    return _union(A,B)
  end
  ST = typejoin(S,T)
  AA = type_convert(ST,A)
  BB = type_convert(ST,B)
  return _union(AA,BB)
end

"""
For multisets, `A|B` is `union(A,B)`. See also `+` which behaves differently.
"""
(|)(A::Multiset,B::Multiset) = union(A,B)

"""
`A+B` for multisets is the disjoint union, i.e., a new multiset in which the
multiplicity of `x` is `A[x]+B[x]`. This can be abbreviated `A|B`.
"""
function (+)(A::Multiset{S}, B::Multiset{T}) where {S,T}
  ST = typejoin(S,T)
  M = Multiset{ST}()
  for (x,v) in A.data
    push!(M,x,v)
  end
  for (x,v) in B.data
    push!(M,x,v)
  end
  return M
end

# private helper for A-B
function _minus(A::Multiset{S},B::Multiset{S}) where S
  M = Multiset{S}()
  for (x,v) in A.data
    M.data[x] = max(v-B[x],0)
  end
  return M
end

"""
`A-B` for multisets is the multiset difference, i.e., a new multiset
in which the multiplicity of `x` is `A[x]-B[x]` unless this goes
below `0`, in which case the multiplicity is 0.
"""
function (-)(A::Multiset{S}, B::Multiset{T}) where {S,T}
  if S==T
    return _minus(A,B)
  end
  ST = typejoin(S,T)
  AA = type_convert(ST,A)
  BB = type_convert(ST,B)
  return _minus(AA,BB)
end

# private helper for intersect
function _inter(A::Multiset{S}, B::Multiset{S}) where S
  M = Multiset{S}()
  for (x,v) in A.data
    push!(M,x,min(v,B[x]))
  end
  return M
end

"""
`intersect(A,B)` for multisets creates a new multiset in which the
multiplicity of `x` is `min(A[x],B[x])`.
This may be abbreviated `A&B`.
"""
function intersect(A::Multiset{S}, B::Multiset{T}) where {S,T}
  if S==T
    return _inter(A,B)
  end
  ST = typejoin(S,T)
  AA = type_convert(ST,A)
  BB = type_convert(ST,B)
  return _inter(AA,BB)
end

"""
`A&B` for multisets is `intersect(A,B)`.
"""
(&)(A::Multiset,B::Multiset) = intersect(A,B)

"""
`A*B` for the Cartesian product of multisets `A` and `B`.
"""
function (*)(A::Multiset{S}, B::Multiset{T}) where {S,T}
  ST = Tuple{S,T}
  M = Multiset{ST}()
  for (a,v) in A.data
    for (b,w) in B.data
      M.data[(a,b)] = v*w
    end
  end
  return M
end

"""
`n*A` is the scalar multiple of a multiset in which the multiplicity of
`x` is `n*A[x]`. Of course, we require `n >= 0`.
"""
function (*)(n::Int, A::Multiset{T}) where T
  @assert n>=0 "Scalar multiplication of a multiset must be by a nonnegative integer"
  M = Multiset{T}()
  for (x,v) in A.data
    M.data[x] = n*v
  end
  return M
end

# private helper for issubset
function _sub(A::Multiset{S}, B::Multiset{S}) where S
  for (x,v) in A.data
    if v > B[x]
      return false
    end
  end
  return true
end

function issubset(A::Multiset{S}, B::Multiset{T}) where {S,T}
  if S==T
    return _sub(A,B)
  end
  ST = typejoin(S,T)
  AA = type_convert(ST,A)
  BB = type_convert(ST,B)
  return _sub(AA,BB)
end

(<=)(A::Multiset, B::Multiset) = issubset(A,B)
(<)(A::Multiset, B::Multiset) = (length(A)<length(B)) && (A<=B)
(>)(A::Multiset, B::Multiset) = B<A
(>=)(A::Multiset, B::Multiset) = B<=A

(==)(A::Multiset, B::Multiset) = (length(A)==length(B)) && issubset(A,B)

function hash(M::Multiset, h::UInt = UInt(0))
  clean!(M)
  return hash(M.data,h)
end

function Set(M::Multiset{T}) where T
  iter = (x for x in keys(M.data) if M.data[x]>0)
  return Set{T}(iter)
end

function first(M::Multiset{T})::T where T
    if length(M)==0
        error("Multiset must be nonempty")
    end
    clean!(M)
    return first(first(M.data))
end

end #end of Module


text:
addprocs(1)

using BaseBenchmarks
using BenchmarkTools
using Base.Test

BaseBenchmarks.loadall!()

@test begin
    run(BaseBenchmarks.SUITE, verbose = true, samples = 1,
        evals = 2, gctrial = false, gcsample = false);
    true
end


text:
# Solvers
abstract type Solver1D end

struct SolverInfo
  maxEvals::Int
  lowerBoundEnforced::Bool
  upperBoundEnforced::Bool
  lowerBound::Float64
  upperBound::Float64
end

# solver functions
function solve(solver::Solver1D, f::Function, accuracy::Float64, guess::Float64, step::Float64)
  ## This method returns the 0 of a function determined by a given accuracy.  This method using bracketing
  ## routine to which an intial guess must be supplied as well as a step
  growth_factor = 1.6
  flipflop = -1

  root = guess
  fxMax = f(root)

  # monotonically crescent bias, as in optionValue(volatility)
  if is_close(fxMax, 0.0)
    return root
  elseif fxMax > 0.0
    xMin = enforced_bounds(solver, root - step)
    fxMin = f(xMin)
    xMax = root
  else
    xMin = root
    fxMin = fxMax
    xMax = enforced_bounds(solver, root + step)
    fxMax = f(xMax)
  end

  eval_num = 2

  while eval_num < solver.solverInfo.maxEvals
    if fxMin * fxMax <= 0.0
      if is_close(fxMin, 0.0)
        return xMin
      end
      if is_close(fxMax, 0.0)
        return xMax
      end
      root = (xMax + xMin) / 2.0
      return _solve(solver, f, accuracy, xMin, xMax, fxMin, fxMax, root, eval_num)
    end

    if abs(fxMin) < abs(fxMax)
      xMin = enforced_bounds(solver, xMin + growth_factor * (xMin - xMax))
      fxMin = f(xMin)
    elseif abs(fxMin) > abs(fxMax)
      xMax = enforced_bounds(solver, xMax + growth_factor * (xMax - xMin))
      fxMax = f(xMax)
    elseif flipflop == -1
      xMin = enforced_bounds(solver, xMin + growth_factor * (xMin - xMax))
      fxMin = f(xMin)
      eval_num += 1
      flipflop = 1
    elseif flipflop == 1
      xMax = enforced_bounds(solver, xMax + growth_factor * (xMax - xMin))
      fxMax = f(xMax)
      flipflop = -1
    end

    eval_num += 1
  end

  return root
end

function solve(solver::Solver1D, f::Function, accuracy::Float64, guess::Float64, xMin::Float64, xMax::Float64)
  fxMin = f(xMin)
  if is_close(fxMin, 0.0)
    return xMin
  end

  fxMax = f(xMax)
  if is_close(fxMax, 0.0)
    return xMax
  end

  eval_num = 2

  return _solve(solver, f, accuracy, xMin, xMax, fxMin, fxMax, guess, eval_num)
end

function enforced_bounds(solver::Solver1D, x::Float64)
  if solver.solverInfo.lowerBoundEnforced && x < solver.solverInfo.lowerBound
    return solver.solverInfo.lowerBound
  end

  if solver.solverInfo.upperBoundEnforced && x > solver.solverInfo.upperBound
    return solver.solverInfo.upperBound
  end

  return x
end


text:
using JLD, HDF5, GZip, StatsBase, LinearAlgebra

function mylogsumexp(x)
	X = maximum(x, dims = 2)
	return log.(sum(exp.(broadcast(-, x, X)), dims = 2)) + X
end

function convert_protein(prot_seq)
	prot_len = length(prot_seq)
	big_vector = zeros(Float32, prot_len * 21)
	#and fill it up like a very sparse row vector. One-hot.
	#Order is ARNDCQEGHILKMFPSTWYV-
	prot_fill = Dict{Char, Int64}('A'=>1, 'R'=>2, 'N'=>3, 'D'=>4, 'C'=>5, 'Q'=>6, 'E'=>7, 'G'=>8, 'H'=>9, 'I'=>10,
										'L'=>11, 'K'=>12, 'M'=>13, 'F'=>14, 'P'=>15, 'S'=>16, 'T'=>17, 'W'=>18, 'Y'=>19, 'V'=>20, '-'=>21, 'X'=>21, '*'=>21)
	chr_count = 0
	for chr in prot_seq
		big_vector[(21 * chr_count) + prot_fill[chr]] = Float32(1.0)
		chr_count += 1
	end
	return big_vector' #more convenient for it to be a row vector later.
end

function compute_energies(prot_mat, w1, w2)
	println("Computing energies")
	energies = prot_mat * w1 + diag(prot_mat * w2 * prot_mat')
	return energies #maybe this works?
end

function compute_psls(prot_mat, w1, w2, nNodes, nProt)
	println("Computing pseudolikelihoods")
	my_pv_w2 = prot_mat * w2
	my_ull = my_pv_w2 * prot_mat'
	my_pv_w1 = prot_mat * w1
	numerator = my_pv_w1 + diag(my_ull)
	partitions = reshape((w1' .+ my_pv_w2)', 21, nNodes * nProt)
	logZs = sum(reshape(mylogsumexp(partitions'), nNodes, nProt), dims = 1)'
	all_scores = numerator - logZs
	return -all_scores
end

function save_energies(computed_energies, gene_name)
    local out_file
    try
	    out_file = open("../main/energies/energy_$(gene_name).txt", "w")
    catch err
        if isa(err, SystemError)
            out_file = open("energy_$(gene_name).txt", "w")
            println("Warning: writting energy output to working directory")
        else
            rethrow()
        end
    finally
	    for energy in computed_energies
		    write(out_file, "$energy\n")
	    end
	    close(out_file)
    end
end

function save_psls(computed_psls, gene_name)
    local out_file
    try
	    out_file = open("../main/psls/psls_$(gene_name).txt", "w")
    catch err
        if isa(err, SystemError)
            out_file = open("psls_$(gene_name).txt", "w")
            println("Warning: writting pseudolhood output to working directory")
        else
            rethrow()
        end
    finally
	    for psl in computed_psls
		    write(out_file, "$psl\n")
	    end
	    close(out_file)
    end
end

function run(gene_name, jld_file, msa_file)
	println("Reading msa file")
	in_file = open(msa_file)
	in_read = readlines(in_file)
	close(in_file)

	seqs = [uppercase(strip(xx)) for xx in filter(x -> x[1] != '>', in_read)]
	seqs = filter(seq -> !('B' in seq) && !('J' in seq) && !('O' in seq) && !('U' in seq) && !('Z' in seq), seqs) #try to avoid annoying sequences in alignment that have non-aas. X's are sorta okay.

	num_prot = length(seqs)
	num_aa = length(seqs[1]) #assume MSA is well formed.

	prot_mat = zeros(Float32, num_prot, num_aa * 21)
	i_count = 0
	for seq in seqs
		seq = uppercase(seq)
		i_count += 1
		prot_mat[i_count, :] = convert_protein(seq)
	end

	println("Reading jld file")
	mrf = load(jld_file)
	w1 = mrf["w1"]
	w2 = mrf["w2"]

	@time computed_energies = compute_energies(prot_mat, w1, w2)
	@time computed_psls = compute_psls(prot_mat, w1, w2, num_aa, num_prot)

	#Write outputs to main directory.
	println("Done. Saving energies/psls.")

    save_energies(computed_energies, gene_name)
    save_psls(computed_psls, gene_name)
end

function main()
	prot_name = ARGS[end - 2]
	prot_jld = ARGS[end - 1]
	prot_msa = ARGS[end]

	if isfile(prot_msa)
		if isfile(prot_jld) && endswith(prot_jld, ".jld")
			run(prot_name, prot_jld, prot_msa)
		else
			println("Something wrong with $prot_jld ... not a file or doesn't end with .jld")
		end
	else
		println("$prot_msa is not a file.")
	end
end

main()

max_stars_repo_name: khyox/CAMEOS

max_issues_repo_name: BiosecSFA/CAMEOS

max_forks_repo_name: BiosecSFA/CAMEOS


text:
function variablecost(phy::JumpVariable, devices::Array{T}) where T <: PowerSystems.HydroGen

    cost = JuMP.AffExpr()

    for (ix, name) in enumerate(pre.axes[1])
        if name == devices[ix].name
                JuMP.add_to_expression!(cost,precost(phy[name,:], devices[ix]))
        else
            @error "Bus name in Array and variable do not match"
        end
    end

    return cost

end

function precost(X::JuMP.VariableRef, device::Union{PowerSystems.RenewableCurtailment,PowerSystems.RenewableFullDispatch})

    return cost = sum(device.econ.curtailcost*(-X))
end


text:
# Code for working with UncertainValues interpreted as measures of composition.

"""
    MaterialLabel

The abstract type associated with `Label`s with `material` and `element` members.
"""
abstract type MaterialLabel <: Label end

"""
    materiallabels(ty::Type{<:MaterialLabel}, material::String, uvs::Union{LabeledValues,UncertainValues})

Useful for extracting the MaterialLabel[s] associated with a single material.
"""
materiallabels(
    ty::Type{<:MaterialLabel},
    material::String,
    luvs::Union{LabeledValues,UncertainValues},
) = filter(l -> (l isa ty) && (l.material == material), labels(luvs))

struct MassFractionLabel <: MaterialLabel
    material::String
    element::Element
end
prefix(mfl::MassFractionLabel) = "C"

struct NormMassFractionLabel <: MaterialLabel
    material::String
    element::Element
end
prefix(mfl::NormMassFractionLabel) = "N"

struct AtomicFractionLabel <: MaterialLabel
    material::String
    element::Element
end
prefix(mfl::AtomicFractionLabel) = "A"

struct AtomicWeightLabel <: MaterialLabel
    material::String
    element::Element
end
prefix(awl::AtomicWeightLabel) = "Az"

Base.isequal(el1::MaterialLabel, el2::MaterialLabel) =
    isequal(typeof(el1), typeof(el2)) &&
    isequal(prefix(el1), prefix(el2)) && #
    isequal(el1.element, el2.element) &&
    isequal(el1.material, el2.material)

Base.isless(el1::MaterialLabel, el2::MaterialLabel) =
    isequal(prefix(el1), prefix(el2)) ?
    (
        isequal(el1.material, el2.material) ? isless(el1.element, el2.element) :
        isless(el1.material, el2.material)
    ) : isless(prefix(el1), prefix(el2))

Base.show(io::IO, el::MaterialLabel) =
    print(io, "$(prefix(el))[$(el.element.symbol),$(el.material)]")

"""
    MaterialFraction

This label represents the amount of a material `constituent` in `material` in a material
defined as the mixture of other materials.
"""
struct MaterialFractionLabel <: Label
    material::String
    constituent::String
end

Base.show(io::IO, mf::MaterialFractionLabel) =
    print(io, "f[", mf.material, ",", mf.constituent, "]")

struct MFtoAF <: MeasurementModel
    material::String
end

function NeXLUncertainties.compute(
    ma::MFtoAF,
    inputs::LabeledValues,
    withJac::Bool,
)::MMResult
    mfls = materiallabels(MassFractionLabel, ma.material, inputs)
    awls = map(mfl -> AtomicWeightLabel(ma.material, mfl.element), mfls)
    scoa = sum(i -> inputs[mfls[i]] / inputs[awls[i]], eachindex(mfls))
    results = map(i -> (inputs[mfls[i]] / inputs[awls[i]]) / scoa, eachindex(mfls))
    outputs = map(mfl -> AtomicFractionLabel(ma.material, mfl.element), mfls)
    vals = LabeledValues(outputs, results)
    jac = withJac ? zeros(Float64, length(outputs), length(inputs)) : missing
    if withJac
        for i in eachindex(outputs)
            for j in eachindex(mfls)
                if i == j
                    @assert outputs[i].element == mfls[j].element
                    jac[i, indexin(inputs, mfls[i])] =
                        results[i] * (1.0 - results[i]) / inputs[mfls[i]]
                    jac[i, indexin(inputs, awls[i])] =
                        results[i] * (results[i] - 1.0) / inputs[awls[i]]
                else
                    jac[i, indexin(inputs, mfls[j])] =
                        -results[i] * results[j] / inputs[mfls[j]]
                    jac[i, indexin(inputs, awls[j])] =
                        results[i] * results[j] / inputs[awls[j]]
                end
            end
        end
    end
    return (vals, jac)
end

"""
    MFtoNMF

Mass fraction to normalized mass fraction measurement model.
"""
struct MFtoNMF <: MeasurementModel
    material::String
end

function NeXLUncertainties.compute(
    ma::MFtoNMF,
    inputs::LabeledValues,
    withJac::Bool,
)::MMResult
    # Extract the elements for which there is mass-fraction data associated with ma.material
    elms = map(
        il -> il.element,
        filter(
            l -> (l isa MassFractionLabel) && isequal(ma.material, l.material),
            labels(inputs),
        ),
    )
    outputs = map(elm -> NormMassFractionLabel(ma.material, elm), elms)
    mfls = map(elm -> MassFractionLabel(ma.material, elm), elms)
    sc = sum(i -> inputs[mfls[i]], eachindex(elms))
    results = map(i -> inputs[mfls[i]] / sc, eachindex(elms))
    jac = withJac ? zeros(Float64, length(outputs), length(inputs)) : missing
    if withJac
        for i in eachindex(outputs)
            for j in eachindex(elms)
                if i == j
                    @assert outputs[i].element == elms[j]
                    jac[i, indexin(inputs, mfls[i])] =
                        results[i] * (1.0 - results[i]) / inputs[mfls[i]]
                else
                    jac[i, indexin(inputs, mfls[j])] =
                        -results[i] * results[j] / inputs[mfls[j]]
                end
            end
        end
    end
    return (LabeledValues(outputs, results), jac)
end

"""
    AFtoNMF

Converts atomic fraction into mass fraction - since the results in by necessity normalized, both the
MassFractionLabel and NormMassFractionLabel versions are populated with identical information.
"""
struct AFtoNMF <: MeasurementModel
    material::String
end

function NeXLUncertainties.compute(
    atm::AFtoNMF,
    inputs::LabeledValues,
    withJac::Bool,
)::MMResult
    # Extract the elements for which there is atomic-fraction data associated with atm.material
    elms = map(
        il -> il.element,
        filter(
            l -> (l isa AtomicFractionLabel) && isequal(atm.material, l.material),
            labels(inputs),
        ),
    )
    outputs = vcat(
        map(elm -> NormMassFractionLabel(atm.material, elm), elms),
        map(elm -> MassFractionLabel(atm.material, elm), elms),
    )
    afls = map(elm -> AtomicFractionLabel(atm.material, elm), elms)
    awls = map(elm -> AtomicWeightLabel(atm.material, elm), elms)
    s = sum(i -> inputs[afls[i]] * inputs[awls[i]], eachindex(elms))
    results = map(i -> (inputs[afls[i]] * inputs[awls[i]]) / s, eachindex(elms))
    results = vcat(results, results)
    jac = withJac ? zeros(Float64, length(outputs), length(inputs)) : missing
    if withJac
        off = length(elms)
        for i in eachindex(elms)
            for j in eachindex(elms)
                if i == j
                    @assert outputs[i].element == elms[j]
                    ii = indexin(inputs, afls[i])
                    jac[i+off, ii] =
                        (jac[i, ii] = results[i] * (1.0 - results[i]) / inputs[afls[i]])
                    ii = indexin(inputs, awls[i])
                    jac[i+off, ii] =
                        (jac[i, ii] = results[i] * (1.0 - results[i]) / inputs[awls[i]])
                else
                    ij = indexin(inputs, afls[j])
                    jac[i+off, ij] =
                        (jac[i, ij] = -results[i] * results[j] / inputs[afls[j]])
                    ij = indexin(inputs, awls[j])
                    jac[i+off, ij] =
                        (jac[i, ij] = -results[i] * results[j] / inputs[awls[j]])
                end
            end
        end
    end
    return (LabeledValues(outputs, results), jac)
end

"""
    MatStats

Computes the mean atomic number (MeanZ) and mean atomic weight (MeanAz).
"""
struct MatStats <: MeasurementModel
    material::String
end

struct MeanZ <: Label
    material::String
end
Base.show(io::IO, mz::MeanZ) = print(io, "MeanZ[", mz.material, "]")

struct MeanAz <: Label
    material::String
end
Base.show(io::IO, maz::MeanAz) = print(io, "MeanAz[$(maz.material)]")

struct AnalyticalTotal <: Label
    material::String
end
Base.show(io::IO, maz::AnalyticalTotal) = print(io, "Σ[$(maz.material)]")

const MatStatTypes = (MeanZ, MeanAz, AnalyticalTotal)

function NeXLUncertainties.compute(
    ma::MatStats,
    inputs::LabeledValues,
    withJac::Bool,
)::MMResult
    # Extract the elements for which there is atomic-fraction data associated with ma.material
    elms = map(
        il -> il.element,
        filter(
            l -> (l isa MassFractionLabel) && isequal(ma.material, l.material),
            labels(inputs),
        ),
    )
    outputs = [MeanZ(ma.material), MeanAz(ma.material), AnalyticalTotal(ma.material)]
    mfls = map(elm -> MassFractionLabel(ma.material, elm), elms)
    awls = map(elm -> AtomicWeightLabel(ma.material, elm), elms)
    results = [
        sum(inputs[mfls[i]] * elms[i].number for i in eachindex(elms)),
        sum(inputs[mfls[i]] * inputs[awls[i]] for i in eachindex(elms)),
        sum(inputs[mfls[i]] for i in eachindex(elms)),
    ]
    jac = withJac ? zeros(Float64, length(outputs), length(inputs)) : missing
    if withJac
        for j in eachindex(elms)
            jac[1, indexin(inputs, mfls[j])] = elms[j].number
            jac[2, indexin(inputs, mfls[j])] = inputs[awls[j]]
            jac[2, indexin(inputs, awls[j])] = inputs[mfls[j]]
            jac[3, indexin(inputs, mfls[j])] = 1.0
        end
    end
    return (LabeledValues(outputs, results), jac)
end

"""
    mf2comp(material::String, mfs::UncertainValues)::UncertainValues
    mf2comp(mat::Material)::UncertainValues

Converts a material composition expressed in the `mfs` UncertainValues struct into a handful
of common representations including normalized mass fraction, atomic fraction, mean Z and
mean atomic number.  The second form converts a Material into UncertainValues form.
"""
mf2comp(material::String, mfs::UncertainValues)::UncertainValues =
    propagate(AllInputs() | MFtoAF(material) | MFtoNMF(material) | MatStats(material), mfs)

function NeXLCore.mf2comp(mat::Material)::UncertainValues
    dlu = Dict{Label,UncertainValue}()
    for elm in keys(mat)
        dlu[MassFractionLabel(name(mat), elm)] = mat[elm]
        dlu[AtomicWeightLabel(name(mat), elm)] = a(elm, mat)
    end
    return mf2comp(name(mat), uvs(dlu))
end

af2comp(material::String, afs::UncertainValues)::UncertainValues =
    propagate((MatStats(material) | AllInputs()) ∘ (AFtoNMF(material) | AllInputs()), afs)

function af2comp(material::String, stoic::Dict{Element,Int})::UncertainValues
    lbls, vals = Vector{Label}(), Vector{Float64}()
    for (elm, n) in stoic
        push!(lbls, AtomicFractionLabel(material, elm))
        push!(vals, convert(Float64, n))
        push!(lbls, AtomicWeightLabel(material, elm))
        push!(vals, a(elm))
    end
    vals[1:2:end] /= sum(vals[1:2:end])
    return af2comp(material, uvs(lbls, vals, zeros(length(vals))))
end

"""
    ElementByDifference

Computes one element as the different between a sum of 1.0 mass fraction.  If the sum of the other element's
mass fraction is already 1.0 or larger, returns zero.
"""
struct ElementByDifference <: MeasurementModel
    material::String
    element::Element
end

function NeXLUncertainties.compute(
    ebd::ElementByDifference,
    inputs::LabeledValues,
    withJac::Bool,
)::MMResult
    # Extract the elements for which there is mass-fraction data associated with ebd.material
    mfls = filter(
        l -> (l isa MassFractionLabel) && isequal(ebd.material, l.material),
        labels(inputs),
    )
    outputs = [MassFractionLabel(ebd.material, ebd.element)]
    s = sum(mf -> inputs[mf] for mf in mfls)
    results = [s < 1.0 ? 1.0 - s : 0.0]
    jac = withJac ? zeros(Float64, 1, length(inputs)) : missing
    if withJac
        if s < 1.0
            foreach(mfl -> jac[1, indexin(inputs, mfl)] = -1.0, mfls)
        end
    end
    return (LabeledValues(outputs, results), jac)
end

"""
    ElementByStoichiometry

Computes the mass-fraction of an element which is related by stoichiometric rules (valence-rules) to the other
elements in the material.
"""
struct ElementByStoichiometry <: MeasurementModel
    material::String
    element::Element
    valence::Any

    ElementByStoichiometry(mat::String, elm::Element = n"O", valences = valence) =
        new(mat, elm, valences)
end

function NeXLUncertainties.compute(
    ebs::ElementByStoichiometry,
    inputs::LabeledValues,
    withJac::Bool,
)::MMResult
    v(elm) = ebs.valence[elm.number]
    # Extract the elements for which there is mass-fraction data associated with ebd.material
    elms = map(
        il -> il.element,
        filter(
            l -> (l isa MassFractionLabel) && isequal(ebs.material, l.material),
            labels(inputs),
        ),
    )
    mfls = map(elm -> MassFractionLabel(ebs.material, elm), elms)
    awls = map(elm -> AtomicWeightLabel(ebs.material, elm), elms)
    outputs = [MassFractionLabel(ebd.material, ebs.element)]
    awe = AtomicWeightLabel(ebs.material, ebs.element)
    aN, vN = inputs[awe], v(ebs.element)
    results = [
        (-aN / vN) *
        sum((inputs[mfls[i]] * v(elms[i])) / inputs[awls[i]] for i in eachindex(elms)),
    ]
    jac = withJac ? zeros(Float64, 1, length(inputs)) : missing
    if withJac
        jac[1, indexin(inputs, awe)] = results[1] / aN
        for j in eachindex(elms)
            aj, vj = inputs[awls[j]], v(elms[j])
            jac[1, indexin(inputs, mfls[i])] = -(aN * vj) / (aj * vN)
            jac[1, indexin(inputs, awls[i])] = (aN * vj) / (vN * Aj^2) * inputs[mfls[j]]
        end
    end
    return (LabeledValues(outputs, results), jac)
end

struct MaterialMixture <: MeasurementModel
    material::String
end

function NeXLUncertainties.compute(
    mm::MaterialMixture,
    inputs::LabeledValues,
    withJac::Bool,
)::MMResult
    isconstitmf(mat, lbl) = (lbl isa MassFractionLabel) && isequal(mat, lbl.material)
    isconstitmf(mat, lbl, elm) = isconstitmf(mat, lbl) && (lbl.element == elm)
    isconstit(mat, lbl) = (lbl isa MaterialFractionLabel) && isequal(mat, lbl.material)
    # Extract the material names for which there is material-fraction data associated with mm.material
    constitlbls = filter(lbl -> isconstit(mm.material, lbl), labels(inputs))
    resmf, den = Dict{Element,Float64}(), Dict{Element,Float64}()
    for constitlbl in constitlbls
        # How much of this constituent???
        f = inputs[constitlbl]
        # Convert this into element mass fractions
        for mfl in filter(lbl -> isconstitmf(constitlbl.constituent, lbl), labels(inputs))
            @assert mfl isa MassFractionLabel
            resmf[mfl.element] = get(resmf, mfl.element, 0.0) + f * inputs[mfl]
            awl = AtomicWeightLabel(constitlbl.constituent, mfl.element)
            den[mfl.element] = get(den, mfl.element, 0.0) + f * inputs[mfl] / inputs[awl]
        end
    end
    elms = sort(collect(keys(resmf))) # Put them in z order for convenience...
    # Construct `outputs` and `results`
    outputs, results =
        Array{Label}(undef, 2 * length(elms)), zeros(Float64, 2 * length(elms))
    jac = withJac ? zeros(Float64, 2 * length(elms), length(inputs)) : missing
    for (i, elm) in enumerate(elms)
        mfi, awi = 2i - 1, 2i
        outputs[mfi] = MassFractionLabel(mm.material, elm)
        results[mfi] = resmf[elm]
        outputs[awi] = AtomicWeightLabel(mm.material, elm)
        aMz = (results[awi] = resmf[elm] / den[elm]) # Eqn 23
        if withJac
            for constitlbl in constitlbls
                Mjz = inputs[constitlbl]
                mf = MassFractionLabel(constitlbl.constituent, elm)
                if haskey(inputs, mf) # Does this constituent have this element?
                    awl = AtomicWeightLabel(mf.material, elm)
                    cjz, ajz, cMz = inputs[mf], inputs[awl], resmf[elm]
                    jac[mfi, indexin(inputs, mf)] = Mjz # Ok
                    jac[mfi, indexin(inputs, constitlbl)] = cjz # Ok
                    jac[awi, indexin(inputs, constitlbl)] =
                        (cjz / cMz) * (aMz - aMz^2 / ajz) # 23.1
                    jac[awi, indexin(inputs, mf)] = (Mjz / cMz) * (aMz - aMz^2 / ajz)
                    jac[awi, indexin(inputs, awl)] = (cMz * (aMz^2)) / ajz
                end
            end
        end
    end
    return (LabeledValues(outputs, results), jac)
end

"""
    mixture(mat::String, mix::Pair{UncertainValues,UncertainValue}...)::UncertainValues
    mixture(mat::String, mix::Pair{Material, UncertainValue}...)::UncertainValues

Converts a mixture of materials into a material.
In the first case, the input materials are in UncertainValues representation.
In the second case, the Material instances are converted to the UncertainValues representation.
In each case, the first item in the Pair represents a material and the second the mass fraction of that material in the mixture.
"""
function mixture(mat::String, mix::Pair{Material,UncertainValue}...)
    mixture(mat, (mf2comp(mat) => uv for (mat, uv) in mix)...)
end

function mixture(mat::String, mix::Pair{UncertainValues,UncertainValue}...)
    function nameofmat(uvs)
        lbls = labels(uvs)
        lbl = lbls[findfirst(lbl -> lbl isa MassFractionLabel, labels(uvs))]
        @assert all(l -> l.material == lbl.material, lbls)
        return lbl.material
    end
    mixes = uvs(Dict(MaterialFractionLabel(mat, nameofmat(uvs)) => uv for (uvs, uv) in mix))
    return propagate(
        AllInputs() | MaterialMixture(mat),
        cat(mixes, (uvs for (uvs, uv) in mix)...),
    )
end

struct μoρElementLabel <: Label
    element::Element
    xray::CharXRay
end
Base.show(io::IO, l::μoρElementLabel) = print(io, "[μ/ρ][$(l.element),$(l.xray)]")

struct μoρLabel <: MaterialLabel
    material::String
    xray::CharXRay
end
Base.show(io::IO, l::μoρLabel) = print(io, "[μ/ρ][$(l.material),$(l.xray)]")

struct μoρMaterial <: MeasurementModel
    material::String
    xray::CharXRay
end

function NeXLUncertainties.compute(
    mm::μoρMaterial,
    inputs::LabeledValues,
    withJac::Bool,
)::MMResult
    mfls = materiallabels(MassFractionLabel, mm.material, inputs)
    emacls = [μoρElementLabel(mfl.element, mm.xray) for mfl in mfls]
    emacs, mfs = [inputs[emacl] for emacl in emacls], [inputs[mfl] for mfl in mfls]
    vals = LabeledValues(
        [μoρLabel(mm.material, mm.xray)],
        [sum(emacs[i] * mfs[i] for i in eachindex(mfls))],
    )
    jac = withJac ? zeros(Float64, length(vals), length(inputs)) : missing
    if withJac
        for i in eachindex(mfls)
            jac[1, indexin(inputs, mfls[i])] = emacs[i]
            jac[1, indexin(inputs, emacls[i])] = mfs[i]
        end
    end
    return (vals, jac)
end


text:
struct EvolutionDescription
  numberOfRates::Int
  rateTimes::Vector{Float64}
  evolutionTimes::Vector{Float64}
  relevanceRates::Vector{Pair{Int, Int}}
  rateTaus::Vector{Float64}
  firstAliveRate::Vector{Int}
end

# Constructors #
EvolutionDescription() = EvolutionDescription(-1, Vector{Float64}(), Vector{Float64}(), Vector{Pair{Int, Int}}(), Vector{Float64}(), Vector{Int}())

function EvolutionDescription(rateTimes::Vector{Float64},
                              evolutionTimes::Vector{Float64},
                              relevanceRates::Vector{Pair{Int, Int}} = Vector{Pair{Int, Int}}())

  numberOfRates = isempty(rateTimes) ? 0 : length(rateTimes) - 1
  evolutionTimes = isempty(evolutionTimes) && ~isempty(rateTimes) ? rateTimes[1:end-1] : evolutionTimes

  firstAliveRate = zeros(length(evolutionTimes))

  rateTaus = check_increasing_times_and_calculate_taus(rateTimes)
  check_increasing_times(evolutionTimes)
  numberOfSteps = length(evolutionTimes)

  evolutionTimes[end] <= rateTimes[end - 1] || error("The last evolution time is past the last fixing time")

  if isempty(relevanceRates)
    relevanceRates = fill(Pair(1, numberOfRates+1), numberOfSteps)
  else
    length(relevanceRates) == numberOfSteps || error("relevanceRates / evolutionTimes mismatch")
  end

  currentEvolutionTime = 0.0
  firstAliveRateNum = 1
  @inbounds @simd for j in eachindex(firstAliveRate)
    while rateTimes[firstAliveRateNum] <= currentEvolutionTime
      firstAliveRateNum += 1
    end

    firstAliveRate[j] = firstAliveRateNum
    currentEvolutionTime = evolutionTimes[j]
  end

  return EvolutionDescription(numberOfRates, rateTimes, evolutionTimes, relevanceRates, rateTaus, firstAliveRate)
end

number_of_steps(evol::EvolutionDescription) = length(evol.evolutionTimes)

clone(ed::EvolutionDescription) = EvolutionDescription(ed.numberOfRates, copy(ed.rateTimes), copy(ed.evolutionTimes), copy(ed.relevanceRates),
                                  copy(ed.rateTaus), copy(ed.firstAliveRate))

## Misc functions ##
function money_market_plus_measure(ev::EvolutionDescription, offset::Int)
  rateTimes = ev.rateTimes
  maxNumeraire = length(rateTimes)

  offset <= maxNumeraire || error("offset is greater than the max allowed value for numeraire")
  evolutionTimes = ev.evolutionTimes
  n = length(evolutionTimes)
  numeraires = Vector{Int}(undef, n)
  j = 1
  @inbounds @simd for i = 1:n
    while rateTimes[j] < evolutionTimes[i]
      j += 1
    end
    numeraires[i] = min(j + offset, maxNumeraire)
  end

  return numeraires
end

money_market_measure(evol::EvolutionDescription) = money_market_plus_measure(evol, 0)

function check_compatibility(evol::EvolutionDescription, numeraires::Vector{Int})
  evolutionTimes = evol.evolutionTimes
  n = length(evolutionTimes)

  length(numeraires) == n || error("size mismatch between numeraires and evolution times")

  rateTimes = evol.rateTimes
  @inbounds @simd for i = 1:n-1
    rateTimes[numeraires[i]] >= evolutionTimes[i] || error("$(i+1) step, evolution time $(evolutionTimes[i]): the numeraire ($(numeraires[i]))
                                                        corresponding to the rate time $(rateTimes[numeraires[i]]) is expired")
  end
end


text:
# SPDX-License-Identifier: X11
# 2020-11-20
# Knapsack 1 (100pt)

function main()
  N,W = parse.(Int, (split ∘ readline)())
  dp = zeros(Int, N + 1, W + 1)

  for i ∈ 2:N + 1
    w,v = parse.(Int, (split ∘ readline)())

    for j ∈ 2:W + 1
      if j <= w
        dp[i, j] = dp[i - 1, j]
      else
        dp[i, j] = max(dp[i - 1, j],
                       dp[i - 1, j - w] + v)
      end
    end
  end

  println(dp[N + 1, W + 1])
end

main()

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
module projectile
using ReachabilityAnalysis, SparseArrays

## system matrix
A = sparse([1, 3], [2, 4], [0.5, 0.7], 4, 4)

## affine term
b = sparsevec([4], [-9.81], 4)

function model(X0)
    S = @system(x' = Ax + b)
    return IVP(S, X0)
end
end # module

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
using BookBayes
using LightGraphs

C = Variable(:c, 3);
S = Variable(:s, 3);
V = Variable(:v, 3);

# P(C)
φC = Factor([C], FactorTable(
    (c=1,) => 0.80,
    (c=2,) => 0.19,
    (c=3,) => 0.01,
    ));
# P(S|C)
φSC = Factor([S,C], FactorTable(
    (s=1,c=1,) => 0.001,
    (s=1,c=2,) => 0.200,
    (s=1,c=3,) => 0.800,
    (s=2,c=1,) => 0.009,
    (s=2,c=2,) => 0.750,
    (s=2,c=3,) => 0.199,
    (s=3,c=1,) => 0.990,
    (s=3,c=2,) => 0.050,
    (s=3,c=3,) => 0.001,
    ));

# P(V|C)
φVC = Factor([V,C], FactorTable(
    (v=1,c=1,) => 0.2,
    (v=1,c=2,) => 0.5,
    (v=1,c=3,) => 0.4,
    (v=2,c=1,) => 0.2,
    (v=2,c=2,) => 0.4,
    (v=2,c=3,) => 0.4,
    (v=3,c=1,) => 0.6,
    (v=3,c=2,) => 0.1,
    (v=3,c=3,) => 0.2,
    ));

g = SimpleDiGraph(3)

add_edge!(g, 1, 2)
add_edge!(g, 1, 3);

bn = BayesianNetwork([C,S,V], [φC, φSC, φVC], g);

m_weighted = LikelihoodWeightedSampling(100000);
m_var_elim = VariableElimination(3:1);

r_exact= infer(ExactInference(), bn, (:c,), (s=2, v=1))
r_weighted = infer(m_weighted, bn, (:c,), (s=2,v=1,))
r_var_elim = infer(m_var_elim, bn, (:c,), (s=2,v=1,))


text:
#
# Stuff to read FRS and HBAI files into frames and merge them
# a lot of this is just to overcome the weird way I have these stored
# and would need adaption by anyone with a sensible layout who had all the
# files in the one format
#
using DataFrames
using StatFiles
using IterableTables
using IteratorInterfaceExtensions
using TableTraits
using CSV

global HBAI_DIR = "/mnt/data/hbai/UKDA-5828-stata11/stata11/"

global HBAIS = Dict(
        2015=>"hbai1516_g4.dta",
        2014=>"hbai1415_g4.dta",
        2013=>"hbai1314_g4.dta",
        2012=>"hbai1213_g4.dta",
        2011=>"hbai1112_g4.dta",
        2010=>"hbai1011_g4.dta",
        2009=>"hbai0910_g4.dta",
        2008=>"hbai0809_g4.dta",
        2007=>"hbai0708_g4.dta",
        2006=>"hbai0607_g4.dta",
        2005=>"hbai0506_g4.dta",
        2004=>"hbai0405_g4.dta",
        2003=>"hbai0304_g4.dta"
)

global FRS_DIR = "/mnt/data/frs/"
global DIRS = Dict(
        2015=>"/2015/UKDA-8171-stata11_se/stata11_se/",
        2014=>"/2014/UKDA-8013-stata11/stata11/",
        2013=>"/2013/UKDA-7753-stata11/stata11/",
        2012=>"/2012/UKDA-7556-stata11/stata11/",
        2011=>"/2011/UKDA-7368-stata9/stata9/",
        2010=>"/2010/UKDA-7085-spss/spss/spss14/",
        2009=>"/2009/UKDA-6886-spss/spss/spss14/",
        2008=>"/2008/tab/",
        2007=>"/2007/tab/",
        2006=>"/2006/tab/",
        2005=>"/2005/tab/",
        2004=>"/2004/tab/",
        2003=>"/2003/tab/"
)

function extn(year)
        if( year < 2011) && ( year > 2008)
                return ".sav"
        elseif( year <= 2008 )
                return ".tab"
        end
        ".dta"
end

function dfname( fname, year)
        string( FRS_DIR, DIRS[year],fname,extn(year))
end

function remapRegion( r :: Integer )
        if r == 112000001 # North East
                r = 1
        elseif r == 112000002 # North West
                r = 2
        # no region 3
        elseif r == 112000003 # Yorks and the Humber
                r = 4
        elseif r == 112000004 # East Midlands
                r = 5
        elseif r == 112000005 # West Midlands
                r = 6
        elseif r == 112000006
                r = 7  # East of England
        elseif r == 112000007 # | London
                r = 8
        elseif r == 112000008 # South East
                r = 9
        elseif r == 112000009 # South West
                r = 10
        elseif r == 399999999 #  wales
                r = 11
        elseif r == 299999999
                r = 12        # Scotland
        elseif r == 499999999 # Northern Ireland
                r = 13
        end
        return r
end

function initialiseMerged( n )

# .. example check
# select value,count(value),label from dictionaries.enums where dataset='frs' and tables='househol' and variable_name='hhcomps' group by value,label;

        merged = DataFrame(
                frs_year    = Vector{Union{Int64,Missing}}(missing,n),
                sernum  = Vector{Union{Int64,Missing}}(missing,n),
                benunit = Vector{Union{Int8,Missing}}(missing,n),
                person  = Vector{Union{Int8,Missing}}(missing,n),
                interview_year  = Vector{Union{Int64,Missing}}(missing,n),
                interview_month  = Vector{Union{Int8,Missing}}(missing,n),
                gvtregn  = Vector{Union{Int8,Missing}}(missing,n), # f reg change 2012
                tenure  = Vector{Union{Int8,Missing}}(missing,n), # f enums OK
                hhld_income = Vector{Union{Real,Missing}}(missing,n),
                gbhscost = Vector{Union{Real,Missing}}(missing,n),
                hhcomps = Vector{Union{Int8,Missing}}(missing,n), # full - check for differences (OK)
                num_children= Vector{Union{Int8,Missing}}(missing,n),
                num_adults = Vector{Union{Int8,Missing}}(missing,n),
                age_in_bands = Vector{Union{Int8,Missing}}(missing,n), # age80 - 2005-2012 ; iagegr2 full
                empstati = Vector{Union{Int8,Missing}}(missing,n),
                sex = Vector{Union{Int8,Missing}}(missing,n),     # full
                marital = Vector{Union{Int8,Missing}}(missing,n), # full
                ethgrp  = Vector{Union{Int8,Missing}}(missing,n), #
                ftwk =  Vector{Union{Int64,Missing}}(missing,n), #
                tea = Vector{Union{Int64,Missing}}(missing,n), #
                nssec = Vector{Union{Real,Missing}}(missing,n), # FIXME Really a Decimal
                hi2qual = Vector{Union{Int8,Missing}}(missing,n), # 2008-11
                sic = Vector{Union{Int64,Missing}}(missing,n),
                soc2010 = Vector{Union{Int64,Missing}}(missing,n),
                spcreg1 = Vector{Union{Int8,Missing}}(missing,n), # full years
                spcreg2 = Vector{Union{Int8,Missing}}(missing,n),
                spcreg3 = Vector{Union{Int8,Missing}}(missing,n),
                disdif1 = Vector{Union{Int8,Missing}}(missing,n), # 2003-11 only for these
                disdif2 = Vector{Union{Int8,Missing}}(missing,n),
                disdif3 = Vector{Union{Int8,Missing}}(missing,n),
                disdif4 = Vector{Union{Int8,Missing}}(missing,n),
                disdif5 = Vector{Union{Int8,Missing}}(missing,n),
                disdif6 = Vector{Union{Int8,Missing}}(missing,n),
                disdif7 = Vector{Union{Int8,Missing}}(missing,n),
                disdif8 = Vector{Union{Int8,Missing}}(missing,n),
                disdifp1 = Vector{Union{Int8,Missing}}(missing,n),
                hourcare = Vector{Union{Real,Missing}}(missing,n), # full

                #hourcb = Vector{Union{Real,Missing}}(missing,n), # full
                #hourch = Vector{Union{Real,Missing}}(missing,n), # full
                #hourcl = Vector{Union{Real,Missing}}(missing,n), # full
                #hourfr = Vector{Union{Real,Missing}}(missing,n), # full
                #hourot = Vector{Union{Real,Missing}}(missing,n), # f
                #hourre = Vector{Union{Real,Missing}}(missing,n), # f

                hourtot = Vector{Union{Real,Missing}}(missing,n), # f
                incseo2  = Vector{Union{Real,Missing}}(missing,n), # f
                tothours  = Vector{Union{Real,Missing}}(missing,n), # f
                indinc   = Vector{Union{Real,Missing}}(missing,n), # f
                indisben  = Vector{Union{Real,Missing}}(missing,n), # f
                inearns   = Vector{Union{Real,Missing}}(missing,n), # f
                ininv    = Vector{Union{Real,Missing}}(missing,n),  # f
                inirben   = Vector{Union{Real,Missing}}(missing,n), # f
                innirben  = Vector{Union{Real,Missing}}(missing,n), # f
                inothben  = Vector{Union{Real,Missing}}(missing,n),  # f
                inpeninc  = Vector{Union{Real,Missing}}(missing,n), # f
                inrinc   = Vector{Union{Real,Missing}}(missing,n), # f
                jobtype  = Vector{Union{Int8,Missing}}(missing,n),  # f - job
                etype = Vector{Union{Int8,Missing}}(missing,n), # 07-08 missing - job
                earnings = Vector{Union{Real,Missing}}(missing,n), # 07-08 missing
                usual_hours = Vector{Union{Real,Missing}}(missing,n), # 07-08 missing
                actual_hours = Vector{Union{Real,Missing}}(missing,n), # 07-08 missing
                hourly_wage = Vector{Union{Real,Missing}}(missing,n),
                jobsect = Vector{Union{Int8,Missing}}(missing,n), # 2010-2015
                pencont = Vector{Union{Real,Missing}}(missing,n), # f
                pays_union_contrib  = Vector{Union{Int8,Missing}}(missing,n), # 2010-2015
                cpi = Vector{Union{Real,Missing}}(missing,n),
                cpih = Vector{Union{Real,Missing}}(missing,n),
                gdpdeflator = Vector{Union{Real,Missing}}(missing,n) # f
                )
        merged
end

function mergeAll( year, hhld, adult, job, penprov, hbai, prices, gdpdef )
        @assert isiterabletable( hhld )
        @assert isiterabletable( adult )
        @assert isiterabletable( job )
        @assert isiterabletable( penprov )
        @assert isiterabletable( hbai )
        # another approach is ::AbstractDataFrame
        # .. and so on
        numadults = size(adult)[1]
        merged = initialiseMerged( numadults )
        adno = 0
        for ao in 1:numadults
                ad   = adult[ao,:]
                asn  = ad.sernum
                abu  = ad.benunit
                apn  = ad.person

                hh   = hhld[( hhld.sernum .== asn ),:]

                hba  = hbai[(( hbai.sernum .== asn ) .&
                             ( hbai.benunit .== abu )), :]

                nhh = size( hh )[1]
                nhba = size( hba )[1]
                @assert nhba <= 1
                @assert nhh  == 1
                gvtregn = remapRegion(hh[1,:gvtregn])
                if( nhba == 1) && ( nhh == 1 ) && ( gvtregn != 13 )  # only those with HBAI record, not in NI
                        adno += 1
                        ajob = job[(( job.sernum .== asn ) .&
                                    ( job.benunit .== abu ) .&
                                    ( job.person .== apn )),:]

                        apen = penprov[(( penprov.sernum .== asn ) .&
                                        ( penprov.benunit .== abu ) .&
                                        ( penprov.person .== apn )),:]
                        njobs = size( ajob )[1]
                        npens = size( apen )[1]
                        nhba  = size( hba )[1]
                        dd = split( hh[1,:intdate], "/" )
                        interview_year = parse( Int64, dd[3])
                        interview_month  = parse( Int8, dd[1])
                        quarter = div( interview_month-1, 3 ) + 1
                        merged[adno,:interview_year] = interview_year
                        merged[adno,:interview_month] = interview_month
                        # print( "interview_month $interview_month quarter $quarter\n" )

                        aprice = prices[((prices.year .== interview_year) .&
                                         (prices.month .== interview_month )), : ]
                        agdp = gdpdef[((gdpdef.year .== interview_year) .&
                                (gdpdef.q .== quarter )), : ]

                        merged[adno,:frs_year] = year
                        merged[adno,:sernum] = asn
                        merged[adno,:benunit] = abu
                        merged[adno,:person] = apn
                        merged[adno,:gvtregn] = gvtregn
                        merged[adno,:tenure] = hh[1,:tenure]
                        merged[adno,:hhcomps] = hh[1,:hhcomps]
                        merged[adno,:num_children] = hba[1,:depchldh]
                        merged[adno,:num_adults] = hba[1,:adulth]
                        if adult[ao,:ftwk] > 0
                                merged[adno,:ftwk] = adult[ao,:ftwk]
                        end
                        if( year >= 2005 )
                                merged[adno,:age_in_bands] = adult[ao,:age80]
                        else
                                iagegr2 = adult[ao,:iagegr2]
                                age = 0
                                if iagegr2 == 4
                                   age = 20 # Age 16 to 24
                                elseif iagegr2 == 5
                                   age = 30 # Age 25 to 34
                                elseif iagegr2 == 6
                                   age = 40 # Age 35 to 44
                                elseif iagegr2 == 7
                                   age = 50 # Age 45 to 54
                                elseif iagegr2 == 8
                                   age = 58 # Age 55 to 59
                                elseif iagegr2 == 9
                                   age = 62 # Age 60 to 64
                                elseif iagegr2 == 10
                                   age = 70 # Age 65 to 74
                                elseif iagegr2 == 11
                                   age = 80 # Age 75 to 84
                                elseif iagegr2 == 12
                                   age = 85 # Age 85 or over
                                else
                                        @assert "unmapped iagegr2=$iagegr2"
                                end # age
                                @assert( age >= 16 ) && (age <= 85)
                                merged[adno,:age_in_bands] = age
                        end
                        merged[adno,:sex] = adult[ao,:sex]
                        merged[adno,:marital] = adult[ao,:marital]
                        merged[adno,:tea] = max(0,adult[ao,:tea])
                        merged[adno,:nssec] = max(0,adult[ao,:nssec])
                        merged[adno,:empstati] = adult[ao,:empstati]

                        ethgr3 = -1
                        ethgrp = -1;
                        if year < 2011
                                ethgrp = adult[ao,:ethgrp]
                                if ethgrp in ( 1,2,3 )
                                        ethgr3 = 1 # white
                                elseif ethgrp in ( 4,5,6,7 )
                                        ethgr3 = 2 # mixed
                                elseif ethgrp in ( 8,9,10,11,15 )
                                        ethgr3 = 3 # asian inc chinese
                                elseif ethgrp in ( 12,13,14 )
                                        ethgr3 = 4 # black
                                elseif ethgrp in ( 16 )
                                        ethgr3 = 5 # other
                                end
                        else
                                ethgr3 = adult[ao,:ethgr3]
                        end
                        @assert (ethgr3 in ( 1:5 )) "failed to map $ethgr3 for person $ao ethgrp $ethgrp year $year"
                        merged[adno,:ethgrp]  = ethgr3 # merged white,oth,asian,black,oth

                        if year >= 2008 && year <= 2011
                                merged[adno,:hi2qual] = adult[ao,:hi2qual] # 2008-11
                        end # fixme merge in others ..
                        merged[adno,:sic] = adult[ao,:sic] # recoded in 2009 Vector{Union{Int8,Missing}}(missing,n),
                        if year  >= 2011
                                merged[adno,:soc2010] = adult[ao,:soc2010] # 2011-15
                        end
                        if adult[ao,:spcreg1] > 0
                                merged[adno,:spcreg1] = adult[ao,:spcreg1]
                        end
                        if adult[ao,:spcreg2] > 0
                                merged[adno,:spcreg2] = adult[ao,:spcreg2]
                        end
                        if adult[ao,:spcreg3] > 0
                                merged[adno,:spcreg3] = adult[ao,:spcreg3]
                        end
                        if year >= 2003 && year <= 2011
                                merged[adno,:disdif1] = adult[ao,:disdif1] == 1 ? 1 : 0
                                merged[adno,:disdif2] = adult[ao,:disdif2] == 1 ? 1 : 0
                                merged[adno,:disdif3] = adult[ao,:disdif3] == 1 ? 1 : 0
                                merged[adno,:disdif4] = adult[ao,:disdif4] == 1 ? 1 : 0
                                merged[adno,:disdif5] = adult[ao,:disdif5] == 1 ? 1 : 0
                                merged[adno,:disdif6] = adult[ao,:disdif6] == 1 ? 1 : 0
                                merged[adno,:disdif7] = adult[ao,:disdif7] == 1 ? 1 : 0
                                merged[adno,:disdif8] = adult[ao,:disdif8] == 1 ? 1 : 0
                        elseif year > 2011
                                merged[adno,:disdifp1] = adult[ao,:disdifp1]  == 1 ? 1 : 0
                        end
                        merged[adno,:hourcare] = max(0,adult[ao,:hourcare])
                        ht = adult[ao,:hourtot]
                        hourtot = 0.0
                        ht = adult[ao,:hourtot]
                        if ht ==  0
                           hourtot = 0 # 0 hours per week
                        elseif ht ==  1
                           hourtot = 2 # 0-4 hours per week
                        elseif ht ==  2
                           hourtot = 7 # 5-9 hours per week
                        elseif ht ==  3
                           hourtot = 15 # 10-19 hours per week
                        elseif ht ==  4
                           hourtot = 28 # 20-34 hours per week
                        elseif ht ==  5
                           hourtot = 42 # 35-49 hours per week
                        elseif ht ==  6
                           hourtot = 75 # 50-99 hours per week
                        elseif ht ==  7
                           hourtot = 100 # 100 or more hours per week
                        elseif ht ==  8
                           hourtot = 10 # Varies - under 20 hours per week
                        elseif ht ==  9
                           hourtot = 27 # Varies - 20-34 hours per week
                        elseif ht ==  10
                           hourtot = 50 # Varies - 35 hours a week or more
                        end
                        merged[adno,:hourtot] = hourtot

                        merged[adno,:incseo2 ] = max(0,adult[ao,:incseo2 ])
                        merged[adno,:tothours ] = max(0,adult[ao,:tothours ])
                        merged[adno,:indinc  ] = adult[ao,:indinc  ]
                        merged[adno,:indisben ] = adult[ao,:indisben ]
                        merged[adno,:inearns  ] = adult[ao,:inearns  ]
                        merged[adno,:ininv ] = adult[ao,:ininv   ]
                        merged[adno,:inirben  ] = adult[ao,:inirben  ]
                        merged[adno,:innirben ] = adult[ao,:innirben ]
                        merged[adno,:inothben ] = adult[ao,:inothben ]
                        merged[adno,:inpeninc ] = adult[ao,:inpeninc ]
                        merged[adno,:inrinc  ] = adult[ao,:inrinc  ]
                        penamt = 0.0
                        for p in 1:npens
                                if( apen[p,:penamt] > 0 )
                                        if apen[p,:penamtpd]==95
                                                penamt += apen[p,:penamt]/52.0
                                        else
                                                penamt += apen[p,:penamt]
                                        end
                                end
                        end
                        earnings = 0.0
                        actual_hours = 0.0
                        usual_hours = 0.0
                        pays_union_contrib = false
                        for j in 1:njobs
                                if j == 1 # take 1st record job for all of these
                                        merged[adno,:etype] =  ajob[j,:etype]
                                        merged[adno,:jobtype] =  ajob[j,:jobtype]
                                        if( year >= 2010 )
                                                merged[adno,:jobsect] = ajob[j,:jobsect]
                                        end
                                end
                                if ajob[j,:jobhours] > 0
                                        actual_hours += ajob[j,:jobhours]
                                end
                                if ! (year in (2007,2008))
                                        if ajob[j,:dvushr] > 0
                                                usual_hours += ajob[j,:dvushr]
                                        end
                                        if( year >= 2012)
                                                if( ajob[j,:othded03] == 1)
                                                        pays_union_contrib = true
                                                end
                                        else
                                                if( ajob[j,:othded3] == 1)
                                                        pays_union_contrib = true
                                                end
                                        end
                                end

                                addBonus = false
                                if ajob[j,:ugross] > 0.0 # take usual when last not usual
                                        earnings += ajob[j,:ugross]
                                        addBonus = true
                                elseif ajob[j,:grwage] > 0.0 # then take last
                                        earnings += ajob[j,:grwage]
                                        addBonus = true
                                elseif ajob[j,:ugrspay] > 0.0 # then take total pay, but don't add bonuses
                                        earnings += ajob[j,:ugrspay]
                                end
                                if addBonus
                                        for i in 1:6
                                                bon = Symbol( string("bonamt",i))
                                                tax = Symbol( string("bontax",i))
                                                if ajob[j,bon] > 0.0
                                                        bon = ajob[j,bon]
                                                        if  ajob[j,tax] == 2
                                                                bon /= (1-0.22) # fixme hack basic rate
                                                        end
                                                        earnings += bon/52.0 # fixwme weeks per year
                                                end
                                        end # bonuss loop
                                end # add bonuses
                        end # jobs loop
                        @assert earnings >= 0.0
                        @assert actual_hours >= 0.0
                        @assert usual_hours >= 0.0
                        merged[adno,:earnings] = earnings
                        merged[adno,:actual_hours] = actual_hours
                        merged[adno,:usual_hours] = usual_hours
                        if pays_union_contrib
                                merged[adno,:pays_union_contrib] = 1
                        else
                                merged[adno,:pays_union_contrib] = 0
                        end
                        if actual_hours > 0
                                merged[adno,:hourly_wage] = earnings/actual_hours # inc bonuses
                        end
                        merged[adno,:pencont] = penamt
                        merged[adno,:hhld_income] = hba[1,:egrinchh] #FRS extended - gross income for the household
                        merged[adno,:gbhscost] = hba[1,:ehcost]

                        merged[adno,:cpi] = aprice[1,:D7BT]
                        merged[adno,:cpih] = aprice[1,:L522]
                        merged[adno,:gdpdeflator] = agdp[1,:L8GG]
                end # add this person - is in HBAI
                # print( "$asn $abu $apn $njobs $npens \n")
        end # loop round Adult
        merged[1:adno,:]
end

function loadone( which, year )
        fname = string( "/mnt/data/frs/", year, "/tab/", which, ".tab")
        loadtoframe( fname )
end

# parse.(Int64,split.(h95[:intdate],"/")[:][1])

global MONTHS = Dict(
        "JAN"=>1,
        "FEB"=>2,
        "MAR"=>3,
        "APR"=>4,
        "MAY"=>5,
        "JUN"=>6,
        "JUL"=>7,
        "AUG"=>8,
        "SEP"=>9,
        "OCT"=>10,
        "NOV"=>11,
        "DEC"=>12 )

function loadGDPDeflator( name )
        prices = DataFrame(CSV.File( name, separator=',' ))
        np = size(prices)[1]
        prices[:year] = zeros(Int64,np) #Union{Int64,Missing},np)
        prices[:q] = zeros(Int8,np) #zeros(Union{Int64,Missing},np)
        dp = r"([0-9]{4}) Q([1-4])"
        for i in 1:np
                rc = match( dp, prices[i,:CDID] )
                if( rc !== nothing )
                        prices[i,:year] = parse(Int64,rc[1])
                        prices[i,:q] = parse( Int8,rc[2])
                #else
        #                prices[i,:year] = 0 # missing
        #                prices[i,:month] = 0 # missing
                end
        end
        prices
end

function loadPrices( name )
        prices = DataFrame(CSV.File( name, separator=',' ))
        np = size(prices)[1]
        prices[:year] = zeros(Int64,np) #Union{Int64,Missing},np)
        prices[:month] = zeros(Int8,np) #zeros(Union{Int64,Missing},np)
        dp = r"([0-9]{4}) ([A-Z]{3})"
        for i in 1:np
                rc = match( dp, prices[i,:CDID] )
                if( rc !== nothing )
                        prices[i,:year] = parse(Int64,rc[1])
                        prices[i,:month] = MONTHS[rc[2]]
                #else
        #                prices[i,:year] = 0 # missing
        #                prices[i,:month] = 0 # missing
                end
        end
        prices
end

"""
Load one frame using [] and jam all names to lower case (FRS changes case by year)
"""
function loadtoframe( name )
        (pname,ext) = splitext( name )
        if ext == ".tab"
                df = CSV.File( name, delim='\t') |> DataFrame
        else
                df = DataFrame(load(name))
        end
        # all names as lc strings
        lcnames = Symbol.(lowercase.(string.(names(df))))
        names!(df,lcnames)
        df
end

function loadoneold( name, year )
        name = dfname( name, year )
        return loadtoframe( name )
end

function makeAll( startyear, endyear )
        merged = initialiseMerged(0)
        prices = loadPrices( "/mnt/data/prices/mm23/mm23_edited.csv" );
        gdpdef = loadGDPDeflator( "/mnt/data/prices/gdpdef.csv" )
        for year in startyear:endyear
                print( "on year $year " )
                hh = loadone("househol", year);
                ad = loadone("adult", year);
                pen = loadone("penprov", year );
                job = loadone("job", year );
                hbai = loadtoframe( string( HBAI_DIR, HBAIS[year]))
                ds = mergeAll( year, hh, ad, job, pen, hbai, prices, gdpdef )
                append!( merged, ds )
        end
        merged
end

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
module RocketErrorIfEmptyOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: error_if_empty()" begin

    println("Testing: operator error_if_empty()")

    run_proxyshowcheck("ErrorIfEmpty", error_if_empty(0))

    run_testset([
        (
            source      = from(1:5) |> error_if_empty("Empty"),
            values      = @ts([ 1:5, c ]),
            source_type = Int
        ),
        (
            source      = completed(Int) |> error_if_empty("Empty"),
            values      = @ts(e("Empty")),
            source_type = Int
        ),
        (
            source      = faulted(Int, "e") |> error_if_empty("Empty"),
            values      = @ts(e("e")),
            source_type = Int
        ),
        (
            source      = never(Int) |> error_if_empty("Empty"),
            values      = @ts(),
            source_type = Int
        ),
        (
            source      = completed(Int) |> error_if_empty("Empty") |> catch_error((d, obs) -> of(1)),
            values      = @ts([ 1, c ]),
            source_type = Int
        )
    ])

end

end


text:
#=
  https://dtai.cs.kuleuven.be/problog/tutorial/basic/03_dice.html
  """
  The following example illustrates the use of a logic program with recursion and lists.
  We start by rolling the first die. Every roll determines the next die to roll, but we stop
  if we have used that die before. We query for the possible sequences of rolled dice.
  We use three-sided dice instead of the regular six-sided ones simply to restrict the number
  of possible outcomes (and thus inference time).
  """

  Cf ~/blog/rolling_dice5.blog
     ~/webppl/rolling_dice5.wppl
=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function rolling_dice5(n=3)

    function roll(a)
        # t ~ DiscreteUniform(1,n)
        t = rand(DiscreteUniform(1,n))
        if t in a
            # We have used this die before -> end this sequence
            # Note: In the WebPPL model the duplicate value is included
            # return vcat(a,t)
            return a
        else
            return roll(vcat(a,t))
        end
    end

    a = roll([])
    len ~ DiscreteUniform(1,n)
    len = length(a)
    # println("a:$a len:$len")

    # In the example, the first roll is 1.
    true ~ Dirac(a[1] == 1)

    return len
    # return a

end

function run_model(len=0)
    model = rolling_dice5(len)

    num_chains = 4

    # chains = sample(model, Prior(), 10_000)

    # chains = sample(model, MH(), MCMCThreads(), 100_000, num_chains)
    # chains = sample(model, MH(), MCMCThreads(), 10_000, num_chains)
    # chains = sample(model, MH(), 10_000)
    chains = sample(model, MH(), 1_000)

    # chains = sample(model, PG(15), MCMCThreads(), 1_000, num_chains)

    # chains = sample(model, SMC(1000), MCMCThreads(), 10_000, num_chains)
    # chains = sample(model, SMC(1000), 10_000)
    # chains = sample(model, IS(), 10_000)
    #
    # display(chains)
    # show_var_dist_pct(chains,:len,1000)

    println("prob return value:")
    genq = generated_quantities(model, chains)
    show_var_dist_pct(genq,1000)
    println("mean: $(mean(genq))")
    return genq
end

genq = undef
for val in 1:30
    global genq
    println("\nval:$val")
    @time genq = run_model(val)
end

println("\nval=10")
println("mean:", mean(run_model(10)))

println("\nval=100")
println("mean:", mean(run_model(100)))

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

module Meshes

using Tables
using StaticArrays
using CircularArrays
using RecipesBase
using LinearAlgebra
using Random

using IterTools: ivec
using StatsBase: Weights
using SpecialFunctions: gamma
using Distances: Euclidean, Mahalanobis
using ReferenceFrameRotations: angle_to_dcm
using NearestNeighbors: KDTree, BallTree, knn, inrange

import Tables
import Random
import Base: values, ==, +, -
import StatsBase: sample
import NearestNeighbors: MinkowskiMetric

# numerical tolerances
include("tolerances.jl")

# basic types
include("vectors.jl")
include("points.jl")
include("angles.jl")

# type traits
include("traits/domain.jl")
include("traits/data.jl")

# point sets
include("pointsets.jl")

# geometries and meshes
include("geometries.jl")
include("connectivities.jl")
include("conventions.jl")
include("mesh.jl")

# domain and data views
include("views.jl")

# neighborhoods and searches
include("neighborhoods.jl")
include("neighborsearch.jl")

# partitions
include("partitions.jl")

# algorithms
include("sampling.jl")
include("partitioning.jl")
include("discretization.jl")
include("boundingboxes.jl")

# utilities
include("utils.jl")

# plot recipes
include("plotrecipes/domain.jl")
include("plotrecipes/data.jl")
include("plotrecipes/geometries.jl")
include("plotrecipes/pointsets.jl")
include("plotrecipes/cartesiangrids.jl")
include("plotrecipes/mesh.jl")

export
  # points
  Point, Point1, Point2, Point3, Point1f, Point2f, Point3f,
  embeddim, coordtype, coordinates,

  # vectors
  Vec, Vec1, Vec2, Vec3, Vec1f, Vec2f, Vec3f,

  # angles
  ∠,

  # traits
  Domain, Data,
  domain, values, asarray,

  # geometries
  Geometry,
  embeddim, paramdim, coordtype,
  measure, area, volume, boundary,
  centroid,

  # primitives
  Primitive,
  Box, Ball, Sphere, Cylinder,
  center, radius, height, sides,
  measure, diagonal,

  # polytopes
  Polytope, Polygon, Polyhedron,
  Segment, Triangle, Quadrangle,
  Pyramid, Tetrahedron, Hexahedron,
  Chain, PolyArea,
  vertices, nvertices,
  faces, facets,
  intersecttype,
  windingnumber, chains,
  isclosed, issimple, hasholes,
  angles, innerangles, close!, open!,
  orientation, bridge,

  # connectivities
  Connectivity,
  polytopetype, connect, materialize,

  # ordering conventions
  OrderingConvention,
  GMSH,
  connectivities,

  # rotation conventions
  RotationConvention,
  GSLIB, Leapfrog, Datamine,
  TaitBryanExtr, TaitBryanIntr,
  EulerExtr, EulerIntr,
  axesseq, orientation, angleunits,
  mainaxis, isextrinsic, rotmat,

  # point sets
  PointSet,
  coordinates, coordinates!,
  nelements,

  # meshes
  Mesh,
  CartesianGrid, UnstructuredMesh,
  faces, elements, spacing,
  coordinates, coordinates!,
  nelements,

  # views
  DomainView, DataView,

  # neighborhoods
  NormBall, Ellipsoid,
  metric, radius,

  # neighbordhood search
  NeighborSearchMethod,
  BoundedNeighborSearchMethod,
  NeighborhoodSearch,
  KNearestSearch,
  KBallSearch,
  BoundedSearch,
  search!, search,

  # partitions
  Partition,
  subsets, metadata,

  # sampling
  SamplingMethod,
  RegularSampling,
  UniformSampling,
  WeightedSampling,
  BallSampling,
  sample,

  # partitioning
  PartitionMethod,
  PredicatePartitionMethod,
  SPredicatePartitionMethod,
  RandomPartition,
  FractionPartition,
  BlockPartition,
  BisectPointPartition,
  partition,

  # discretization
  DiscretizationMethod,
  FIST,
  discretize,

  # bounding boxes
  boundingbox,

  # utililities
  signarea,
  sideof

end # module


text:
using IntervalConstraintProgramming
using ValidatedNumerics

using FactCheck

facts("Utilities") do
    @fact IntervalConstraintProgramming.unify_variables([:a, :c], [:c, :b]) -->
        ([:a,:b,:c], [1,3], [3,2], [1,0,2], [0,2,1])
end

facts("Separators") do
    II = -100..100
    X = IntervalBox(II, II)
    S = @constraint x^2 + y^2 <= 1

    inner, outer = S(X)
    @fact inner --> (-1..1, -1..1)
    @fact outer --> (II, II)

    X = IntervalBox(-∞..∞, -∞..∞)
    inner, outer = S(X)
    @fact inner --> (-1..1, -1..1)
    @fact outer --> (-∞..∞, -∞..∞)
end

facts("pave") do
    S1a = @constraint x > 0
    S1b = @constraint y > 0

    S1 = S1a ∩ S1b
    paving = pave(S1, IntervalBox(-3..3, -3..3), 0.1)
    @fact paving.inner --> [IntervalBox(1.5..3, 0..3), IntervalBox(0..1.5, 0..3)]
    @fact isempty(paving.boundary) --> true

    S2 = S1a ∪ S1b
    paving = pave(S2, IntervalBox(-3..3, -3..3), 0.1)
    @fact paving.inner --> [IntervalBox(-3..0, 0..3), IntervalBox(0..3, -3..3)]
    @fact isempty(paving.boundary) --> true

    S3 = @constraint x^2 + y^2 <= 1
    X = IntervalBox(-Inf..Inf, -Inf..Inf)
    paving = pave(S3, X, 1)

    @fact paving.inner --> [IntervalBox(Interval(0.0, 0.5), Interval(0.0, 0.8660254037844386)),
                    IntervalBox(Interval(0.0, 0.5), Interval(-0.8660254037844386, 0.0)),
                    IntervalBox(Interval(-0.5, 0.0), Interval(0.0, 0.8660254037844386)),
                    IntervalBox(Interval(-0.5, 0.0), Interval(-0.8660254037844386, 0.0))]

    @fact paving.boundary --> [ IntervalBox(Interval(0.5, 1.0), Interval(0.0, 0.8660254037844387)),
                        IntervalBox(Interval(0.0, 0.5), Interval(0.8660254037844386, 1.0)),
                        IntervalBox(Interval(0.5, 1.0), Interval(-0.8660254037844387, 0.0)),
                        IntervalBox(Interval(0.0, 0.5), Interval(-1.0, -0.8660254037844386)),
                        IntervalBox(Interval(-0.5, 0.0), Interval(0.8660254037844386, 1.0)),
                        IntervalBox(Interval(-1.0, -0.5), Interval(0.0, 0.8660254037844387)),
                        IntervalBox(Interval(-0.5, 0.0), Interval(-1.0, -0.8660254037844386)),
                        IntervalBox(Interval(-1.0, -0.5), Interval(-0.8660254037844387, 0.0))]
end

facts("Constants") do
    a = 3
    x = y = -∞..∞
    X = IntervalBox(x, y)
    S4 = @constraint x^2 + y^2 <= $a
    paving = pave(S4, X)

    @fact paving.ϵ --> 0.01
    @fact length(paving.inner) --> 1532
    @fact length(paving.boundary) --> 1536
end

facts("Volume") do
    x = 3..5
    @fact Vol(x).bounds --> 2

    V = Vol(IntervalBox(-1..1.5, 2..3.5))
    @fact typeof(V) --> IntervalConstraintProgramming.Vol{2, Float64}
    @fact V.bounds --> 3.75

end


text:
include(normpath(joinpath(@__DIR__, "IOM", "src", "models", "EMOM", "EMOM.jl")))
using NCDatasets
using MPI
using Formatting

MPI.Init()

println("Processing data...")

using TOML
config = TOML.parsefile("data/config.toml")

init_POP_file = "hist/paper2021_POP2_CTL.pop.h.daily.0002-01-01.nc"
domain_file = config["MODEL_CORE"]["domain_file"]

init_file  = config["MODEL_MISC"]["init_file"]
topo_file  = config["MODEL_CORE"]["topo_file"]

Nz = length(config["MODEL_CORE"]["z_w"]) - 1 # Layers used. Thickness ≈ 503m

Dataset(domain_file, "r") do ds
    global Nx = ds.dim["ni"]
    global Ny = ds.dim["nj"]
    global mask = ds["mask"][:]
end

println("Initial profile uses the first record of $(init_POP_file)")
Dataset(init_POP_file, "r") do ds
    global TEMP  = permutedims(nomissing(ds["TEMP"][:, :, 1:Nz, 1],  NaN), [3, 1, 2])
    global SALT  = permutedims(nomissing(ds["SALT"][:, :, 1:Nz, 1],  NaN), [3, 1, 2])
end

Nx, Ny = size(TEMP)[2:3]
Nz_bot = zeros(Int64, Nx, Ny)
mask_T = zeros(Float64, Nz, Nx, Ny)
mask_T[isfinite.(TEMP)] .= 1.0
valid_idx = mask_T .== 1.0
Nz_bot .= sum(mask_T, dims=1)[1, :, :]

println(format("Output file: {}", topo_file))

Dataset(topo_file, "c") do ds

    defDim(ds, "Nx", Nx)
    defDim(ds, "Ny", Ny)
    defVar(ds, "Nz_bot", Nz_bot, ("Nx", "Ny", ), ; attrib = Dict(
        "long_name" => "z-grid idx of the deepest cell",
        "units"     => "none",
    ))

end

println("Create a model to save initial condition.")

ev = EMOM.Env(config["MODEL_CORE"])
mb = EMOM.ModelBlock(ev; init_core=false)

mb.fi.sv[:TEMP][valid_idx] .= TEMP[valid_idx]
mb.fi.sv[:SALT][valid_idx] .= SALT[valid_idx]

println(format("Output file: {}", init_file))

EMOM.takeSnapshot(DateTimeNoLeap(1,1,1), mb, init_file)

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM