text:
################################################################################
#  Copyright 2020, Marta Vanin, Tom Van Acker                                  #
################################################################################
# PowerModelsDistributionStateEstimation.jl                                    #
# An extention package of PowerModelsDistribution.jl for Static Power System   #
# State Estimation.                                                            #
################################################################################

## CSV to measurement parser
function dataString_to_array(input::String)::Array
    if size(split(input, ","))[1] > 1
        rm_brackets = input[2:(end-1)]
        raw_string = split(rm_brackets, ",")
        float_array = [parse(Float64, ss) for ss in raw_string]
    else
        float_array = [parse(Float64, input)]
    end
    return float_array
end
"""
    sample_fake_measurements(meas_row_dst, sorted_args, seed)
Called within read_measurement!. Samples measurement errors from a distribution.
"""
function sample_fake_measurements(meas_row_dst, sorted_args, seed)
    pkg_id = meas_row_dst == "ExtendedBeta" ? _PMDSE : _DST
    distr = [getfield(pkg_id, Symbol(meas_row_dst))(tuple(sa...)...) for sa in sorted_args]
    randRNG = [_RAN.seed!(seed+i) for i in 1:length(sorted_args)]
    fake_meas = [_RAN.rand(randRNG[i], distr[i]) for i in 1:length(sorted_args)]
    return [[fake_meas[i], last.(sorted_args)[i]] for i in 1:length(sorted_args)]
end
"""
    read_measurement!(data, meas_row; sample_fake_meas, seed)
Adds measurement data from separate CSV file to the PowerModelsDistribution data
dictionary, called from within `add_measurements!()`.
"""
function read_measurement!(data::Dict, meas_row::_DFS.DataFrameRow, sample_fake_meas, seed)
    dst_params = []
    for nr in 1:4
        if "par_$(nr)" ∈ names(meas_row) && !ismissing(meas_row[Symbol("par_",nr)]) && meas_row[Symbol("par_",nr)] != "missing"
            par_array = isa(meas_row[Symbol("par_",nr)], String) ? dataString_to_array(meas_row[Symbol("par_",nr)]) : [meas_row[Symbol("par_",nr)]]
            push!(dst_params, par_array)
        end
    end
    @assert all([length(dst_params[1])==length(dst_params[i]) for i in 1:length(dst_params)]) "The lengths of the distribution parameters must be equal, e.g, all three-phase"
    sorted_args = []
    for ph in 1:length(dst_params[1])
        sorted_arg = [param[ph] for param in dst_params]
        push!(sorted_args, sorted_arg)
    end
    if sample_fake_meas
        sorted_args = sample_fake_measurements(meas_row[:dst], sorted_args, seed)
    end
    meas_row[:dst] == "ExtendedBeta" ? pkg_id = _PMDSE : pkg_id = _DST
    data["dst"] = [getfield(pkg_id, Symbol(meas_row[:dst]))(tuple(sa...)...) for sa in sorted_args]
end
"""
    add_measurements!(data::Dict, meas_file::String; actual_meas::Bool = false, seed::Int=0)

Add measurement data from separate CSV file to the PowerModelsDistribution data
dictionary. To fully understand how this function works, it is recommended to first read
the documentation section that describes the CSV measurement file format.
# Arguments
-   `data`: MATHEMATICAL data dictionary in a format usable with PowerModelsDistribution
-   `meas_file`: path to and name of file with measurement data
-   `actual_meas`: default is false. When applied to non-normal distributions,
        the effect is overruled to that of `true`. For normal distributions, the following applies:
    *   `false`: the "par_1" in meas_file are not actual measurements, but, e.g.,
        error-free powerflow results. Then, a fake measurement is built, extracting
        a value with an error from the given normal distribution.
    *   `true`: the "par_1" values in meas_file are actual measurement values,
        and the "par_2" are the σs of the measurements' distributions. These are
        directly used as input of the state estimator without further processing.
    *   if a "parse" column is present in the CSV file, the `true` or `false` is associated
        to each individual row (i.e., measurement), and overrules whatever the actual_meas
        input of add_measurements!() itself is.
-   `seed`: random seed value to make the results reproducible and explore different
    Monte Carlo scenarios when sampling measurement with errors from a probability distribution.
"""
function add_measurements!(data::Dict, meas_file::String; actual_meas::Bool=false, seed::Int = 0)
    meas_df = _CSV.read(meas_file)
    data["meas"] = Dict{String, Any}()
    [data["meas"]["$m_id"] = Dict{String,Any}() for m_id in meas_df[!,:meas_id]]
    for row in 1:size(meas_df)[1]
        sample_fake_meas = meas_df[row,:dst] ∈ ["Normal"] ? !actual_meas : false
        data["meas"]["$(meas_df[row,:meas_id])"] = Dict{String,Any}()
        data["meas"]["$(meas_df[row,:meas_id])"]["cmp"] = Symbol(meas_df[row,:cmp_type])
        data["meas"]["$(meas_df[row,:meas_id])"]["cmp_id"] = meas_df[row,:cmp_id]
        data["meas"]["$(meas_df[row,:meas_id])"]["var"] = Symbol(meas_df[row,:meas_var])
        if "crit" ∈ names(meas_df)
            data["meas"]["$(meas_df[row,:meas_id])"]["crit"] = meas_df[row, :crit]
        end
        if "sample" ∈ names(meas_df) && !ismissing(meas_df[row, :sample]) && meas_df[row,:dst] ∈ ["Normal"]
            sample_fake_meas = meas_df[row, :sample]
        end
        read_measurement!(data["meas"]["$(meas_df[row,:meas_id])"], meas_df[row,:], sample_fake_meas, seed)
    end
end

## Powerflow results to CSV
repeated_measurement(df::_DFS.DataFrame, cmp_id::String, cmp_type::String, phases) =
    sum(.&(df[!,:cmp_id].==cmp_id,df[!,:cmp_type].==cmp_type,df[!,:phase].==string(phases))) > 0

function get_measures(model::DataType, cmp_type::String)
    if model <: _PMD.AbstractUnbalancedACPModel
        if cmp_type == "bus"  return ["vm"] end
        if cmp_type == "gen"  return ["pg","qg"] end
        if cmp_type == "load" return ["pd","qd"] end
    elseif model  <: _PMD.AbstractUnbalancedIVRModel
        # NB: IVR is a subtype of ACR, therefore it should preceed ACR
        if cmp_type == "bus"  return ["vr","vi"] end
        if cmp_type == "gen"  return ["crg","cig"] end
        if cmp_type == "load" return ["crd_bus","cid_bus"] end
    elseif model <: _PMD.AbstractUnbalancedACRModel
        if cmp_type == "bus"  return ["vr","vi"] end
        if cmp_type == "gen"  return ["pg","qg"] end
        if cmp_type == "load" return ["pd","qd"] end
    elseif model <: _PMD.SDPUBFPowerModel
        #if cmp_type == "bus"  return ["vr","vi"] end
        if cmp_type == "gen"  return ["pg","qg"] end
        if cmp_type == "load" return ["pd","qd"] end
    elseif model <: _PMD.LinDist3FlowModel
        if cmp_type == "bus"  return ["w"] end
        if cmp_type == "gen"  return ["pg","qg"] end
        if cmp_type == "load" return ["pd","qd"] end
    end
    return []
end
function reduce_name(meas_var::String)
    if meas_var == "crd_bus" return "crd" end
    if meas_var == "cid_bus" return "cid" end
    return meas_var
end
function get_sigma(σ::Float64, meas_var::String,phases)
    sigma = meas_var in ["vm","vr","vi"] ? σ/3 : σ/5/3 ;
    return length(phases) == 1 ? sigma : sigma.*ones(length(phases)) ;
end
get_configuration(cmp_type::String, cmp_data::Dict{String,Any}) = "G"
init_measurements() =
    _DFS.DataFrame(meas_id=Int64[], cmp_type=String[], cmp_id=String[],
                         meas_type=String[], meas_var=String[], phase=String[],
                         dst=String[], par_1=Union{String, Missing}[], par_2=Union{String, Missing}[],
                         par_3=Union{String, Missing}[], par_4=Union{String, Missing}[], crit=Union{String, Missing}[])
"""
   Function to write the dataframe row for a Normal component
"""
function write_cmp_measurement!(df::_DFS.DataFrame, model::Type, cmp_id::String, cmp_type::String, cmp_data::Dict{String,Any},
                                        cmp_res::Dict{String,Any}, phases, very_basic_case::Bool; exclude::Vector{String}=String[], σ::Float64)

    if length(phases) == 1
        phases = phases[1]
        ph = 1
    else
        ph = [1,2,3]
    end

    if !repeated_measurement(df, cmp_id, cmp_type, phases)
        config = get_configuration(cmp_type, cmp_data)
        for meas_var in get_measures(model, cmp_type) if !(meas_var in exclude)
            par_2 = very_basic_case ? string(get_sigma(σ, meas_var,phases)) : string(length(phases) == 1 ? σ : σ.*ones(length(phases)))
            push!(df, [length(df.meas_id)+1,                  # meas_id
                       cmp_type,                              # cmp_type
                       cmp_id,                                # cmp_id
                       config,                                # config
                       reduce_name(meas_var),                 # meas_var
                       string(phases),                        # phase
                       "Normal",                              # dst
                       string(cmp_res[meas_var][ph]),         # par_1
                       par_2,                                 # par_2
                       missing, missing, missing             # par_3,4,crit
                      ])
    end end end
end
function write_cmp_measurements!(df::_DFS.DataFrame, model::Type, cmp_type::String,
                                 data::Dict{String,Any}, pf_results::Dict{String,Any}, very_basic_case::Bool;
                                 exclude::Vector{String}=String[], σ::Union{Float64, Dict})
    for (cmp_id, cmp_res) in pf_results["solution"][cmp_type]
        # write the properties for the component
        cmp_data = data[cmp_type][cmp_id]
        phases = cmp_data["connections"]
        σ_cmp = isa(σ, Dict) ? σ[cmp_type] : σ
        write_cmp_measurement!(df, model, cmp_id, cmp_type, cmp_data, cmp_res, phases, very_basic_case, exclude = exclude, σ = σ_cmp)
        # write the properties for its bus
        cmp_id = string(cmp_data["$(cmp_type)_bus"])
        cmp_data = data["bus"][cmp_id]
        cmp_res = pf_results["solution"]["bus"][cmp_id]
        σ_cmp = isa(σ, Dict) ? σ["bus"] : σ
        write_cmp_measurement!(df, model, cmp_id, "bus", cmp_data, cmp_res, cmp_data["terminals"], very_basic_case, exclude = exclude, σ = σ_cmp)
    end
end
"""
    write_measurements!(model::Type, data::Dict, pf_results::Dict, path::String; exclude::Vector{String}=String[], σ::Float64=0.005)

Function to write a csv file with measurements, to be used to run state estimation calculations. The
file is built starting from power flow results from PowerModelsDistribution.jl (or any dictionary with
the same format). The measurements consist of voltage and power/current injections in correspondence
of all generators and loads. The exact measurement type depends on the chosen power flow formulation,
e.g., with the AC Polar formulation, these are voltage magnitude and active and reactive power.
# Arguments
-   `model`: power flow type of the generated measurements, e.g., ACPUPowerModel.
             If it does not match the power flow model of the `pf_results`, it might not work.
             `pf_results` can be post-processed, e.g., polar results can be converted in rectangular
             and viceversa, to make the result dictionary compatible.
-   `data`: MATHEMATICAL data dictionary in a format usable with PowerModelsDistribution
-   `pf_results`: PowerModelsDistribution solution dictionary or similar format
-   `path`: path where the csv file will be generated and stored
-   `exclude`: select quantities from the `pf_results` dictionary to be excluded from the measurement
               generation. For example, to ignore generator results with ACPUPowerModel,
               set exclude = ["pg", "qg"].
-   `σ`: standard deviation of demand/generation measurements. If a float, their stdev of the relative
        bus voltage measurements is taken with `get_sigma()`. If a dictionary (recommended option), the
        individual σ of each measurement quantity is declared.
"""
function write_measurements!(model::Type, data::Dict, pf_results::Dict, path::String; exclude::Vector{String}=String[], σ::Union{Dict, Float64}=0.005)
    df = init_measurements()
    for cmp_type in ["gen", "load"]
        σ_cmp = isa(σ, Dict) ? σ[cmp_type] : σ
        very_basic_case = isa(σ, Dict) ? false : true
        write_cmp_measurements!(df, model, cmp_type, data, pf_results, very_basic_case, exclude = exclude, σ = σ_cmp)
    end
    _CSV.write(path, df)
end
"""
    add_voltage_measurement!(model::Type, data::Dict, pf_results::Dict, path::String)

This function can be run after add_measurements! for the cases in which only power and/or
    current measurements are generated.
"""
function add_voltage_measurement!(data::Dict, pf_result::Dict, sigma::Float64)

    first_key =  first(keys(pf_result["solution"]["bus"]))
    if !(haskey(pf_result["solution"]["bus"][first_key], "vm") || haskey(pf_result["solution"]["bus"][first_key], "w"))
        vm = sqrt.(pf_result["solution"]["bus"][first_key]["vi"].^2+pf_result["solution"]["bus"][first_key]["vr"].^2)
        voltage_meas_idx = string(maximum(parse(Int64,i) for i in keys(data["meas"]) ) + 1)
        bus_idx = parse(Int64, first_key)
        data["meas"][voltage_meas_idx] = Dict{String, Any}("var"=>:vm,"cmp"=>:bus,
                                        "dst"=>[Distributions.Normal{Float64}(v_ph, sigma) for v_ph in vm],
                                        "cmp_id"=>bus_idx)
    end
end
"""
    write_measurements_and_pseudo!(model::Type, data::Dict, pf_results::Dict, path::String; exclude::Vector{String}=String[], distribution_info::String, σ::Float64=0.005)

Helper function to write a csv file with a combination of measurements and pseudo-measurements.
Works similarly to write_measurements!() with additional support for non-Normal distributions
for the pseudo-measurements. In order to use this function, the load data for pseudo measurements
data["load"] need to point to an external file where information on the probability distribution is stored.
An example of such a file is distr_example.csv in test/extra/measurements.
The arguments of the function are the same as write_measurements!(), with the addition of
`distribution_info`: the path to external csv file for pseudo-measurements distributions.
"""
function write_measurements_and_pseudo!(model::Type, data::Dict, pf_results::Dict, path::String; exclude::Vector{String}=String[], distribution_info::String, σ::Float64=0.005)

    df = _DFS.DataFrame(meas_id=Int64[], cmp_type=String[], cmp_id=String[],meas_type=String[], meas_var=String[],
             phase=String[], dst=String[], par_1=Union{String, Missing}[], par_2=Union{String, Missing}[],
             par_3=Union{String, Missing}[], par_4=Union{String, Missing}[], crit=Union{String, Missing}[] )
    di_df = _CSV.read(distribution_info)
    very_basic_case = true # TO DO: make it an argument?
    for cmp_type in ["gen", "load"]
        for (cmp_id, cmp) in data[cmp_type]
            phases = data[cmp_type][cmp_id]["connections"]
            if haskey(cmp,"pseudo")
                write_cmp_pseudo!(df, model, cmp_type, cmp_id, data, di_df, phases, exclude)
            else
                write_cmp_measurement!(df, model, cmp_id, cmp_type, data[cmp_type][cmp_id], pf_results["solution"][cmp_type][cmp_id],
                                                    phases, very_basic_case, exclude = exclude, σ = σ)
                #if this is an actual measurement, also add voltage measurement
                bus_id = string(data[cmp_type][cmp_id]["$(cmp_type)_bus"])
                phases = data["bus"][bus_id]["terminals"]
                if !repeated_measurement(df, bus_id, "bus", phases)
                    write_cmp_measurement!(df, model, bus_id, "bus", data["bus"][bus_id], pf_results["solution"]["bus"][bus_id],
                                                    phases, very_basic_case,exclude = exclude, σ = σ)
                end
            end
        end
    end
    _CSV.write(path, df)
end
"""
    write_cmp_pseudo!(df::_DFS.DataFrame, model::Type, cmp_type::String, cmp_id::String, data::Dict{String,Any}, distr_info::_DFS.DataFrame, phases)
Used within write_measurements_and_pseudo!() to generate pseudo-measurements
"""
function write_cmp_pseudo!(df::_DFS.DataFrame, model::Type, cmp_type::String, cmp_id::String,
                                 data::Dict{String,Any}, distr_info::_DFS.DataFrame, phases, exclude)
    cosϕ = distr_info[1, :PF]
    if !repeated_measurement(df, cmp_id, cmp_type, phases)
        config = get_configuration(cmp_type, data[cmp_type][cmp_id])
        row_idx = find_row(distr_info, data[cmp_type][cmp_id]["pseudo"]["time_step"], data[cmp_type][cmp_id]["pseudo"]["cluster"],data[cmp_type][cmp_id]["pseudo"]["day"])
        for meas_var in get_measures(model, cmp_type) if !(meas_var in exclude)
            if distr_info[row_idx, :distr][1] == "ExtendedBeta" && !distr_info[row_idx, :per_unit][1]
                per_unit_div = data["settings"]["sbase_default"]
            else
                per_unit_div = 1
            end
            occursin("q", string(meas_var)) ? per_unit_div/=tan(acos(cosϕ)) : per_unit_div/=1
            #push active power pseudo
            push!(df, [length(df.meas_id)+1,              # meas_id
                   cmp_type,                              # cmp_type
                   cmp_id,                                # cmp_id
                   config,                                # meas_id
                   reduce_name(meas_var),                 # meas_var
                   string(phases),                        # phase
                   distr_info[row_idx, :distr][1],           # dst
                   string(has_par(distr_info, row_idx, 1)[1]),       # par_1
                   string(has_par(distr_info, row_idx, 2)[1]),       # par_2
                   string(has_par(distr_info, row_idx, 3)[1]/per_unit_div),  # par_3
                   string(has_par(distr_info, row_idx, 4)[1]/per_unit_div),  # par_4
                   "mle"                                  # crit
                  ])
    end end end
end
function find_row(df::DataFrames.DataFrame, time_step::Int64, cluster::Int64, day::Int64)
    bitarray_sum = ( df[!, :day].==day) + (df[!, :time_step].==time_step) + (df[!, :cluster].==cluster)
    row_index = findall(x->x==maximum(bitarray_sum), bitarray_sum)
    @assert length(row_index) == 1 "There is more than one entry that matches time step and cluster"
    return row_index
end
"""
checks if load/gen data has a given distribution parameter
"""
has_par(df, row, par_num) =
      "par_$(par_num)" ∈ names(df) ? df[row, Symbol("par_$(par_num)") ] : missing


text:
x = 12
println(typeof(x))
y = convert(UInt8, x)
println(y)
println(typeof(y))

# This will flag an error because Uint8 holds up to 2^8 - 1 = 255
# convert(UInt8, 256)

struct Point{T<:Real}
    x::T
    y::T
end

Base.convert(::Type{Point{T}}, x::Array{T, 1}) where {T<:Real} = Point(x...)
println(convert(Point{Int64}, [1, 2]))

a, b = 1.5, 2
println("$(typeof(a)), $(typeof(b))")
c, d = promote(a, b)
println("$(typeof(c)), $(typeof(d))")
# To add functionality for promote, add definition for promote_rule.
# Eg. promote_rule(::Type{Float64}, ::Type{Int32}) = Float64

ext: jl


text:
include("action_selectors/abstract_action_selector.jl")
include("agents/abstract_agent.jl")
include("approximators/abstract_approximator.jl")
include("buffers/abstract_turn_buffer.jl")
include("environments/abstract_environment.jl")
include("environment_models/abstract_environment_model.jl")
include("learners/abstract_learner.jl")
include("policies/abstract_policy.jl")
include("spaces/abstract_space.jl")


text:
using OhMyREPL

timenow = Dates.today()
welcome = "Yes, master! Today is $timenow, and my version is $VERSION."

println(welcome)

max_stars_repo_name: moldray/dotfiles

max_issues_repo_name: moldray/dotfiles

max_forks_repo_name: moldray/dotfiles


text:
# This file contains the function that manipulate the genotype files.
"""
    function pack_G(src, out; g_majored = true, skip = 0)
Pack the genotypes in `src` into the packed genotype format.

## Source file
I presume that the source file `src` are in the following format.
That is each line starts with an ID name, either of an animal or a SNP,
and then the genotypes, separated with space. And nothing else.
The number of header lines to be skipped is defined by `skip`.
The lines started with `#`s are also ignored. `skip` has priority.

## Packed genotypes
The binary file has a 24-byte header, which contains 3 Int64 number.
The first Int64 defines whether the genotypes are genotype-majored.
If `1`, then the genotypes of one ID are stored continuously
in the binary file `out`.  Otherwise, it is ID majored.  The genotypes
of one locus are stored continuously.  These 4 bytes can store other
information later.

The next 16-byte store the dimensions of the genotype matrix.
If `g_majored`, then it is `Nloci`×`Nid`.  Otherwise it is
`Nid`×`Nloci`.  For GBLUP, it is better to save `g-majored` to
speed up calculation.  For SNP-BLUP, a transposed storage is
more suitable.  The rest of the file store the genotypes of Int8.

Note that Julia matrix is column majored.
"""
function pack_G(src, out; g_majored = true, skip = 0)
    m, n, gdm = 0, 0, zeros(Int64, 3)
    open(out, "w") do io
        write(io, gdm)
        open(src, "r") do gf
            [readline(gf) for _ in 1:skip]
            for line in eachline(src)
                line[1] == '#' && continue
                s = parse.(Int8, split(line)[2:end])
                n += 1          # number of lines
                m += length(s)  # number of columns
                write(io, s)
            end
        end
    end            # genotypes written
    # A simple check point to see if the file is right.
    (m % n == 0) || @error "Genotypes not rectanglar"
    open(out, "r+") do io
        gdm = mmap(io, Vector{Int64}, 3)
        gdm[1], gdm[2], gdm[3] = g_majored, m/n, n # auto to Int64
        Mmap.sync!(gdm)
    end
end

"""
    function transpose_G(src, out)
---
Transpose **G** matrix in file `src`, and save the results to file `out`.
"""
function transpose_G(src, out)
    open(out, "w") do foo
        open(src, "r") do fin
            header = zeros(Int64, 3)
            read!(fin, header)
            write(foo, 1 - header[1], header[3], header[2])
            gt = mmap(fin, Matrix{Int8}, Tuple(header[2:3]), 24)
            write(foo, gt')
        end
    end
end


text:
using Compat
using FunctionalCollections
using FactCheck

import FunctionalCollections.KVPair

facts("Persistent Maps") do

    context("KVPairs") do
        @fact KVPair(1, 1) --> (1, 1)
        @fact (1, 1) --> KVPair(1, 1)
    end

end

typealias PAM PersistentArrayMap

facts("Persistent Array Maps") do

    context("construction") do
        @fact length(PAM{Int, Int}().kvs) --> 0
        @fact length(PAM((1, 1), (2, 2)).kvs) --> 2

        @fact length(PAM((1, 1))) --> 1
        @fact length(PAM((1, 1), (2, 2))) --> 2
    end

    context("accessing") do
        m = PAM((1, "one"), (2, "two"), (3, "three"))
        @fact m[1] --> "one"
        @fact m[2] --> "two"
        @fact m[3] --> "three"

        @fact get(m, 1) --> "one"
        @fact get(m, 1, "foo") --> "one"
        @fact get(m, 4, "foo") --> "foo"
    end

    context("haskey") do
        m = PAM((1, "one"), (2, "two"), (3, "three"))
        @fact haskey(m, 1) --> true
        @fact haskey(m, 2) --> true
        @fact haskey(m, 3) --> true
        @fact haskey(m, 4) --> false
    end

    context("assoc") do
        m = PAM{Int, AbstractString}()
        @fact assoc(m, 1, "one") --> (m) -> m[1] == "one"
        @fact try m[1]; false catch e true end --> true

        m = PAM((1, "one"))
        @fact assoc(m, 1, "foo") --> (m) -> m[1] == "foo"
    end

    context("dissoc") do
        m = PAM((1, "one"))
        m = dissoc(m, 1)
        @fact try m[1]; false catch e true end --> true

        m = PAM((1, "one"), (2, "two"))
        @fact dissoc(m, 1) --> PAM((2, "two"))
    end

    context("iterating") do
        m = PAM((1, "one"), (2, "two"), (3, "three"))
        @fact [v for (k, v) in m] --> ["one", "two", "three"]
    end

    context("isempty") do
        @fact PAM{Int, Int}() --> isempty
        @fact PAM((1, "one")) --> not(isempty)
    end

    context("equality") do
        @fact PAM((1, "one")) --> PAM((1, "one"))
        @fact PAM((1, "one"), (2, "two")) --> PAM((2, "two"), (1, "one"))
        @fact isequal(PAM((1, "one")), PAM((1, "one"))) --> true

        @fact PAM((1, "one")) --> not(PAM((2, "two")))
    end

    context("kwargs construction") do
        @fact PAM(x=1, y=2, z=3) --> PAM((:x, 1), (:y, 2), (:z, 3))
    end

    context("map") do
        m = PAM((1, 1), (2, 2), (3, 3))

        @fact map((kv) -> (kv[1], kv[2]+1), m) --> PAM((1, 2), (2, 3), (3, 4))
    end
end

typealias PHM PersistentHashMap

facts("Persistent Hash Maps") do

    context("constructor") do
        hashmap = PHM{Int, Int}()
        @fact length(hashmap) --> 0
        @fact length(PHM((1, 1), (2, 2), (3, 3))) --> 3
        @fact length(PHM(x=1, y=2, z=3)) --> 3
    end

    context("equality") do
        @fact PHM{Int, Int}() --> PHM{Int, Int}()
        @fact PHM{Int, Int}() --> PHM{AbstractString, AbstractString}()

        m1 = PHM{Int, Int}()
        m2 = PHM{Int, Int}()
        @fact assoc(m1, 1, 100) --> assoc(m2, 1, 100)
        @fact assoc(m1, 1, 100) --> not(assoc(m2, 1, 200))
        @fact assoc(m1, 1, 100) --> not(assoc(m2, 2, 100))

        m3 = PHM((1, 10), (2, 20), (3, 30))
        m4 = PHM((3, 30), (2, 20), (1, 10))
        @fact m3 --> m4
        @fact m3 --> not(m1)

        @fact m3 --> (@compat Dict(1 => 10, 2 => 20, 3 => 30))
    end

    context("assoc") do
        m = PHM{Int, AbstractString}()
        @fact assoc(m, 1, "one") --> (m) -> m[1] == "one"
        @fact try m[1]; false catch e true end --> true

        m = PHM{Int, AbstractString}()
        m = assoc(m, 1, "one")
        @fact assoc(m, 1, "foo") --> (m) -> m[1] == "foo"
    end

    context("covariance") do
        m = PHM{Any, Any}()
        @fact assoc(m, "foo", "bar") --> (@compat Dict("foo" => "bar"))
    end

    context("dissoc") do
        m = PAM((1, "one"))
        m = dissoc(m, 1)
        @fact try m[1]; false catch e true end --> true

        m = PHM((1, "one"), (2, "two"))
        @fact dissoc(m, 1) --> PHM((2, "two"))
    end

    context("get") do
        m = PHM{Int, AbstractString}()
        @fact get(m, 1, "default") --> "default"
        m = assoc(m, 1, "one")
        @fact get(m, 1, "default") --> "one"
        m = assoc(m, 1, "newone")
        @fact get(m, 1, "default") --> "newone"
        @fact try get(m, 2); false catch e true end --> true
    end

    context("haskey") do
        m = PHM{Int, AbstractString}()
        @fact haskey(m, 1) --> false
        m = assoc(m, 1, "one")
        @fact haskey(m, 1) --> true
        @fact haskey(m, 2) --> false
    end

    context("map") do
        m = PHM((1, 1), (2, 2), (3, 3))
        @fact map((kv) -> (kv[1], kv[2]+1), m) --> PHM((1, 2), (2, 3), (3, 4))
    end

    context("filter") do
        @fact filter((kv) -> iseven(kv[2]), PHM((1, 1), (2, 2))) --> PHM((2, 2))
    end

    context("merge") do
        @fact merge(PHM((1, 1), (2, 2)), PHM((2, 3), (3, 4))) -->
                PHM((1,1), (2, 3), (3, 4))
    end
end


text:
using Base: Float64
import ProjAssistant
const PjAss = ProjAssistant
PjAss.@quickactivate

@time begin
    import Chemostat_EColi
    const ChE = Chemostat_EColi

    import Chemostat_Kayser2005
    const ChK = Chemostat_Kayser2005
    const KiJR = ChK.iJR904
    const Kd = ChK.KayserData

    import ChemostatPlots
    const ChP = ChemostatPlots

    import Chemostat
    import Chemostat.LP.MathProgBase
    const Ch = Chemostat
    const ChU = Ch.Utils
    const ChLP = Ch.LP
    const ChEP = Ch.MaxEntEP
    const ChT = Ch.Test
    using ProgressMeter
    using Base.Threads
    using SparseArrays

    using Serialization
    PjAss.set_verbose(false)

    using Base.Threads
    using Statistics
    using Plots
    import GR
    !isinteractive() && GR.inline("png")

    import SimTools
    import SimTools: gd_value, grad_desc!
    const ST = SimTools

end

## --------------------------------------------------------------------
# z_beta vs z_av
function beta_vs_av(model, mov_ider, mov_betas, fix_ider, fix_betas;
        # MaxEnt
        alpha = Inf,
        maxiter = 5000,
        verbose = false,
        epsconv = 1e-5,
    )

    M, N = size(model)
    mov_idx = ChU.rxnindex(model, mov_ider)
    fix_idx = ChU.rxnindex(model, fix_ider)

    p = plot(; xlabel = string(mov_ider, " beta"), ylabel = string(mov_ider, " av"))

    lk = ReentrantLock()
    @threads for fix_beta in collect(fix_betas)

        @info("At", fix_beta)

        epout = nothing
        avs = Float64[]
        betas = Float64[]
        beta_vec = zeros(N)
        for mov_beta in mov_betas

            try
                beta_vec[mov_idx] = mov_beta
                beta_vec[fix_idx] = fix_beta
                epout = ChEP.maxent_ep(model;
                    alpha, beta_vec, maxiter, verbose, epsconv,
                    solution = epout
                )
                av = ChU.av(model, epout, mov_idx)
                (isnan(av) || isinf(av)) && break

                push!(avs, av)
                push!(betas, mov_beta)
            catch err; break; end

        end

        isempty(avs) && continue
        lock(lk) do
            plot!(p, betas, avs; label = "", lw = 3)
        end
    end

    ChE.sfig(p,
        "2.0", "beta_vs_av", (;mov_ider, fix_ider), "png"
    )
end

## --------------------------------------------------------------------
# ToyModel
let
    model = ChT.toy_model();
    mov_betas = range(0.0, 1e3; length = 500)
    fix_betas = range(0.0, 1e3; length = 50)

    mov_ider = "biom"
    fix_ider = "gt"

    beta_vs_av(model, mov_ider, mov_betas, fix_ider, fix_betas)

    mov_ider = "gt"
    fix_ider = "biom"

    beta_vs_av(model, mov_ider, mov_betas, fix_ider, fix_betas)
end

## --------------------------------------------------------------------
# beta_z vs bate_ug
function make_heatmap(f::Function, model,
        ider1, ider1_betas, ider2, ider2_betas;
        alpha = Inf,
        maxiter = 5000,
        verbose = false,
        epsconv = 1e-5,
    )

    M, N = size(model)

    idx1 = ChU.rxnindex(model, ider1)
    idx2 = ChU.rxnindex(model, ider2)

    mat = zeros(length(ider1_betas), length(ider2_betas))
    @threads for (ider1_betai, ider1_beta) in collect(enumerate(ider1_betas))

        beta_vec = zeros(N)
        @info("At", ider1_beta)
        epout = nothing

        beta_vec[idx1] = ider1_beta

        for (ider2_betai, ider2_beta) in enumerate(ider2_betas)
            beta_vec[idx2] = ider2_beta
            epout = ChEP.maxent_ep(model;
                alpha, beta_vec, maxiter, verbose, epsconv,
                solution = epout
            )
            mat[ider1_betai, ider2_betai] =  f(epout, ider1_beta, ider2_beta)
        end
    end

    p = plot(;xlabel = string(ider1, " beta"), ylabel = string(ider2, " beta"))
    heatmap!(p, ider1_betas, ider2_betas, mat'; label = "")

end

## --------------------------------------------------------------------
# "z_valid", "betas_heatmap"
let
    model = ChT.toy_model();

    Ds = range(0.0, 3.0; length = 12)

    bins = 50
    ider1 = "biom"
    ider1_betas = range(-1e2, 1e3; length = bins)
    ider2 = "gt"
    ider2_betas = range(-1e3, 1e3; length = bins)

    ps = map(Ds) do D

        z_valid_th = max(D * 0.1, 0.01)

        function get_valid(epout, args...)
            ider1_av = ChU.av(model, epout, ider1)
            ider2_av = ChU.av(model, epout, ider2)
            (isnan(ider1_av) || isinf(ider1_av)) && return NaN
            (isnan(ider2_av) || isinf(ider2_av)) && return NaN

            is_valid = (abs(ider1_av - D) <= z_valid_th)
            # is_valid = (ider2_av <= cgDX)
            is_valid ? ider1_av : 0.0
        end

        p =  make_heatmap(get_valid, model, ider1, ider1_betas, ider2, ider2_betas)
        plot!(p; title = string("D: ", round(D; sigdigits = 3)))
        return p
    end

    ChE.sfig(ps,
        "2.0", "z_valid", "betas_heatmap", (;ider1, ider2), "png"
    )

end

## --------------------------------------------------------------------
# "ug_valid", "betas_heatmap"
let
    model = ChT.toy_model();

    cgDXs = range(0.0, 10.0; length = 12)

    bins = 50
    ider1 = "biom"
    ider1_betas = range(-1e2, 1e3; length = bins)
    ider2 = "gt"
    ider2_betas = range(-1e3, 1e3; length = bins)

    ps = map(cgDXs) do cgDX

        function get_valid(epout, args...)
            ider1_av = ChU.av(model, epout, ider1)
            ider2_av = ChU.av(model, epout, ider2)
            (isnan(ider1_av) || isinf(ider1_av)) && return NaN
            (isnan(ider2_av) || isinf(ider2_av)) && return NaN

            is_valid = (ider2_av <= cgDX)
            is_valid ? ider2_av : 0.0
        end

        p =  make_heatmap(get_valid, model, ider1, ider1_betas, ider2, ider2_betas)
        plot!(p; title = string("cgDX: ", round(cgDX; sigdigits = 3)))
        return p
    end

    ChE.sfig(ps,
        "2.0", "ug_valid", "betas_heatmap", (;ider1, ider2), "png"
    )

end

## --------------------------------------------------------------------
# ## --------------------------------------------------------------------
# # EColi core
# let
#     model = ChT.ecoli_core_model();
#     mov_betas = range(0.0, 1e4; length = 500)
#     fix_betas = range(0.0, 1e2; length = 15)

#     mov_ider = "Biomass_Ecoli_core_w_GAM"
#     fix_ider = "EX_glc(e)"

#     beta_vs_av(model, mov_ider, mov_betas, fix_ider, fix_betas)
# end


text:
include("FEM1D.jl")

function B_lL_DD(x::AbstractArray{<:Real, 1}, ψ::AbstractArray{<:Real, 1}, χ::AbstractArray{<:Real, 1})
    N = length(x)+1
    m = length(ψ)
    n = Int64(N/m)
    y = Array{Float64, 1}(undef, n-1)
    for i = 1:(n-2)
        y[i] = dot(x[(1+(i-1)*m):(i*m)], χ) + dot(x[(1+i*m):((i+1)*m)], ψ)
    end
    y[end] = dot(x[(1+(n-2)*m):((n-1)*m)], χ) + dot(x[(1+(n-1)*m):end], ψ[1:(end-1)])
    return y
end

function B_Ll_DD!(x::Array{<:Real, 1}, y::Array{<:Real, 1}, ψ::Array{<:Real, 1}, χ::Array{<:Real, 1})
    n = length(y)
    m = length(ψ)
    x .= 0
    for i = 1:(n-1)
        x[(1+(i-1)*m):(i*m)] += χ*y[i]
        x[(1+i*m):((i+1)*m)] += ψ*y[i]
    end
    x[(1+(n-1)*m):(n*m)] += χ*y[end]
    x[(1+n*m):end] += ψ[1:(end-1)]*y[end]
    return nothing
end

function BPX_DD!(x::Array{<:Real, 1}, ψ::Array{<:Real, 1}, χ::Array{<:Real, 1})
    y = B_lL_DD(x, ψ, χ)
    B_Ll_DD!(x, y, ψ, χ)
    return nothing
end

function B_Ll_DD(y::AbstractArray{<:Real, 1}, ψ::AbstractArray{<:Real, 1}, χ::AbstractArray{<:Real, 1})
    n = length(y)
    m = length(ψ)
    x = Array{Float64, 1}(undef, (n+1)*m-1)
    x .= 0
    for i = 1:(n-1)
        x[(1+(i-1)*m):(i*m)] += χ*y[i]
        x[(1+i*m):((i+1)*m)] += ψ*y[i]
    end
    x[(1+(n-1)*m):(n*m)] += χ*y[end]
    x[(1+n*m):end] += ψ[1:(end-1)]*y[end]
    return x
end

function BPX_DD(x::AbstractArray{<:Real}, ψ::AbstractArray{<:Real, 1}, χ::AbstractArray{<:Real, 1})
    L = Int64(log2(length(x[:, 1])+1))
    l = L - Int64(log2(length(ψ)))
    n = 2^l
    Id = spdiagm(0 => ones(n))
    S = spdiagm(-1 => ones(n-1))
    B = (kron(Id, ψ) + kron(S, χ))[1:(end-1), 1:(end-1)]
    y = B*(transpose(B)*x)
end

function BPX_DN(x::AbstractArray{<:Real}, ψ::AbstractArray{<:Real, 1}, χ::AbstractArray{<:Real, 1})
    L = Int64(log2(length(x[:, 1])))
    l = L - Int64(log2(length(ψ)))
    n = 2^l
    Id = spdiagm(0 => ones(n))
    S = spdiagm(-1 => ones(n-1))
    B = (kron(Id, ψ) + kron(S, χ))
    y = B*(transpose(B)*x)
end

function BPX_ND(x::AbstractArray{<:Real}, ψ::AbstractArray{<:Real, 1}, χ::AbstractArray{<:Real, 1})
    L = Int64(log2(length(x[:, 1])))
    l = L - Int64(log2(length(ψ)))
    n = 2^l
    Id = spdiagm(0 => ones(n))
    S = spdiagm(+1 => ones(n-1))
    B = (kron(Id, ψ) + kron(S, χ))
    y = B*(transpose(B)*x)
end

function BPX_NN(x::AbstractArray{<:Real}, ψ::AbstractArray{<:Real, 1}, χ::AbstractArray{<:Real, 1})
    L = Int64(log2(length(x[:, 1])-1))
    l = L - Int64(log2(length(ψ)))
    n = 2^l
    Id = spdiagm(0 => ones(n))
    S = spdiagm(-1 => ones(n-1))
    e1 = zeros(n); e1[1] = 1
    e2 = zeros((1, n+1)); e2[1] = 1
    B = vcat(e2, hcat(kron(e1, χ), (kron(Id, ψ) + kron(S, χ))))
    y = B*(transpose(B)*x)
end

function BPX_layer_operator(x::AbstractArray{<:Real}, ψ::AbstractArray{<:Real, 1}, χ::AbstractArray{<:Real, 1}, BCs::AbstractArray{String, 1})
    if (BCs[1] == "N") & (BCs[2] == "N")
        return BPX_NN(x, ψ, χ)
    elseif (BCs[1] == "D") & (BCs[2] == "N")
        return BPX_DN(x, ψ, χ)
    elseif (BCs[1] == "N") & (BCs[2] == "D")
        return BPX_ND(x, ψ, χ)
    elseif (BCs[1] == "D") & (BCs[2] == "D")
        return BPX_DD(x, ψ, χ)
    end
end

function BPX_DD_matrix(L::Integer, ψ::AbstractArray{<:Real, 1}, χ::AbstractArray{<:Real, 1})
    l = L - Int64(log2(length(ψ)))
    n = 2^l
    Id = spdiagm(0 => ones(n))
    S = spdiagm(-1 => ones(n-1))
    B = (kron(Id, ψ) + kron(S, χ))[1:(end-1), 1:(end-1)]
    return B*transpose(B)
end

function BPX_DN_matrix(L::Integer, ψ::AbstractArray{<:Real, 1}, χ::AbstractArray{<:Real, 1})
    l = L - Int64(log2(length(ψ)))
    n = 2^l
    Id = spdiagm(0 => ones(n))
    S = spdiagm(-1 => ones(n-1))
    B = (kron(Id, ψ) + kron(S, χ))
    y = B*transpose(B)
end

function BPX_ND_matrix(L::Integer, ψ::AbstractArray{<:Real, 1}, χ::AbstractArray{<:Real, 1})
    l = L - Int64(log2(length(ψ)))
    n = 2^l
    Id = spdiagm(0 => ones(n))
    S = spdiagm(+1 => ones(n-1))
    B = (kron(Id, ψ) + kron(S, χ))
    y = B*transpose(B)
end

function BPX_NN_matrix(L::Integer, ψ::AbstractArray{<:Real, 1}, χ::AbstractArray{<:Real, 1})
    l = L - Int64(log2(length(ψ)))
    n = 2^l
    Id = spdiagm(0 => ones(n))
    S = spdiagm(-1 => ones(n-1))
    e1 = zeros(n); e1[1] = 1
    e2 = zeros((1, n+1)); e2[1] = 1
    B = vcat(e2, hcat(kron(e1, χ), (kron(Id, ψ) + kron(S, χ))))
    y = B*transpose(B)
end

function BPX_layer_matrix(L::Integer, ψ::AbstractArray{<:Real, 1}, χ::AbstractArray{<:Real, 1}, BCs::AbstractArray{String, 1})
    if (BCs[1] == "N") & (BCs[2] == "N")
        return BPX_NN_matrix(L, ψ, χ)
    elseif (BCs[1] == "D") & (BCs[2] == "N")
        return BPX_DN_matrix(L, ψ, χ)
    elseif (BCs[1] == "N") & (BCs[2] == "D")
        return BPX_ND_matrix(L, ψ, χ)
    elseif (BCs[1] == "D") & (BCs[2] == "D")
        return BPX_DD_matrix(L, ψ, χ)
    end
end

function BPX_DD_half_matrix(L::Integer, ψ::AbstractArray{<:Real, 1}, χ::AbstractArray{<:Real, 1})
    l = L - Int64(log2(length(ψ)))
    n = 2^l
    Id = spdiagm(0 => ones(n))
    S = spdiagm(-1 => ones(n-1))
    B = (kron(Id, ψ) + kron(S, χ))[1:(end-1), 1:(end-1)]
    return B
end

function BPX_DN_half_matrix(L::Integer, ψ::AbstractArray{<:Real, 1}, χ::AbstractArray{<:Real, 1})
    l = L - Int64(log2(length(ψ)))
    n = 2^l
    Id = spdiagm(0 => ones(n))
    S = spdiagm(-1 => ones(n-1))
    B = (kron(Id, ψ) + kron(S, χ))
    y = B
end

function BPX_ND_half_matrix(L::Integer, ψ::AbstractArray{<:Real, 1}, χ::AbstractArray{<:Real, 1})
    l = L - Int64(log2(length(ψ)))
    n = 2^l
    Id = spdiagm(0 => ones(n))
    S = spdiagm(+1 => ones(n-1))
    B = (kron(Id, ψ) + kron(S, χ))
    y = B
end

function BPX_NN_half_matrix(L::Integer, ψ::AbstractArray{<:Real, 1}, χ::AbstractArray{<:Real, 1})
    l = L - Int64(log2(length(ψ)))
    n = 2^l
    Id = spdiagm(0 => ones(n))
    S = spdiagm(-1 => ones(n-1))
    e1 = zeros(n); e1[1] = 1
    e2 = zeros((1, n+1)); e2[1] = 1
    B = vcat(e2, hcat(kron(e1, χ), (kron(Id, ψ) + kron(S, χ))))
    y = B
end

function BPX_layer_half_matrix(L::Integer, ψ::AbstractArray{<:Real, 1}, χ::AbstractArray{<:Real, 1}, BCs::AbstractArray{String, 1})
    if (BCs[1] == "N") & (BCs[2] == "N")
        return BPX_NN_half_matrix(L, ψ, χ)
    elseif (BCs[1] == "D") & (BCs[2] == "N")
        return BPX_DN_half_matrix(L, ψ, χ)
    elseif (BCs[1] == "N") & (BCs[2] == "D")
        return BPX_ND_half_matrix(L, ψ, χ)
    elseif (BCs[1] == "D") & (BCs[2] == "D")
        return BPX_DD_half_matrix(L, ψ, χ)
    end
end

function BPX_weights(L::Integer, l::Integer)
    ψ = Array{Float64, 1}(undef, 2^(L-l))
    χ = Array{Float64, 1}(undef, 2^(L-l))
    ψ .= range(1, 2^(L-l), step=1)/2^(L-l)
    χ .= 1 .- ψ
    return ψ, χ
end

function BPX_weights(L::Integer, l::Integer, BCs::AbstractArray{String, 1})
    ψ = Array{Float64, 1}(undef, 2^(L-l))
    χ = Array{Float64, 1}(undef, 2^(L-l))
    ψ .= range(1, 2^(L-l), step=1)/2^(L-l)
    if (BCs[1] == "N") & (BCs[2] == "D")
        ψ = reverse(ψ)
    end
    χ .= 1 .- ψ
    return ψ, χ
end

function BPX_DD_full(L::Integer)
    ind = []
    c = 0
    for l=1:(L-1)
        vec_size = 2^(L-l)
        push!(ind, [c+1, c+vec_size])
        push!(ind, [c+vec_size+1, c+2*vec_size])
        push!(ind, [c+2*vec_size+1, c+2*vec_size+1])
        c += 2*vec_size+1
    end
    parameters = Array{Float64, 1}(undef, 2^(L+1)+L-5+1)
    l = 1
    for k=1:3:length(ind)
        ψ, χ = BPX_weights(L, l)
        parameters[ind[k][1]:ind[k][2]] .= ψ
        parameters[ind[k+1][1]:ind[k+1][2]] .= χ
        parameters[ind[k+2][1]:ind[k+2][2]] .= 1
        l += 1
    end
    parameters[end] = 1
    B = parameters[end]*I + sum([parameters[ind[k+2][1]:ind[k+2][2]].*BPX_DD_matrix(L, parameters[ind[k][1]:ind[k][2]], parameters[ind[k+1][1]:ind[k+1][2]]) for k=1:3:length(ind)])
end

function BPX_matrix(L::Integer, BCs::AbstractArray{String, 1})
    ind = []
    c = 0
    for l=1:(L-1)
        vec_size = 2^(L-l)
        push!(ind, [c+1, c+vec_size])
        push!(ind, [c+vec_size+1, c+2*vec_size])
        push!(ind, [c+2*vec_size+1, c+2*vec_size+1])
        c += 2*vec_size+1
    end
    parameters = Array{Float64, 1}(undef, 2^(L+1)+L-5+1)
    l = 1
    for k=1:3:length(ind)
        ψ, χ = BPX_weights(L, l, BCs)
        parameters[ind[k][1]:ind[k][2]] .= ψ
        parameters[ind[k+1][1]:ind[k+1][2]] .= χ
        parameters[ind[k+2][1]:ind[k+2][2]] .= 1
        l += 1
    end
    parameters[end] = 1
    B = parameters[end]*I + sum([parameters[ind[k+2][1]:ind[k+2][2]].*BPX_layer_matrix(L, parameters[ind[k][1]:ind[k][2]], parameters[ind[k+1][1]:ind[k+1][2]], BCs) for k=1:3:length(ind)])
end

function Schwarz_matrix(A::SparseMatrixCSC{Float64, Int64}, L::Integer, BCs::AbstractArray{String, 1})
    B = spdiagm(0 => zeros(size(A)[1]))
    for l=1:(L-1)
        ψ, χ = BPX_weights(L, l, BCs)
        P = BPX_layer_half_matrix(L, ψ, χ, BCs)
        if l == 1
            A_P = sparse(inv(Array(transpose(P)*A*P) + 1e-6*I)) # spdiagm(0 => diag(transpose(P)*A*P).^-1)
        else
            A_P = spdiagm(0 => diag(transpose(P)*A*P).^-1)
        end
        B += P*A_P*transpose(P)
    end
    B += spdiagm(0 => diag(A).^-1)
    return B
end

using Gadfly, Cairo, Fontconfig, DelimitedFiles

function test_Laplace(BCs::AbstractArray{String, 1}, verbose=false)
    Ls = [5, 6, 7, 8, 9, 10, 11]
    CondA, CondAS, CondAL, CondSz = [], [], [], []
    for L = Ls
        F = [x -> -one(x), x -> zero(x), x -> zero(x), x -> x^2]
        BCs_ = [[1, BCs[1]], [2, BCs[2]]]
        A, _, _ = BVP_1D(L, F, BCs_)
        B = BPX_matrix(L, BCs)
        Sz = Schwarz_matrix(A, L, BCs)
        A_L, A_S, A_Sz = Array(B*A), Array(B*A*B), Array(Sz*A)
        singular_values_A = svd(Array(A)).S
        values_A_L = sort(abs.(eigvals(A_L)))
        values_A_S = sort(abs.(eigvals(A_S)))
        values_A_Sz = sort(abs.(eigvals(A_Sz)))
        shift = Int64((BCs[1] == "N")*(BCs[2] == "N"))
        push!(CondA, singular_values_A[1]/singular_values_A[end-shift])
        push!(CondAL, values_A_L[end]/values_A_L[1+shift])
        push!(CondAS, values_A_S[end]/values_A_S[1+shift])
        push!(CondSz, values_A_Sz[end]/values_A_Sz[1+shift])
        if verbose
            println("L = ", L)
            println("Condition number of A = ", CondA[end])
            println("Condition number of A_L = ", CondAL[end])
            println("Condition number of A_S = ", CondAS[end])
            println("Condition number of A_Sz = ", CondSz[end])
            println("=======================")
        end
    end
    # p = plot(Scale.y_log10, Guide.manual_color_key("", ["A", "left BPX", "symmetric BPX", "Schwarz"], ["red", "black", "blue", "green"]), Guide.xlabel("level"), Guide.ylabel("condition number"))
    # push!(p, layer(x=Ls, y=CondA, Geom.line, Geom.point, color=[colorant"red"]))
    # push!(p, layer(x=Ls, y=CondAL, Geom.line, Geom.point, color=[colorant"black"]))
    # push!(p, layer(x=Ls, y=CondAS, Geom.line, Geom.point, color=[colorant"blue"]))
    # push!(p, layer(x=Ls, y=CondSz, Geom.line, Geom.point, color=[colorant"green"]))
    # results = hcat(Ls, CondA, CondAS, CondAL, CondSz)
    # writedlm("precs_comparison/Laplace_"*BCs[1]*BCs[2]*".txt", results)
    # draw(SVG("precs_comparison/Laplace_"*BCs[1]*BCs[2]*".svg", 6inch, 4inch), p)
    return nothing
end

function test_Helmholtz(BCs::AbstractArray{String, 1}, k::Real, verbose=false)
    Ls = [5, 6, 7, 8, 9, 10, 11]
    CondA, CondAS, CondAL, CondSz = [], [], [], []
    for L = Ls
        F = [x -> -one(x), x -> zero(x), x -> -k^2*one(x), x -> x^2]
        BCs_ = [[1, BCs[1]], [2, BCs[2]]]
        A, _, _ = BVP_1D(L, F, BCs_)
        B = BPX_matrix(L, BCs)
        Sz = Schwarz_matrix(A, L, BCs)
        A_L, A_S, A_Sz = Array(B*A), Array(B*A*B), Array(Sz*A)
        singular_values_A = svd(Array(A)).S
        values_A_L = sort(abs.(eigvals(A_L)))
        values_A_S = sort(abs.(eigvals(A_S)))
        values_A_Sz = sort(abs.(eigvals(A_Sz)))
        shift = 0 # Int64((BCs[1] == "N")*(BCs[2] == "N"))
        push!(CondA, singular_values_A[1]/singular_values_A[end-shift])
        push!(CondAL, values_A_L[end]/values_A_L[1+shift])
        push!(CondAS, values_A_S[end]/values_A_S[1+shift])
        push!(CondSz, values_A_Sz[end]/values_A_Sz[1+shift])
        if verbose
            println("L = ", L)
            println("Condition number of A = ", CondA[end])
            println("Condition number of A_L = ", CondAL[end])
            println("Condition number of A_S = ", CondAS[end])
            println("Condition number of A_Sz = ", CondSz[end])
            println("=======================")
        end
    end
    p = plot(Scale.y_log10, Guide.manual_color_key("", ["A", "left BPX", "symmetric BPX", "Schwarz"], ["red", "black", "blue", "green"]), Guide.xlabel("level"), Guide.ylabel("condition number"))
    push!(p, layer(x=Ls, y=CondA, Geom.line, Geom.point, color=[colorant"red"]))
    push!(p, layer(x=Ls, y=CondAL, Geom.line, Geom.point, color=[colorant"black"]))
    push!(p, layer(x=Ls, y=CondAS, Geom.line, Geom.point, color=[colorant"blue"]))
    push!(p, layer(x=Ls, y=CondSz, Geom.line, Geom.point, color=[colorant"green"]))
    results = hcat(Ls, CondA, CondAS, CondAL, CondSz)
    writedlm("precs_comparison/Helmholtz_"*string(k)*"_"*BCs[1]*BCs[2]*".txt", results)
    draw(SVG("precs_comparison/Helmholtz_"*string(k)*"_"*BCs[1]*BCs[2]*".svg", 6inch, 4inch), p)
    return nothing
end

function test_Convection_Diffusion(BCs::AbstractArray{String, 1}, v::Real, verbose=false)
    Ls = [5, 6, 7, 8, 9, 10, 11]
    CondA, CondAS, CondAL, CondSz = [], [], [], []
    for L = Ls
        F = [x -> -one(x), x -> -v*one(x), x -> zero(x), x -> x^2]
        BCs_ = [[1, BCs[1]], [2, BCs[2]]]
        A, _, _ = BVP_1D(L, F, BCs_)
        B = BPX_matrix(L, BCs)
        Sz = Schwarz_matrix(A, L, BCs)
        A_L, A_S, A_Sz = Array(B*A), Array(B*A*B), Array(Sz*A)
        singular_values_A = svd(Array(A)).S
        values_A_L = sort(abs.(eigvals(A_L)))
        values_A_S = sort(abs.(eigvals(A_S)))
        values_A_Sz = sort(abs.(eigvals(A_Sz)))
        shift = Int64((BCs[1] == "N")*(BCs[2] == "N"))
        if (BCs[1] == "D")*(BCs[2] == "N")
            shift = 1
        end
        push!(CondA, singular_values_A[1]/singular_values_A[end-shift])
        push!(CondAL, values_A_L[end]/values_A_L[1+shift])
        push!(CondAS, values_A_S[end]/values_A_S[1+shift])
        push!(CondSz, values_A_Sz[end]/values_A_Sz[1+shift])
        if verbose
            println("L = ", L)
            println("Condition number of A = ", CondA[end])
            println("Condition number of A_L = ", CondAL[end])
            println("Condition number of A_S = ", CondAS[end])
            println("Condition number of A_Sz = ", CondSz[end])
            println("=======================")
        end
    end
    p = plot(Scale.y_log10, Guide.manual_color_key("", ["A", "left BPX", "symmetric BPX", "Schwarz"], ["red", "black", "blue", "green"]), Guide.xlabel("level"), Guide.ylabel("condition number"))
    push!(p, layer(x=Ls, y=CondA, Geom.line, Geom.point, color=[colorant"red"]))
    push!(p, layer(x=Ls, y=CondAL, Geom.line, Geom.point, color=[colorant"black"]))
    push!(p, layer(x=Ls, y=CondAS, Geom.line, Geom.point, color=[colorant"blue"]))
    push!(p, layer(x=Ls, y=CondSz, Geom.line, Geom.point, color=[colorant"green"]))
    results = hcat(Ls, CondA, CondAS, CondAL, CondSz)
    writedlm("precs_comparison/convection_diffusion_"*string(v)*"_"*BCs[1]*BCs[2]*".txt", results)
    draw(SVG("precs_comparison/convection_diffusion_"*string(v)*"_"*BCs[1]*BCs[2]*".svg", 6inch, 4inch), p)
    return nothing
end

function test_Disc_Diffusion(BCs::AbstractArray{String, 1}, σ::Real, verbose=false)
    Ls = [5, 6, 7, 8, 9, 10, 11]
    CondA, CondAS, CondAL, CondSz = [], [], [], []
    for L = Ls
        F = [x -> -one(x)-σ*(x>0.5), x -> zero(x), x -> zero(x), x -> x^2]
        BCs_ = [[1, BCs[1]], [2, BCs[2]]]
        A, _, _ = BVP_1D(L, F, BCs_)
        B = BPX_matrix(L, BCs)
        Sz = Schwarz_matrix(A, L, BCs)
        A_L, A_S, A_Sz = Array(B*A), Array(B*A*B), Array(Sz*A)
        singular_values_A = svd(Array(A)).S
        values_A_L = sort(abs.(eigvals(A_L)))
        values_A_S = sort(abs.(eigvals(A_S)))
        values_A_Sz = sort(abs.(eigvals(A_Sz)))
        shift = Int64((BCs[1] == "N")*(BCs[2] == "N"))
        push!(CondA, singular_values_A[1]/singular_values_A[end-shift])
        push!(CondAL, values_A_L[end]/values_A_L[1+shift])
        push!(CondAS, values_A_S[end]/values_A_S[1+shift])
        push!(CondSz, values_A_Sz[end]/values_A_Sz[1+shift])
        if verbose
            println("L = ", L)
            println("Condition number of A = ", CondA[end])
            println("Condition number of A_L = ", CondAL[end])
            println("Condition number of A_S = ", CondAS[end])
            println("Condition number of A_Sz = ", CondSz[end])
            println("=======================")
        end
    end
    p = plot(Scale.y_log10, Guide.manual_color_key("", ["A", "left BPX", "symmetric BPX", "Schwarz"], ["red", "black", "blue", "green"]), Guide.xlabel("level"), Guide.ylabel("condition number"))
    push!(p, layer(x=Ls, y=CondA, Geom.line, Geom.point, color=[colorant"red"]))
    push!(p, layer(x=Ls, y=CondAL, Geom.line, Geom.point, color=[colorant"black"]))
    push!(p, layer(x=Ls, y=CondAS, Geom.line, Geom.point, color=[colorant"blue"]))
    push!(p, layer(x=Ls, y=CondSz, Geom.line, Geom.point, color=[colorant"green"]))
    results = hcat(Ls, CondA, CondAS, CondAL, CondSz)
    writedlm("precs_comparison/disc_diffusion"*string(σ)*"_"*BCs[1]*BCs[2]*".txt", results)
    draw(SVG("precs_comparison/disc_diffusion"*string(σ)*"_"*BCs[1]*BCs[2]*".svg", 6inch, 4inch), p)
    return nothing
end

function test_variable_Convection_Diffusion(BCs::AbstractArray{String, 1}, verbose=false)
    Ls = [5, 6, 7, 8, 9, 10, 11]
    CondA, CondAS, CondAL, CondSz = [], [], [], []
    for L = Ls
        F = [x -> -one(x)-cos(3*pi*x)^2, x -> sin(3*pi*x)^2, x -> zero(x), x -> x^2]
        BCs_ = [[1, BCs[1]], [2, BCs[2]]]
        A, _, _ = BVP_1D(L, F, BCs_)
        B = BPX_matrix(L, BCs)
        Sz = Schwarz_matrix(A, L, BCs)
        A_L, A_S, A_Sz = Array(B*A), Array(B*A*B), Array(Sz*A)
        singular_values_A = svd(Array(A)).S
        values_A_L = sort(abs.(eigvals(A_L)))
        values_A_S = sort(abs.(eigvals(A_S)))
        values_A_Sz = sort(abs.(eigvals(A_Sz)))
        shift = Int64((BCs[1] == "N")*(BCs[2] == "N"))
        push!(CondA, singular_values_A[1]/singular_values_A[end-shift])
        push!(CondAL, values_A_L[end]/values_A_L[1+shift])
        push!(CondAS, values_A_S[end]/values_A_S[1+shift])
        push!(CondSz, values_A_Sz[end]/values_A_Sz[1+shift])
        if verbose
            println("L = ", L)
            println("Condition number of A = ", CondA[end])
            println("Condition number of A_L = ", CondAL[end])
            println("Condition number of A_S = ", CondAS[end])
            println("Condition number of A_Sz = ", CondSz[end])
            println("=======================")
        end
    end
    p = plot(Scale.y_log10, Guide.manual_color_key("", ["A", "left BPX", "symmetric BPX", "Schwarz"], ["red", "black", "blue", "green"]), Guide.xlabel("level"), Guide.ylabel("condition number"))
    push!(p, layer(x=Ls, y=CondA, Geom.line, Geom.point, color=[colorant"red"]))
    push!(p, layer(x=Ls, y=CondAL, Geom.line, Geom.point, color=[colorant"black"]))
    push!(p, layer(x=Ls, y=CondAS, Geom.line, Geom.point, color=[colorant"blue"]))
    push!(p, layer(x=Ls, y=CondSz, Geom.line, Geom.point, color=[colorant"green"]))
    results = hcat(Ls, CondA, CondAS, CondAL, CondSz)
    writedlm("precs_comparison/variable_conv_diffusion_"*BCs[1]*BCs[2]*".txt", results)
    draw(SVG("precs_comparison/variable_conv_diffusion_"*BCs[1]*BCs[2]*".svg", 6inch, 4inch), p)
    return nothing
end

function test_implicit_Diffusion(BCs::AbstractArray{String, 1}, μ::Real, verbose=false)
    Ls = [5, 6, 7, 8, 9, 10, 11]
    CondA, CondAS, CondAL, CondSz = [], [], [], []
    for L = Ls
        F = [x -> μ*one(x), x -> zero(x), x -> zero(x), x -> x^2]
        BCs_ = [[1, BCs[1]], [2, BCs[2]]]
        A, _, _ = BVP_1D(L, F, BCs_)
        A = I - A
        B = BPX_matrix(L, BCs)
        Sz = Schwarz_matrix(A, L, BCs)
        A_L, A_S, A_Sz = Array(B*A), Array(B*A*B), Array(Sz*A)
        singular_values_A = svd(Array(A)).S
        values_A_L = sort(abs.(eigvals(A_L)))
        values_A_S = sort(abs.(eigvals(A_S)))
        values_A_Sz = sort(abs.(eigvals(A_Sz)))
        shift = Int64((BCs[1] == "N")*(BCs[2] == "N"))
        push!(CondA, singular_values_A[1]/singular_values_A[end-shift])
        push!(CondAL, values_A_L[end]/values_A_L[1+shift])
        push!(CondAS, values_A_S[end]/values_A_S[1+shift])
        push!(CondSz, values_A_Sz[end]/values_A_Sz[1+shift])
        if verbose
            println("L = ", L)
            println("Condition number of A = ", CondA[end])
            println("Condition number of A_L = ", CondAL[end])
            println("Condition number of A_S = ", CondAS[end])
            println("Condition number of A_Sz = ", CondSz[end])
            println("=======================")
        end
    end
    p = plot(Scale.y_log10, Guide.manual_color_key("", ["A", "left BPX", "symmetric BPX", "Schwarz"], ["red", "black", "blue", "green"]), Guide.xlabel("level"), Guide.ylabel("condition number"))
    push!(p, layer(x=Ls, y=CondA, Geom.line, Geom.point, color=[colorant"red"]))
    push!(p, layer(x=Ls, y=CondAL, Geom.line, Geom.point, color=[colorant"black"]))
    push!(p, layer(x=Ls, y=CondAS, Geom.line, Geom.point, color=[colorant"blue"]))
    push!(p, layer(x=Ls, y=CondSz, Geom.line, Geom.point, color=[colorant"green"]))
    results = hcat(Ls, CondA, CondAS, CondAL, CondSz)
    writedlm("precs_comparison/implicit_diffusion"*string(μ)*"_"*BCs[1]*BCs[2]*".txt", results)
    draw(SVG("precs_comparison/implicit_diffusion"*string(μ)*"_"*BCs[1]*BCs[2]*".svg", 6inch, 4inch), p)
    return nothing
end

function tests()
    BCS = [["D", "D"], ["D", "N"], ["N", "D"], ["N", "N"]]
    functions = [test_Laplace, x -> test_Helmholtz(x, 10), x -> test_Helmholtz(x, 100), x -> test_Convection_Diffusion(x, 100), x -> test_Convection_Diffusion(x, 1000), x -> test_Disc_Diffusion(x, 10), x -> test_Disc_Diffusion(x, 1000), test_variable_Convection_Diffusion, x -> test_implicit_Diffusion(x, 0.01), x -> test_implicit_Diffusion(x, 100)]
    for f=functions
        for BCs = BCS
            f(BCs)
        end
    end
end

max_stars_repo_name: VLSF/neuralBPX

max_issues_repo_name: VLSF/neuralBPX

max_forks_repo_name: VLSF/neuralBPX


text:
const libsigma = Pkg.dir("SDE", "deps", "libsigma")

#void fe_mu(double *mu, double *y, int N, int L, int K)
function fe_mu_c(y::Vector{Float64}, L, K)

	n = 2^(L-1)
	mu = zeros(Float64,2n-1 + K)
	N = length(y)
	product = ccall( (:fe_mu_at, libsigma),
                  Void,
                  (Ptr{Float64}, Ptr{Float64}, Int32, Int32),
                  mu, y, N, L)
	return mu;
end

#void fe_Sigma_at(double *S, double *y, int N, double dt, int L)
function fe_Sigma_c(y::Vector{Float64}, dt::Float64, L)
	n = 2^(L-1)
	N = length(y)
	S = zeros(Float64,(2n-1)^2)
	product = ccall( (:fe_Sigma_at,  libsigma),
                  Void,
                  (Ptr{Float64}, Ptr{Float64}, Int32, Float64, Int32),
                  S, y, N, dt, L)
	return reshape(S, 2n-1 , 2n-1)

end
function fe_SigmaB1_c(y::Vector{Float64}, dt::Float64, L)
	n = 2^(L-1)
	N = length(y)
	S = zeros(Float64,(2n)^2)
	product = ccall( (:fe_SigmaB1_at,  libsigma),
                  Void,
                  (Ptr{Float64}, Ptr{Float64}, Int32, Float64, Int32),
                  S, y, N, dt, L)
	return reshape(S, 2n , 2n)

end

fe_mu = fe_mu_c
fe_SigmaB1 = fe_SigmaB1_c
fe_Sigma = fe_Sigma_c


text:
"""Return the list of available libxc functionals as strings"""
function available_functionals()
    n_xc = ccall((:xc_number_of_functionals, libxc), Cint, ())
    max_string_length = ccall((:xc_maximum_name_length, libxc), Cint, ())

    funcnames = Vector{String}(undef, n_xc)
    for i in 1:n_xc
        funcnames[i] = "\0"^(max_string_length + 2)
    end
    ccall((:xc_available_functional_names, libxc), Cvoid, (Ptr{Ptr{UInt8}}, ), funcnames)
    [string(split(funcnames[i], "\0")[1]) for i in 1:n_xc]
end

@enum FunctionalKind begin
    functional_exchange             = 0
    functional_correlation          = 1
    functional_exchange_correlation = 2
    functional_kinetic              = 3
end

@enum FunctionalFamily begin
    family_unknown  = -1
    family_lda      = 1
    family_gga      = 2
    family_mggai    = 4
    family_lca      = 8
    family_oep      = 16
    family_hyb_gga  = 32
    family_hyb_mgga = 64
end

mutable struct Functional
    number::Int
    identifier::Symbol
    kind::FunctionalKind
    family::FunctionalFamily
    n_spin::Int

    # Pointer holding the LibXC representation of this functional
    pointer::Ptr{XCFuncType}
end

"""
    Functional(identifier::Symbol; n_spin::Integer = 1)

Construct a Functional from a libxc `identifier` and the number
of spins `n_spin` to consider. `
"""
function Functional(identifier::Symbol; n_spin::Integer = 1)
    if n_spin != 1 && n_spin != 2
        error("n_spin needs to be 1 or 2")
    end

    number = ccall((:xc_functional_get_number, libxc), Cint, (Cstring, ),
                   string(identifier))
    if number == -1
        error("Functional $identifier is not known.")
    end

    function pointer_cleanup(ptr::Ptr{XCFuncType})
        if ptr != C_NULL
            xc_func_end(ptr)
            xc_func_free(ptr)
        end
    end

    pointer = xc_func_alloc()
    ret = xc_func_init(pointer, number, n_spin)
    if ret != 0
        error("Something went wrong initialising the functional")
    end

    try
        funcinfo = xc_func_get_info(pointer)
        kind = xc_func_info_get_kind(funcinfo)
        family = xc_func_info_get_family(funcinfo)
        flags = xc_func_info_get_flags(funcinfo)
        # TODO Extract references ....

        # Make functional and attach finalizer for cleaning up the pointer
        func = Functional(number, identifier, FunctionalKind(kind),
                          FunctionalFamily(family), n_spin, pointer)
        finalizer(cls -> pointer_cleanup(cls.pointer), func)
        return func
    catch
        pointer_cleanup(pointer)
        rethrow()
    end
end


text:
@testset "Test" begin
    CP.Test.runtests(
        MOI.Utilities.MockOptimizer(
            MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Int}()),
        ),
        MOI.Test.Config(),
        warn_unsupported = true,
    )
end


text:
"""
#### Calculate the total mass of any ContainMassDataSetType:

```julia
msum(dataobject::ContainMassDataSetType; unit::Symbol=:standard, mask::MaskType=[false])

return Float64
```
#### Arguments
##### Required:
- **`dataobject`:** needs to be of type: "ContainMassDataSetType"

##### Optional Keywords:
- **`unit`:** the unit of the result (can be used w/o keyword): :standard (code units)  :Msol, :Mearth, :Mjupiter, :g, :kg  (of typye Symbol) ..etc. ; see for defined mass-scales viewfields(info.scale)
- **`mask`:** needs to be of type MaskType which is a supertype of Array{Bool,1} or BitArray{1} with the length of the database (rows)

"""
function msum(dataobject::ContainMassDataSetType, unit::Symbol; mask::MaskType=[false])
    return msum(dataobject, unit=unit, mask=mask)
end

function msum(dataobject::ContainMassDataSetType; unit::Symbol=:standard, mask::MaskType=[false])
    return sum( getvar(dataobject, :mass, unit=unit, mask=mask) )
end

"""
#### Calculate the center-of-mass of any ContainMassDataSetType:

```julia
center_of_mass(dataobject::ContainMassDataSetType; unit::Symbol=:standard, mask::MaskType=[false])

return Tuple{Float64, Float64, Float64,}
```
#### Arguments
##### Required:
- **`dataobject`:** needs to be of type: "ContainMassDataSetType"

##### Optional Keywords:
- **`unit`:** the unit of the result (can be used w/o keyword): :standard (code units), :Mpc, :kpc, :pc, :mpc, :ly, :au , :km, :cm (of typye Symbol) ..etc. ; see for defined length-scales viewfields(info.scale)
- **`mask`:** needs to be of type MaskType which is a supertype of Array{Bool,1} or BitArray{1} with the length of the database (rows)

"""
function center_of_mass(dataobject::ContainMassDataSetType, unit::Symbol; mask::MaskType=[false])
    return center_of_mass(dataobject, unit=unit, mask=mask)
end

function center_of_mass(dataobject::ContainMassDataSetType; unit::Symbol=:standard, mask::MaskType=[false])
    selected_unit = getunit(dataobject.info, unit)
    return ( average_mweighted(dataobject, :x, mask=mask), average_mweighted(dataobject, :y, mask=mask), average_mweighted(dataobject, :z,  mask=mask) ) .* selected_unit
end

"""
#### Calculate the center-of-mass of any ContainMassDataSetType:

```julia
com(dataobject::ContainMassDataSetType; unit::Symbol=:standard, mask::MaskType=[false])

return Tuple{Float64, Float64, Float64,}
```
#### Arguments
##### Required:
- **`dataobject`:** needs to be of type: "ContainMassDataSetType"

##### Optional Keywords:
- **`unit`:** the unit of the result (can be used w/o keyword): :standard (code units), :Mpc, :kpc, :pc, :mpc, :ly, :au , :km, :cm (of typye Symbol) ..etc. ; see for defined length-scales viewfields(info.scale)
- **`mask`:** needs to be of type MaskType which is a supertype of Array{Bool,1} or BitArray{1} with the length of the database (rows)

"""
function com(dataobject::ContainMassDataSetType, unit::Symbol; mask::MaskType=[false])
    return center_of_mass(dataobject, unit, mask=mask)
end

function com(dataobject::ContainMassDataSetType; unit::Symbol=:standard, mask::MaskType=[false])
    return center_of_mass(dataobject, unit=unit, mask=mask)
end

"""
#### Calculate the joint center-of-mass of any HydroPartType:

```julia
center_of_mass(dataobject::Array{HydroPartType,1}, unit::Symbol; mask::MaskArrayAbstractType=[[false],[false]])

return Tuple{Float64, Float64, Float64,}
```
#### Arguments
##### Required:
- **`dataobject`:** needs to be of type: "Array{HydroPartType,1}""

##### Optional Keywords:
- **`unit`:** the unit of the result (can be used w/o keyword): :standard (code units), :Mpc, :kpc, :pc, :mpc, :ly, :au , :km, :cm (of typye Symbol) ..etc. ; see for defined length-scales viewfields(info.scale)
- **`mask`:** needs to be of type MaskArrayAbstractType which contains two entries with supertype of Array{Bool,1} or BitArray{1} and the length of the database (rows)

"""
function center_of_mass(dataobject::Array{HydroPartType,1}, unit::Symbol; mask::MaskArrayAbstractType=[[false],[false]])
    return  center_of_mass(dataobject; unit=unit, mask=mask)
end

function center_of_mass(dataobject::Array{HydroPartType,1}; unit::Symbol=:standard, mask::MaskArrayAbstractType=[[false],[false]])

    selected_unit = getunit(dataobject[1].info, unit) # assuming both datasets are from same simulation output

    if length(mask[1]) == 1 && length(mask[2]) == 1
        m1 = getvar(dataobject[1], :mass)
        m1_sum = sum(m1)

        m2 = getvar(dataobject[2], :mass)
        m2_sum = sum(m2)

        m_sum = m1_sum + m2_sum

        x_weighted = (sum( getvar(dataobject[1], :x) .* m1 ) + sum( getvar(dataobject[2], :x) .* m2) ) / m_sum
        y_weighted = (sum( getvar(dataobject[1], :y) .* m1 ) + sum( getvar(dataobject[2], :y) .* m2) ) / m_sum
        z_weighted = (sum( getvar(dataobject[1], :z) .* m1 ) + sum( getvar(dataobject[2], :z) .* m2) ) / m_sum
    else
        m1 = getvar(dataobject[1], :mass)[mask[1]]
        m1_sum = sum(m1)

        m2 = getvar(dataobject[2], :mass)[mask[2]]
        m2_sum = sum(m2)

        m_sum = m1_sum + m2_sum

        x_weighted = (sum( getvar(dataobject[1], :x)[mask[1]] .* m1 ) + sum( getvar(dataobject[2], :x)[mask[2]] .* m2) ) / m_sum
        y_weighted = (sum( getvar(dataobject[1], :y)[mask[1]] .* m1 ) + sum( getvar(dataobject[2], :y)[mask[2]] .* m2) ) / m_sum
        z_weighted = (sum( getvar(dataobject[1], :z)[mask[1]] .* m1 ) + sum( getvar(dataobject[2], :z)[mask[2]] .* m2) ) / m_sum

    end
    return ( x_weighted, y_weighted, z_weighted ) .* selected_unit

end

"""
#### Calculate the joint center-of-mass of any HydroPartType:

```julia
com(dataobject::Array{HydroPartType,1}, unit::Symbol; mask::MaskArrayAbstractType=[[false],[false]])

return Tuple{Float64, Float64, Float64,}
```
#### Arguments
##### Required:
- **`dataobject`:** needs to be of type: "Array{HydroPartType,1}""

##### Optional Keywords:
- **`unit`:** the unit of the result (can be used w/o keyword): :standard (code units), :Mpc, :kpc, :pc, :mpc, :ly, :au , :km, :cm (of typye Symbol) ..etc. ; see for defined length-scales viewfields(info.scale)
- **`mask`:** needs to be of type MaskArrayAbstractType which contains two entries with supertype of Array{Bool,1} or BitArray{1} and the length of the database (rows)

"""
function com(dataobject::Array{HydroPartType,1}, unit::Symbol; mask::MaskArrayAbstractType=[[false],[false]])
    return  center_of_mass(dataobject, unit, mask=mask)
end

function com(dataobject::Array{HydroPartType,1}; unit::Symbol=:standard, mask::MaskArrayAbstractType=[[false],[false]])
    return center_of_mass(dataobject, unit=unit, mask=mask)
end

function average_mweighted(dataobject::ContainMassDataSetType, var::Symbol; mask::MaskType=[false])
    return sum( getvar(dataobject, var, mask=mask) .* getvar(dataobject, :mass, mask=mask) ) ./ sum( getvar(dataobject, :mass, mask=mask))
end

"""
#### Calculate the average velocity (w/o mass-weight) of any ContainMassDataSetType:

```julia
bulk_velocity(dataobject::ContainMassDataSetType; unit::Symbol=:standard, weighting::Symbol=:mass, mask::MaskType=[false])

return Tuple{Float64, Float64, Float64,}
```
#### Arguments
##### Required:
- **`dataobject`:** needs to be of type: "ContainMassDataSetType"

##### Optional Keywords:
- **`unit`:** the unit of the result (can be used w/o keyword): :standard (code units)  :km_s, :m_s, :cm_s (of typye Symbol) ..etc. ; see for defined velocity-scales viewfields(info.scale)
- **`weighting`:** use different weightings: :mass (default), :volume (hydro), :no
- **`mask`:** needs to be of type MaskType which is a supertype of Array{Bool,1} or BitArray{1} with the length of the database (rows)

"""
function bulk_velocity(dataobject::ContainMassDataSetType, unit::Symbol; weighting::Symbol=:mass, mask::MaskType=[false])
    return bulk_velocity(dataobject, unit=unit, weighting=weighting, mask=mask)
end

function bulk_velocity(dataobject::ContainMassDataSetType; unit::Symbol=:standard, weighting::Symbol=:mass, mask::MaskType=[false])
    selected_unit = getunit(dataobject.info, unit)
    if weighting == :mass
        return ( average_mweighted(dataobject, :vx, mask=mask), average_mweighted(dataobject, :vy, mask=mask), average_mweighted(dataobject, :vz,  mask=mask) ) .* selected_unit
    elseif weighting == :volume && typeof(dataobject) == HydroDataType
        isamr = checkuniformgrid(dataobject, dataobject.lmax)
        if isamr
            return ( sum( getvar(dataobject, :vx, mask=mask) .* getvar(dataobject, :volume, mask=mask) ) ./ sum( getvar(dataobject, :volume, mask=mask) ),
                     sum( getvar(dataobject, :vy, mask=mask) .* getvar(dataobject, :volume, mask=mask) ) ./ sum( getvar(dataobject, :volume, mask=mask) ),
                     sum( getvar(dataobject, :vz, mask=mask) .* getvar(dataobject, :volume, mask=mask) ) ./ sum( getvar(dataobject, :volume, mask=mask) ) ) .* selected_unit
        else
            return ( mean( getvar(dataobject, :vx, mask=mask) ), mean( getvar(dataobject, :vy, mask=mask) ), mean( getvar(dataobject, :vz,  mask=mask)) ) .* selected_unit
        end
    elseif weighting == :no # for AMR
            return ( mean( getvar(dataobject, :vx, mask=mask) ), mean( getvar(dataobject, :vy, mask=mask) ), mean( getvar(dataobject, :vz,  mask=mask)) ) .* selected_unit
    end
end

"""
#### Calculate the average velocity (w/o mass-weight) of any ContainMassDataSetType:

```julia
average_velocity(dataobject::ContainMassDataSetType; unit::Symbol=:standard, weighting::Symbol=:mass, mask::MaskType=[false])

return Tuple{Float64, Float64, Float64,}
```
#### Arguments
##### Required:
- **`dataobject`:** needs to be of type: "ContainMassDataSetType"

##### Optional Keywords:
- **`unit`:** the unit of the result (can be used w/o keyword): :standard (code units)  :km_s, :m_s, :cm_s (of typye Symbol) ..etc. ; see for defined velocity-scales viewfields(info.scale)
- **`weighting`:** use different weightings: :mass (default), :volume (hydro), :no
- **`mask`:** needs to be of type MaskType which is a supertype of Array{Bool,1} or BitArray{1} with the length of the database (rows)

"""
function average_velocity(dataobject::ContainMassDataSetType, unit::Symbol; weighting::Symbol=:mass, mask::MaskType=[false])
    return bulk_velocity(dataobject, unit, weighting=weighting, mask=mask)
end

function average_velocity(dataobject::ContainMassDataSetType; unit::Symbol=:standard, weighting::Symbol=:mass, mask::MaskType=[false])
    return bulk_velocity(dataobject, unit=unit, weighting=weighting,  mask=mask)
end

"""
#### Calculate statistical values w/o weighting of any Array:

```julia
wstat(array::Array{<:Real,1}; weight::Array{<:Real,1}=[1.], mask::MaskType=[false])

WStatType(mean, median, std, skewness, kurtosis, min, max)
```
#### Arguments
##### Required:
- **`array`:** Array needs to be of type: "<:Real"

##### Optional Keywords:
- **`weight`:** Array needs to be of type: "<:Real" (can be used w/o keyword)
- **`mask`:** needs to be of type MaskType which is a supertype of Array{Bool,1} or BitArray{1} with the length of the Array

"""
function wstat(array::Array{<:Real,1}, weight::Array{<:Real,1}; mask::MaskType=[false])
    return  wstat(array, weight=weight, mask=mask)
end

function wstat(array::Array{<:Real,1}; weight::Array{<:Real,1}=[1.], mask::MaskType=[false])

            if length(mask) > 1
                array=array[mask]
                if length(weight) > 1
                    weight=weight[mask]
                end
            end

        if length(weight)==1
            mean_      = mean(array)
            median_    = median(array)
            std_       = std(array, mean=mean_)
            skewness_  = skewness(array, mean_)
            kurtosis_  = kurtosis(array, mean_)
            min_       = minimum(array)
            max_       = maximum(array)
            return WStatType(mean_, median_, std_, skewness_, kurtosis_, min_, max_)

        else length(weight) > 1
            mean_      = mean( array, weights( weight ))
            median_    = median(array, weights( weight ))
            std_       = std(array, mean=mean_, weights( weight ), corrected=false)
            skewness_  = skewness(array, mean_)
            kurtosis_  = kurtosis(array, mean_)
            min_       = minimum(array)
            max_       = maximum(array)
            return WStatType(mean_, median_, std_, skewness_, kurtosis_, min_, max_)
        end

end


text:
using JLLisp
using Test
eval_(x) = JLLisp.Eval.eval_(x)
read_ = JLLisp.Reader.read
symT = JLLisp.Symbols.symbolT
Nil = JLLisp.Nil

@testset "JLLisp.jl" begin
    @test Atom <: T
    @test Number_ <: Atom
    @test List <: T
    @test Null <: List
    @test Nil == Null()
    @test string(T) == "T"
    @test string(Nil) == "NIL"
    x = JLLisp.Integer__.Integer_(1)
    y = JLLisp.Integer__.Integer_(2)
    @test JLLisp.Integer__.add(x,y) == JLLisp.Integer__.Integer_(3)
    @test string(x) == "1"

    @test isa(x, JLLisp.Atom)
    @test x == JLLisp.Eval.eval_(x) # Atom
    @test isa(Nil, JLLisp.Null)
    @test JLLisp.Eval.eval_(Nil) == Nil # Null

    z = JLLisp.Cons_.Cons(x, JLLisp.Cons_.Cons(y, Nil))
    @test JLLisp.Integer__.Integer_(3) == JLLisp.Function_.funcall(JLLisp.Function_.Add(), z)

    form = JLLisp.Cons_.Cons()
    form.car = JLLisp.Symbols.symboltable["+"]
    form.cdr = JLLisp.Cons_.Cons(x, JLLisp.Cons_.Cons(y, Nil))
    @test eval_(form) == JLLisp.Integer__.Integer_(3)
    @test eval_(read_(ex = "(+ 1 2)")) == JLLisp.Integer__.Integer_(3)
    @test eval_(read_(ex = "(+ -1 2)")) == JLLisp.Integer__.Integer_(1)
    @test eval_(read_(ex = "(- 1 2)")) == JLLisp.Integer__.Integer_(-1)
    @test eval_(read_(ex = "(- 1 -2)")) == JLLisp.Integer__.Integer_(3)
    @test eval_(read_(ex = "(* 2 3)")) == JLLisp.Integer__.Integer_(6)
    @test eval_(read_(ex = "(* 2 -3)")) == JLLisp.Integer__.Integer_(-6)
    @test eval_(read_(ex = "(/ 1 2)")) == JLLisp.Integer__.Integer_(0)
    @test eval_(read_(ex = "(/ 3 2)")) == JLLisp.Integer__.Integer_(1)
    @test eval_(read_(ex = "(/ -3 2)")) == JLLisp.Integer__.Integer_(-1)
    @test eval_(read_(ex = "(+ (+ 1 2) (+ 3 4))")) == JLLisp.Integer__.Integer_(10)
    @test eval_(read_(ex = "(- (+ -1 2) (- 1 2))")) == JLLisp.Integer__.Integer_(2)
    @test eval_(read_(ex = "(* (+ 1 2) (- 2 4))")) == JLLisp.Integer__.Integer_(-6)
    @test eval_(read_(ex = "(/ (+ 100 -2) (- 4 2))")) == JLLisp.Integer__.Integer_(49)
    @test eval_(read_(ex = "(<= 1 1)")) == symT
    @test eval_(read_(ex = "(<= -3 2)")) == symT
    @test eval_(read_(ex = "(>= 1 1)")) == symT
    @test eval_(read_(ex = "(>= -3 2)")) == JLLisp.Nil
    @test eval_(read_(ex = "(< 1 1)")) == JLLisp.Nil
    @test eval_(read_(ex = "(< -3 2)")) == symT
    @test eval_(read_(ex = "(> 1 1)")) == JLLisp.Nil
    @test eval_(read_(ex = "(> 2 -3)")) == symT
    @test eval_(read_(ex = "(= 1 1)")) == symT
    @test eval_(read_(ex = "(= -3 2)")) == JLLisp.Nil
    @test eval_(read_(ex="(setq hoge 10)")) == JLLisp.Integer__.Integer_(10)
    @test eval_(read_(ex="hoge")) == JLLisp.Integer__.Integer_(10)
    @test eval_(read_(ex="(if (< 10 20) 100 200)")) == JLLisp.Integer__.Integer_(100)
    @test eval_(read_(ex="(if (> 10 20) 100 200)")) == JLLisp.Integer__.Integer_(200)
    @test eval_(read_(ex="(type-of 1)")) == JLLisp.Symbols.symbol_("JLLISP.INTEGER__.INTEGER_")
    eval_(read_(ex = "(defun f () (+ 10 20))"))
    @test eval_(read_(ex = "(f)")) == JLLisp.Integer__.Integer_(30)

    eval_(read_(ex = "(defun g (x) (+ 10 x))"))
    @test eval_(read_(ex = "(g 100)")) == JLLisp.Integer__.Integer_(110)

    eval_(read_(ex = "(defun fac (n) (if (< n 1) 1 (* n (fac (- n 1)))))"))
    @test eval_(read_(ex = "(fac 5)")) == JLLisp.Integer__.Integer_(120)
end


text:
"""
    AbstractSymmetry

An abstract types which dictates the symmetry of the setup. That is, the symmetry shared between the incident wave and the shape of the material.
"""
abstract type AbstractSymmetry{Dim} end

"""
    WithoutSymmetry

A type used to describe materials and incident waves which share no common symmetry. This will lead to the most general regular wave expansion for the eignvectors.
"""
struct WithoutSymmetry{Dim} <: AbstractSymmetry{Dim} end

"""
An incident plane-wave and halfspace material will result in all fields being plane-waves.
"""
abstract type AbstractPlanarSymmetry{Dim} <: AbstractSymmetry{Dim} end

"""
    PlanarSymmetry

A type used to describe materials and incident waves that both share a planar symmetry.
"""
struct PlanarSymmetry{Dim} <: AbstractPlanarSymmetry{Dim} end

"""
For spatial dimension > 2, we can consider problems that have azimuthal symmetry. For example, a plane-wave incident on a sphere.
"""
abstract type AbstractAzimuthalSymmetry{Dim} <: AbstractSymmetry{Dim} end

"""
    AzimuthalSymmetry

A type used to describe materials and incident waves in 3 spatial dimensions that share a symmetry about the azimuth.
"""
struct AzimuthalSymmetry{Dim} <: AbstractAzimuthalSymmetry{Dim} end
AzimuthalSymmetry() = AzimuthalSymmetry{3}()

"""
    RadialSymmetry

A type used to describe materials and incident waves in that are both radially symmetric. That is, the material is a sphere, and the incident wave is radially symmetric around the center of this spherical material.
"""
struct RadialSymmetry{Dim} <: AbstractAzimuthalSymmetry{Dim} end

"""
    PlanarAzimuthalSymmetry

A type used to describe a materail and incident wave which have both [`PlanarSymmetry`](@ref) and [`AzimuthalSymmetry`](@ref).
"""
struct PlanarAzimuthalSymmetry{Dim} <: AbstractPlanarSymmetry{Dim} end
PlanarAzimuthalSymmetry() = PlanarAzimuthalSymmetry{3}()

function azimuthalnormal(dim::Int)
    if dim == 2
        return [1,0]
    elseif dim == 3
        return [0,0,1]
    else
        return [zeros(dim-1);1]
    end
end

# If T1 is a shape, source, or material, first calculate its symmetry
"""
    Symmetry(T)

returns the symmetry of the object T. Some possible symmetries include [`PlanarSymmetry`](@ref),  [`RadialSymmetry`](@ref),  and [`WithoutSymmetry`](@ref).
"""

"""
    Symmetry(T1,T2)

returns the combined symmetry of the two objects T1 and T2.
"""
Symmetry(T1,T2) = Symmetry(Symmetry(T1),Symmetry(T2))

Symmetry(sym1::AbstractSymmetry{Dim},sym2::AbstractSymmetry{Dim}) where Dim = WithoutSymmetry{Dim}()
Symmetry(sym1::S,sym2::S) where S<:AbstractSymmetry = sym1

Symmetry(sym1::AbstractPlanarSymmetry{Dim},sym2::AbstractPlanarSymmetry{Dim}) where Dim = PlanarSymmetry{Dim}()
Symmetry(sym1::PlanarAzimuthalSymmetry{Dim},sym2::PlanarAzimuthalSymmetry{Dim}) where Dim = PlanarAzimuthalSymmetry{Dim}()

Symmetry(sym1::AzimuthalSymmetry{Dim},sym2::PlanarAzimuthalSymmetry{Dim}) where Dim = AzimuthalSymmetry{Dim}()
Symmetry(sym1::PlanarAzimuthalSymmetry{Dim},sym2::AzimuthalSymmetry{Dim}) where Dim = AzimuthalSymmetry{Dim}()

Symmetry(sym1::AbstractAzimuthalSymmetry{Dim},sym2::RadialSymmetry{Dim}) where Dim = AzimuthalSymmetry{Dim}()
Symmetry(sym1::RadialSymmetry{Dim},sym2::AbstractAzimuthalSymmetry{Dim}) where Dim = AzimuthalSymmetry{Dim}()

Symmetry(sym1::RadialSymmetry{Dim},sym2::RadialSymmetry{Dim}) where Dim = RadialSymmetry{Dim}()

Symmetry(sym1::PlanarAzimuthalSymmetry{Dim},sym2::RadialSymmetry{Dim}) where Dim = AzimuthalSymmetry{Dim}()
Symmetry(sym1::RadialSymmetry{Dim},sym2::PlanarAzimuthalSymmetry{Dim}) where Dim = AzimuthalSymmetry{Dim}()


text:
# test_convergence.jl

# example problem from On the Robustness of a Multiple Semicoarsened Grid Method, Oosterlee and Wesseling.

function get_problem(ϵ,β)
    n = 256
    m = n
    A = anisotropic2d(ϵ, β*π/180, n, m)
    A = V_cycle(A, (n,m))
    b = ones((n-1)*(m-1))
    return (A,b)
end

@testset "Convergence for example from Oosterlee et. al." begin
    for ϵ in [1 1e-2 1e-8]
        for β in [0 15 30 45]
            (A,b) = get_problem(ϵ,β)
            A.grids[1].b .= b # copy rhs
            for item in Base.Iterators.take(A,15) end # iterate
            if !( ( ϵ == 1e-8 ) && ( β == 45 ) ) # known failure for last case in 15 iters
                @test A.resnorm[end] < 1/256^2
            end
            log(A,ϵ,β)
        end
    end
end


text:
#Helper functions for ratings RatingSystems - duplicating, getting games already played etc

function dupe_for_rating(games)
    #Takes dataframe of games and predictions, reverses the players, result and prediction and appends them to end#

    dupe_backwards = rename(games, :P1 => :P2, :P2 => :P1, :P1_wins => :P2_wins, :P2_wins => :P1_wins)
    duplicated_games = vcat(games, dupe_backwards)
    sort!(duplicated_games, :Period)

    return duplicated_games
end

function brier(predictions)
    #Brier score simple (versus all-wins)
    withoutmissing = filter(!isnan, collect(skipmissing(predictions)))
    return sum((1.0 .- withoutmissing).^2) / length(collect(skipmissing(predictions)))
end

function brier(predictions, true_scores)
    #Brier score with actual scores comparison
    return sum(collect(skipmissing(true_scores .- predictions).^2)) / length(collect(skipmissing(predictions)))
end

function brierss(predictions)
    #Brier skill score -- compares to unskilled = (1-0.5)^2 forecast.
    return 1 - (brier(predictions) / 0.25)
end

function cumcount(m)
    #Cumulative count occurrences
    d = Dict()
    function g(a)
        d[a] = get(d, a, 0) + 1
    end
    return g.(m)
end

function cumcount(m, surface)
    d = Dict()
    for i in 1:maximum(surface)
        d[i] = Dict()
    end

    function g(a, surface)
        d[surface][a] = get(d[surface], a, 0) + 1
    end

    return g.(m, surface)
end

function actual_scores(games)
    function score_safe(x, y)
        if x + y == 0
            return 0.5
        else
            return x / (x+y)
        end
    end
    return score_safe.(games.P1_wins, games.P2_wins)
end

function games_previously_played_count(training_games)
    m = vcat(training_games.P1, training_games.P2)
    flattened = collect(Iterators.flatten(transpose(m)))
    counts = cumcount(flattened)
    return transpose(reshape(counts, 2, :))
end

function games_previously_played_count(training_games, testing_games)
    m = hcat(vcat(training_games.P1, testing_games.P2), vcat(training_games.P2, testing_games.P1))
    flattened = collect(Iterators.flatten(transpose(m)))
    counts = cumcount(flattened)
    return transpose(reshape(counts, 2, :))
end

function filter_today_played_before(training_games, today_names, threshold)
    #Filter predictions for day on threshold of games played in training set
    relevant = games_previously_played_count(training_games, today_names)[size(training_games)[1] + 1: size(training_games)[1] + size(today_names)[1],:]

    today_names.P1_before = relevant[:, 1]
    today_names.P2_before = relevant[:, 2]

    return filter(row -> row.P1_before >= threshold && row.P2_before >= threshold, today_names)
end

function mask_played_before(training_games, testing_games, threshold)
    #Filter testing_games for those played before
    played_before = minimum(games_previously_played_count(training_games, testing_games), dims = 2)
    return played_before[end - size(testing_games, 1) + 1 : end] .>= threshold
end

function recent_players(games::DataFrame, recent_days::Int)
    #Return players who have had a game in the last few days
    latest_day = maximum(games.Day)
    games_played_recently = filter(row -> row.Day >= latest_day - recent_days, games)
    return unique(vcat(games_played_recently.P1, games_played_recently.P2))
end

function log_loss(predicted, actual; tol = 0.001)
    predicted_clean = predicted[clean_nan_missing.(predicted)]
    actual_clean = actual[clean_nan_missing.(predicted)]
    predicted_capped = min.(max.(predicted_clean, tol), 1 .- tol)
	return -sum((actual_clean .* log.(predicted_capped)) .+ ((1 .- actual_clean) .* log.(1 .- predicted_capped))) / length(actual_clean)
end

function clean_nan_missing(m)
    return !(ismissing(m) | isnan(m))
end

function log_loss(predicted; tol = 0.001)
    predicted_clean = predicted[clean_nan_missing.(predicted)]
	log_loss(predicted_clean, ones(length(predicted_clean)); tol = tol)
end

function add_blanks(playerdayratings::Dict)
    days = collect(keys(sort(playerdayratings)))
    first_day = minimum(days)
    last_day = maximum(days)
    all_days = fill(NaN, last_day - first_day + 1)
    for j in days
        all_days[j - first_day + 1] = playerdayratings[j]
    end
    return all_days
end

function order_by_recency(games::DataFrame)
    P1 = vcat(games.P1)
    P2 = vcat(games.P2)
    return reverse(unique(transpose(hcat(P1, P2))))
end

function change_in_surface(training_games::DataFrame, testing_games::DataFrame, surface_training, surface_testing)
    m = hcat(vcat(training_games.P1, testing_games.P2), vcat(training_games.P2, testing_games.P1))
    flattened = collect(Iterators.flatten(transpose(m)))

    s = hcat(vcat(surface_training, surface_testing), vcat(surface_training, surface_testing))
    flattened_surface = collect(Iterators.flatten(transpose(s)))

    last_surface_played = transpose(reshape(last_surface(flattened, flattened_surface), 2, :))
    return .!isequal.(s, last_surface_played)
end

function last_surface(m, s)
    #Return the last surface and update
    d = Dict()
    function g(a, s)
        b = get(d, a, s)
        d[a] = s
        return b
    end
    return g.(m, s)
end

function jumble(matrix)
    #Run through and reverse half the games for ratings assessment
    row_length = size(matrix)[1]
    true_scores = ones(Float64, row_length)
    jumbled_matrix = similar(matrix)

    for row in 1:row_length

        if isodd(row)
            jumbled_matrix[row, :] = 1 .- matrix[row, :]
            true_scores[row] = 0.0
        else
            jumbled_matrix[row, :] = matrix[row, :]
            true_scores[row] = 1.0
        end
    end
    return jumbled_matrix, true_scores
end

function unscrambletuples(arr)
    r1 = zeros(Float64, size(arr)[1])
    r2 = zeros(Float64, size(arr)[1])
    var1 = zeros(Float64, size(arr)[1])
    var2 = zeros(Float64, size(arr)[1])
    for i in 1:size(arr)[1]
        r1[i] = arr[i][1][1]
        r2[i] = arr[i][2][1]
        var1[i] = arr[i][1][2]
        var2[i] = arr[i][2][2]
    end
    return r1, r2, var1, var2
end


text:
using TableReader

# Measure elapsed time.
filename = ARGS[1]
reader = occursin(".csv", filename) ? readcsv :
         occursin(".tsv", filename) ? readtsv :
         error("the extension of the filename must .csv or .tsv")
if length(ARGS) ≥ 2
    kwargs = eval(Meta.parse(ARGS[2]))
else
    kwargs = ()
end
for _ in 1:3
    @time reader(filename; kwargs...)
end

# Profile.
using Profile
@profile reader(filename; kwargs...)
Profile.print(IOContext(stdout, :displaysize => (10000, 10000)))


text:
module Simulation

export PlanktonSimulation
export update!, vel_copy!

using CUDA
using StructArrays
using LinearAlgebra

using PlanktonIndividuals.Architectures
using PlanktonIndividuals.Grids
using PlanktonIndividuals.Fields
using PlanktonIndividuals.Biogeochemistry
using PlanktonIndividuals.Parameters
using PlanktonIndividuals.Plankton
using PlanktonIndividuals.Diagnostics
using PlanktonIndividuals.Model
using PlanktonIndividuals.Output

using PlanktonIndividuals: AbstractMode, CarbonMode, QuotaMode, MacroMolecularMode

import Base: show

include("utils.jl")
include("simulations.jl")
include("update.jl")

end