text:
module TestData
    importall Base # so that we get warnings for conflicts
    using Base.Test
    using DataFrames

    #test_group("DataVector creation")
    dvint = @data([1, 2, NA, 4])
    dvint2 = data([5:8])
    dvint3 = data(5:8)
    dvflt = @data([1.0, 2, NA, 4])
    dvstr = @data(["one", "two", NA, "four"])
    dvdict = DataArray(Dict, 4)    # for issue #199

    #test_group("constructors")
    df1 = DataFrame({dvint, dvstr}, [:Ints, :Strs])
    df2 = DataFrame({dvint, dvstr})
    df3 = DataFrame({dvint})
    df4 = convert(DataFrame, [1:4 1:4])
    df5 = DataFrame({@data([1,2,3,4]), dvstr})
    df6 = DataFrame({dvint, dvint, dvstr}, [:A, :B, :C])
    df7 = DataFrame(x = dvint, y = dvstr)
    @test size(df7) == (4, 2)
    @test isequal(df7[:x], dvint)

    #test_group("description functions")
    @test size(df6, 1) == 4
    @test size(df6, 2) == 3
    @test all(names(df6) .== [:A, :B, :C])
    @test all(names(df2) .== [:x1, :x2])
    @test all(names(df7) .== [:x, :y])

    #test_group("ref")
    @test df6[2, 3] == "two"
    @test isna(df6[3, 3])
    @test df6[2, :C] == "two"
    @test isequal(df6[:B], dvint)
    @test size(df6[[2,3]], 2) == 2
    @test size(df6[2,:], 1) == 1
    @test size(df6[[1, 3], [1, 3]]) == (2, 2)
    @test size(df6[1:2, 1:2]) == (2, 2)
    @test size(head(df6,2)) == (2, 3)
    # lots more to do

    #test_group("combining")

    dfc = hcat(df3, df4)
    @test size(dfc, 2) == 3
    @test all(names(dfc) .== [:x1, :x1_1, :x2])
    @test isequal(dfc[:x1], df3[:x1])

    @test isequal(dfc, [df3 df4])

    dfr = vcat(df4, df4)
    @test size(dfr, 1) == 8
    @test all(names(df4) .== names(dfr))
    @test isequal(dfr, [df4, df4])

    dfr = vcat(df2, df3)
    @test size(dfr) == (8,2)
    @test all(names(df2) .== names(dfr))
    @test isna(dfr[8,:x2])

    #test_group("assign")
    df6[3] = @data(["un", "deux", "troix", "quatre"])
    @test df6[1, 3] == "un"
    df6[:B] = [4, 3, 2, 1]
    @test df6[1,2] == 4
    df6[:D] = [true, false, true, false]
    @test df6[1,4] == true
    delete!(df6, :D)
    @test all(names(df6) .== [:A, :B, :C])
    @test size(df6, 2) == 3

    #test_group("NA handling")
    @test nrow(df5[complete_cases(df5), :]) == 3

    #test_context("SubDataFrames")

    #test_group("constructors")
    # single index is rows
    sdf6a = sub(df6, 1)
    sdf6b = sub(df6, 2:3)
    sdf6c = sub(df6, [true, false, true, false])
    @test size(sdf6a) == (1,3)
    sdf6d = sub(df6, [1,3], :B)
    @test size(sdf6d) == (2,1)

    #test_group("ref")
    @test sdf6a[1,2] == 4

    #test_context("Within")
    #test_group("Associative")

    #test_group("DataFrame")
    srand(1)
    N = 20
    #Cast to int64 as rand() behavior differs between Int32/64
    d1 = pdata(rand(int64(1:2), N))
    d2 = (@pdata ["A", "B", NA])[rand(int64(1:3), N)]
    d3 = data(randn(N))
    d4 = data(randn(N))
    df7 = DataFrame({d1, d2, d3}, [:d1, :d2, :d3])

    #test_group("groupby")
    gd = groupby(df7, :d1)
    @test length(gd) == 2
    # @test isequal(gd[2]["d2"], PooledDataVector["A", "B", NA, "A", NA, NA, NA, NA])
    @test sum(gd[2][:d3]) == sum(df7[:d3][dropna(df7[:d1] .== 2)])

    g1 = groupby(df7, [:d1, :d2])
    g2 = groupby(df7, [:d2, :d1])
    @test sum(g1[1][:d3]) == sum(g2[1][:d3])

    res = 0.0
    for x in g1
        res += sum(x[:d1])
    end
    @test res == sum(df7[:d1])

    # TODO: Don't use symbols here
    df8 = colwise(df7[[1, 3]], :sum)
    @test df8[1, :d1_sum] == sum(df7[:d1])

    df8 = colwise(groupby(df7, :d2), [:sum, :length])
    @test size(df8, 1) == 3
    @test size(df8, 2) == 5
    @test df8[2, :d1_length] == 8

    df9 = df7 |> groupby([:d2]) |> [:sum, :length]
    @test isequal(df9, df8)
    df9 = by(df7, :d2, [:sum, :length])
    @test isequal(df9, df8)

    #test_group("reshape")
    d1 = DataFrame(a = repeat([1:3], inner = [4]),
                   b = repeat([1:4], inner = [3]),
                   c = randn(12),
                   d = randn(12))

    d1s = stack(d1, [:a, :b])
    d1s2 = stack(d1, [:c, :d])
    d1s3 = melt(d1, [:c, :d])
    @test isequal(d1s[1:12, :c], d1[:c])
    @test isequal(d1s[13:24, :c], d1[:c])
    @test all(names(d1s) .== [:variable, :value, :c, :d])
    @test isequal(d1s, d1s3)
    d1s_df = stackdf(d1, [:a, :b])
    @test isequal(d1s[:variable], d1s_df[:variable][:])
    @test isequal(d1s[:value], d1s_df[:value][:])
    @test isequal(d1s[:c], d1s_df[:c][:])
    @test isequal(d1s[1,:], d1s_df[1,:])

    d1s[:idx] = [1:12, 1:12]
    d1s2[:idx] = [1:12, 1:12]
    d1us = unstack(d1s, :variable, :idx, :value)
    d1us2 = unstack(d1s2, :variable, :idx, :value)
    @test isequal(d1us[:a], d1[:a])
    @test isequal(d1us2[:d], d1[:d])

    const letters = convert(Vector{ASCIIString}, split("abcdefghijklmnopqrstuvwxyz", ""))
    const LETTERS = convert(Vector{ASCIIString}, split("ABCDEFGHIJKLMNOPQRSTUVWXYZ", ""))

    #test_group("merge")

    srand(1)
    df1 = DataFrame(a = shuffle!([1:10]),
                    b = [:A,:B][rand(1:2, 10)],
                    v1 = randn(10))

    df2 = DataFrame(a = shuffle!(reverse([1:5])),
                    b2 = [:A,:B,:C][rand(1:3, 5)],
                    v2 = randn(5))

    m1 = join(df1, df2, on = :a)
    @test isequal(m1[:a], @data([1, 2, 3, 4, 5]))
    # TODO: Re-enable
    # m2 = join(df1, df2, on = :a, kind = :outer)
    # @test isequal(m2[:b2], DataVector["A", "B", "B", "B", "B", NA, NA, NA, NA, NA])
    # @test isequal(m2[:b2], DataVector["B", "B", "B", "C", "B", NA, NA, NA, NA, NA])

    df1 = DataFrame(a = [1, 2, 3],
                    b = ["America", "Europe", "Africa"])
    df2 = DataFrame(a = [1, 2, 4],
                    c = ["New World", "Old World", "New World"])

    m1 = join(df1, df2, on = :a, kind = :inner)
    @test isequal(m1[:a], @data([1, 2]))

    m2 = join(df1, df2, on = :a, kind = :left)
    @test isequal(m2[:a], @data([1, 2, 3]))

    m3 = join(df1, df2, on = :a, kind = :right)
    @test isequal(m3[:a], @data([1, 2, 4]))

    m4 = join(df1, df2, on = :a, kind = :outer)
    @test isequal(m4[:a], @data([1, 2, 3, 4]))

    # test with NAs (issue #185)
    df1 = DataFrame()
    df1[:A] = @data(["a", "b", "a", NA])
    df1[:B] = @data([1, 2, 1, 3])

    df2 = DataFrame()
    df2[:A] = @data(["a", NA, "c"])
    df2[:C] = @data([1, 2, 4])

    m1 = join(df1, df2, on = :A)
    @test size(m1) == (3,3)
    @test isequal(m1[:A], @data([NA,"a","a"]))

    m2 = join(df1, df2, on = :A, kind = :outer)
    @test size(m2) == (5,3)
    @test isequal(m2[:A], @data([NA,"a","a","b","c"]))

    srand(1)
    df1 = DataFrame(
        a = [:x,:y][rand(1:2, 10)],
        b = [:A,:B][rand(1:2, 10)],
        v1 = randn(10)
    )

    df2 = DataFrame(
        a = [:x,:y][[1,2,1,1,2]],
        b = [:A,:B,:C][[1,1,1,2,3]],
        v2 = randn(5)
    )
    df2[1,:a] = NA

    # # TODO: Restore this functionality
    # m1 = join(df1, df2, on = [:a,:b])
    # @test isequal(m1[:a], DataArray(["x", "x", "y", "y", fill("x", 5)]))
    # m2 = join(df1, df2, on = ["a","b"], kind = :outer)
    # @test isequal(m2[10,:v2], NA)
    # @test isequal(m2[:a], DataVector["x", "x", "y", "y", "x", "x", "x", "x", "x", "y", NA, "y"])

    srand(1)
    function spltdf(d)
        d[:x1] = map(x -> x[1], d[:a])
        d[:x2] = map(x -> x[2], d[:a])
        d[:x3] = map(x -> x[3], d[:a])
        d
    end
    df1 = DataFrame(
        a = ["abc", "abx", "axz", "def", "dfr"],
        v1 = randn(5)
    )
    df1 = spltdf(df1)
    df2 = DataFrame(
        a = ["def", "abc","abx", "axz", "xyz"],
        v2 = randn(5)
    )
    df2 = spltdf(df2)

    # m1 = join(df1, df2, on = :a)
    # m2 = join(df1, df2, on = [:x1, :x2, :x3])
    # @test isequal(sort(m1[:a]), sort(m2[:a]))

    #test_group("New DataVector constructors")
    dv = DataArray(Int, 5)
    @test all(isna(dv))
    dv = DataArray(Float64, 5)
    @test all(isna(dv))
    dv = @data(zeros(5))
    @test all(dv .== 0.0)
    dv = @data(ones(5))
    @test all(dv .== 1.0)

    # No more NA corruption
    dv = @data(ones(10_000))
    @test !any(isna(dv))

    PooledDataArray(falses(2), falses(2))
    PooledDataArray(falses(2), trues(2))

    # Test vectorized comparisons work for DataVector's and PooledDataVector's
    @data([1, 2, NA]) .== 1
    @pdata([1, 2, NA]) .== 1
    @data(["1", "2", NA]) .== "1"
    @pdata(["1", "2", NA]) .== "1"

    # Test unique()
    #test_group("unique()")
    # TODO: Restore this
    # dv = DataArray(1:4)
    # dv[4] = NA
    # @test (1 in unique(dv))
    # @test (2 in unique(dv))
    # @test (3 in unique(dv))
    # @test (NA in unique(dv))

    #test_group("find()")
    dv = DataArray([true, false, true])
    @test isequal(find(dv), [1, 3])

    pdv = PooledDataArray([true, false, true])
    @test isequal(find(pdv), [1, 3])

    dv[1] = NA
    @test isequal(find(dv), [3])

    pdv[1] = NA
    @test isequal(find(pdv), [3])
end


text:
### Derivative free optimization tests

using LikelihoodProfiler

# The following NLopt derivative-free algorithms are compared
dfo_algorithms = [
    # good
    (algorithm = :LN_NELDERMEAD, skip = []), # Nelder Mead
    # unstable
    (algorithm = :LN_PRAXIS, skip = [:f_3p_1im, :f_4p_3im, :f_5p_3im]), # "PRAXIS" gradient-free local optimization via the "principal-axis method"
    # errors
    (algorithm = :LN_SBPLX, skip = [:f_3p_1im, :f_3p_1im_dep]), # Subplex (a variant of Nelder-Mead that uses Nelder-Mead on a sequence of subspaces)
    (algorithm = :LN_NEWUOA, skip = [:f_1p, :f_3p_1im, :f_4p_3im, :f_5p_3im, :f_2p_1im, :f_3p_im, :f_2p]),
    (algorithm = :LN_BOBYQA, skip = [:f_1p, :f_3p_1im, :f_1p_ex, :f_4p_3im, :f_5p_3im, :f_2p_1im]), # BOBYQA algorithm for bound constrained optimization without derivatives
    #(algorithm = :LN_COBYLA, skip = []), # Constrained Optimization BY Linear Approximations
]

[test_alg(alg; bounds=(-Inf,Inf)) for alg in dfo_algorithms]

ext: jl


text:
function name(p::Penalty)
    s = replace(string(typeof(p)), "PenaltyFunctions." => "")
    # s = replace(s, r"\{(.*)", "")
    f = fieldnames(typeof(p))
    flength = length(f)
    if flength > 0
        s *= "("
        for (i, field) in enumerate(f)
            s *= "$field = $(getfield(p, field))"
            if i < flength
                s *= ", "
            end
        end
        s *= ")"
    end
    s
end

name(p::ArrayPenalty) = replace(string(typeof(p)), "PenaltyFunctions." => "")

Base.show(io::IO, p::Penalty) = print(io, name(p))
Base.show(io::IO, sp::ScaledArrayPenalty) = print(io, "$(sp.λ) * ($(sp.penalty))")


text:
# Records state from a population.
# Data is recorded as a dict of arrays.

mutable struct RecordedState{S,T,R,IS}
    instance::T
    arrays::Dict{Symbol,Array}
    steps::R
    idx::Int
    state::IS
end

"""
Creates a structure that will hold recorded data.

Recording `obj.a` and `obj.b` at timesteps 100 to 1000:

    r = RecordedState(obj, 100:1000, :a, :b)
    for step in 1:N
        record(r, step)
    end

`steps` can be any `Iterable`, eg. a `Range` (`1:10:10000`)
or `Array` of timesteps.

Accessing recorded data:

    r.arrays[:a]

Accessing timestamps:

    timestamps(r, dt)

"""
function RecordedState(instance, steps, fields...)
    arrays = Dict()
    cols = length(steps)
    for sym in fields
        field = getfield(instance, sym)
        # Initialise an array with the same datatype as the field
        # we want to record, and an extra dimension corresponding
        # to time:
        if isa(field, AbstractArray) # array field
            arr = zeros(eltype(field), size(field)..., cols)
        else # scalar field
            arr = zeros(typeof(field), cols)
        end
        arrays[sym] = arr
    end
    # We use the Val{} trick to store the field symbols in the type
    # parameters. This allows the @generated record() function below
    # to specialise for a specific list of fields.
    state = start(steps)
    return RecordedState{Val{fields}, typeof(instance), typeof(steps), typeof(state)}(instance, arrays, steps, 0, state)
end

"""
Resets recorded data.
"""
function reset!(rec::RecordedState)
    rec.idx = 0
    rec.state = start(rec.steps)
    for arr in values(rec.arrays)
        arr .= zero(eltype(arr))
    end
end

"""
Returns the timestamps corresponding to the recorded timesteps.
"""
timestamps(rec::RecordedState, dt) = rec.steps * dt

quoted(expr) = Expr(:quote, expr)

# Helper functions to efficiently perform the copy:
#   a[:,...,:,i] = b[:,...,:]
function copyslice(a::AbstractArray, b::AbstractArray, i)
    R = CartesianRange(size(b))
    @inbounds for I in R
        a[I,i] = b[I]
    end
end
function copyslice(a::AbstractArray{T}, b::T, i) where {T}
    a[i] = b
end

"""
Signals that there is new data to be recorded for timestep `step`.
"""
@generated function record!(data::RecordedState{Val{syms}}, step) where {syms}
    record_statements = []
    for var in syms
        push!(record_statements, quote
            copyslice(data.arrays[$(QuoteNode(var))], data.instance.$var, data.idx)
        end)
    end
    func = quote
        $(Expr(:meta, :inline))
        $(Expr(:meta, :fastmath))
        if !done(data.steps, data.state)
            next_step, next_state = next(data.steps, data.state)
            if step == next_step
                data.idx += 1
                $(record_statements...)
                data.state = next_state
            end
        end
        nothing
    end
    # println(func)
    func
end


text:
module LeafletJS

    import Base.show
    import GeoJSON, JSON
    import Colors
    import UUIDs

    export LeafletMap, Layer, Provider

    include("providers.jl")
    include("layer.jl")
    include("render.jl")

end # module


text:
using DeepForest, Flux, ProgressMeter, Plots, StatsPlots, Random, CSV, DataFrames
using StatsBase: mode, mean
gr(fmt=:svg, size=(1500, 1000))

data = CSV.File("wdbc.data", header=0) |> DataFrame |> Array
x = Float32.(data[:, 3:end]) |> transpose |> Array
y = Int.(data[:, 2] .== "M") .+ 1

perm = shuffle(Array(1:size(x, 2)))
x = x[:, perm]
y = y[perm]

model = Forest(100, 2, 10, 30, 0.25)
forest_train!(2500, model, x, y)

ŷ = predict(model, x);
@show mean(ŷ .== y)

imp = importance(model, x)
plot_importance(imp, x)


text:
function get_indecies(k, n)
    a = collect(1:k-1)
    append!(a, collect(k+1:n))
    return a
end

function update_W!(W, A, B, r, m)
    for k = 1:r
        indecies = get_indecies(k, r)
        ΣWB = sum(i -> W[:, i] * B[i, k], indecies)
        W[:, k] = maximum(hcat(zeros(m, 1), (A[:, k] - ΣWB) / B[k, k]), 2)
    end
end

function update_H!(H, C, D, r, n)
    for k = 1:r
        indecies = get_indecies(k, r)
        ΣDH = sum(i -> D[k, i] * H[i, :], indecies)
        H[k, :] = maximum(hcat(zeros(n, 1), (C[k, :] - ΣDH) / D[k, k]), 2)'
    end
end

function hierarchical_als(X, W, H, r)
    (m, n) = size(X)
    while stopping_criteria
        A = X * H'
        B = H * H'
        update_W!(W, A, B, r, m)
        C = W' * X
        D = W' * W
        update_H!(H, C, D, r, n)
    end
end


text:
using Documenter, GpABC

makedocs(format=:latex)


text:
using Test
using MPI

MPI.Init()
comm = MPI.COMM_WORLD
nnodes = MPI.Comm_size(comm)
ndims = 3
reorder = 1
periods = [0 1 0]
dims = [0 0 0]
MPI.Dims_create!(nnodes, dims)

@test prod(dims) == nnodes

cperiods = Cint.(periods[:])
cdims = Cint.(dims[:])
comm_cart = MPI.Cart_create(comm, ndims, cdims, cperiods, reorder)
@test MPI.Comm_size(comm_cart) == nnodes

comm_cart2 = MPI.Cart_create(comm, dims, periods, reorder)
@test MPI.Comm_size(comm_cart2) == nnodes

comm_sub1 = MPI.Cart_sub(comm_cart, (false, true, true))
@test MPI.Comm_size(comm_sub1) == div(nnodes, dims[1])

comm_sub2 = MPI.Cart_sub(comm_cart, (true, false, true))
@test MPI.Comm_size(comm_sub2) == div(nnodes, dims[2])

comm_sub3 = MPI.Cart_sub(comm_cart, (true, true, false))
@test MPI.Comm_size(comm_sub3) == div(nnodes, dims[3])

comm_cart = comm_cart2 = comm_sub1 = comm_sub2 = comm_sub3 = nothing

GC.gc()

MPI.Finalize()
@test MPI.Finalized()


text:
export Population, Chromosome, imChromosome, Parameters, indiv_result_type, pop_result_type, elt_type
export CompositionalCircuit, CompGate, CompCircuit, Need, CC, Cc, Cg, Need

using Distributions
using DataFrames
using Random
using QuadGK
using CSV
using Statistics
using Distributed
using Printf
import Random.seed!
export MyInt, DIST_TYPE, Population, PopVect, IPopulation, IntRange, Circuit

if !@isdefined(MyInt)  # MyInt should be defined in CGP.jl
  const MyInt = UInt16
end
if !@isdefined(MyFunc)  # MyFunc should be defined in CGP.jl
  const MyFunc = UInt128
end
const elt_type = Union{Int64,MyInt,String,Float64}
const DIST_TYPE = Dict{Any,Float64}
#const IDIST_TYPE = Dict{Int64,Float64}
#const SDIST_TYPE = Dict{String,Float64}
const Population = Union{Vector{Int64},Vector{MyInt},Vector{String},Vector{Float64}}
const PopVect = Union{Vector{Vector{Int64}},Vector{Vector{MyInt}},Vector{Vector{String}},Vector{Vector{Float64}}}
#const IPopulation = Array{Int64,1}
#const MIPopulation = Array{MyInt,1}
#const SPopulation = Array{String,1}
const FPopulation = Array{Float64,1}
#const CPopulation = Vector{Chromosome}   # population of chromosomes
const IntRange = Union{Integer, Bool, AbstractRange{Int64}, AbstractRange{Bool}}
const FloatRange = Union{Float64,AbstractRange{Float64}}
const Goal =  Vector{MyInt}
const GoalList = Vector{Goal}
#Ones = 0x0f
Ones = 0xffff
#Ones = 0xffffffff

mutable struct Parameters
    mu::Integer
    lambda::Integer
    mutrate::Real
    targetfitness::Real
    numinputs::Integer
    numoutputs::Integer
    nodearity::Integer
    numinteriors::Integer   # Number of gates
    numlevelsback::Integer  # For LinCircuits, number of registers
end
abstract type Node end
mutable struct Func
    func::Function
    arity::Integer
    name::AbstractString
end
mutable struct InputNode <: Node
    index::Integer
    active::Bool
    cache::MyInt
end
mutable struct InteriorNode <: Node
    func::Func
    inputs::Vector{Integer}
    active::Bool
    cache::MyInt
end
mutable struct OutputNode <: Node
    input::Integer
end
PredType = Int64
mutable struct Chromosome
    params::Parameters
    inputs::Vector{InputNode}
    interiors::Vector{InteriorNode}
    outputs::Vector{OutputNode}
    fitness::Float64
    robustness::Union{Float64,PredType}
end
mutable struct LinCircuit
    # Here, numinteriors represents number of gates, numlevelsback represents number of computational registers
    # the First numoutputs registers are the output registers
    circuit_vects::Vector{Vector{MyInt}}
    params::Parameters
end
Circuit = Union{Chromosome,LinCircuit}

# indiv_result_type is used to parallelize population evolution in Indiv_evolution.jl
mutable struct indiv_result_type
  numinputs::Int64
  numoutputs::Int64
  numints::Int64
  levelsback::Int64
  ngoals::Int64
  hamming_sel::Bool
  #avgfitness::Bool
  #robust_sel::Bool
  active_only::Bool
  maxsteps::Int64
  #gl_reps::Int64
  #fault_tol::Bool
  #fit_limit::Float64
  hdist::Float64
  steps::Int64
  same::Int64
  worse::Int64
  better::Int64
  nactive::Int64
  redundancy::Float64
  complexity::Float64
  degeneracy::Float64
  sdegeneracy::Float64
end

# indiv_result_type is used to parallelize population evolution in Indiv_evolution.jl
mutable struct geno_pheno_result_type
  goallist::Vector{Vector{MyInt}}
  numinputs::Int64
  numoutputs::Int64
  numints::Int64
  levelsback::Int64
  ngoals::Int64
  hamming_sel::Bool
  active_only::Bool
  maxsteps::Int64
  gl_reps::Int64
  steps::Int64
  logsteps::Float64
  avgfit::Float64
  nactive::Int64
  redund::Float64
  complex::Float64
  gb_complex::Float64
  degen::Float64
  #gb_degen::Float64
  sdegen::Float64
  f_mutinf::Float64
  mutrobust::Float64
  evolvability::Float64
end

# env_result_type is used to parallelize population evolution in Env_evolution.jl
mutable struct env_result_type
  numinputs::Int64
  numoutputs::Int64
  numints::Int64
  levelsback::Int64
  ngoals::Int64
  hamming_sel::Bool
  #robust_sel::Bool
  active_only::Bool
  maxsteps::Int64
  gl_reps::Int64
  num_flip_bits::Int64
  perturb_goal::Bool
  avgfitness::Bool
  perm_heuristic::Bool
  fit_limit::Float64
  steps::Int64
  avgfit::Float64   # average fitness of final chromosome
  same::Int64
  worse::Int64
  better::Int64
  nactive::Int64
  redundancy::Float64
  complexity::Float64
  degeneracy::Float64
  sdegeneracy::Float64
end

# pop_result_type is used to parallelize population evolution in Pop_evolution.jl
mutable struct pop_result_type
  numinputs::Int64
  numoutputs::Int64
  numints::Int64
  levelsback::Int64
  ngoals::Int64
  popsize::Int64
  max_pop_gens::Int64
  max_indiv_steps::Int64
  hamming_sel::Bool
  robust_sel::Bool
  all_max_sel::Bool
  active_only::Bool
  maxfit::Float64
  maxrobust::Float64
  nactive::Int64
  redundancy::Float64
  complexity::Float64
  degeneracy::Float64
  sdegeneracy::Float64
end

# pop_result_type is used to parallelize population evolution in Inf_alleles.jl
mutable struct inf_alleles_result_type
  numinputs::Int64
  numoutputs::Int64
  numints::Int64
  levsback::Int64
  gl::Vector{Goal}
  popsize::Int64
  max_pop_gens::Int64
  tourn_size::Int64
  func_evals::Int64
  fitness::Float64
  gen_finished::Int64
end

mutable struct evo_result_type
  goal::Goal
  nchromes::Int64
  numinputs::Int64
  numoutputs::Int64
  numints::Int64
  levelsback::Int64
  nrepeats::Int64
  maxsteps::Int64
  all_count::Int64
  evolvable_count::Int64
  evo_diff_count::Int64
end

mutable struct evo_pairs_type
  source_g::Goal
  dest_g::Goal
  numinputs::Int64
  numoutputs::Int64
  numints::Int64
  levelsback::Int64
  maxsteps::Int64
  hamming_dist::Float64
  source_g_count::Int64
  dest_g_count::Int64
  steps::Int64
end
#=
abstract type CC end
mutable struct Cc <: CC
  circuits::Vector{CC}
  inputs::Vector{Int64}
  numoutputs::Int64
  Cc( circuits::Vector{CC}, inputs::Vector{Int64}, numoutputs::Int64 ) = ( x=new(); x.inputs=inputs; x.numoutput=numoutputs; x.circuits=x )
  Cc( cv::Vector{CC}, inputs::Vector{Int64}, numoutputs::Int64 ) = fCc( new(), cv, inputs, numoutputs )
end

function fCc( c::Cc, circuits::Vector{CC}, inputs::Vector{Int64}, numoutputs::Int64 )
  c.circuits=circuits
  c.inputs=inputs
  c.numoutputs=numoutputs
  return c
end

mutable struct Cg <: CC
  func::Func
  inputs::Vector{Int64}
  numoutputs::Int64
end

mutable struct Need
  goal::Goal
  circuit_index::Int64   # If nonzero, the index of the element of all_circuits that meets the need
  inputs::Vector{Int64}
end
=#
abstract type CC end

mutable struct Cc <: CC
  circuits::Vector{CC}
  inputs::Vector{Int64}
  numoutputs::Int64
  output::Goal    # cached output when applied to standard for length(inputs)
  key::String     # dictionary key
  cache_valid::Int64      # non-zero implies cache (which is output) is valid
  Cc( circuits::Vector{CC}, inputs::Vector{Int64}, numoutputs::Int64, output::Goal, key::String, cache_valid::Int64 ) = ( x=new(); x.inputs=inputs; x.numoutputs=numoutputs; x.circuits=circuits; x.output=output; x.key=key; x.cache_valid=x.cache_valid )
  Cc( cv::Vector{CC}, inputs::Vector{Int64}, numoutputs::Int64, output::Goal, key::String, cache_valid::Int64 ) = fCc( new(), cv, inputs, numoutputs, output::Goal, key::String, cache_valid::Int64 )
end

function fCc( c::Cc, circuits::Vector{CC}, inputs::Vector{Int64}, numoutputs::Int64, output::Goal, key::String, cache_valid::Int64 )
  c.circuits=circuits
  c.inputs=inputs
  c.numoutputs=numoutputs
  c.output = output  # use as a cache
  c.key = key
  c.cache_valid = cache_valid
  return c
end

mutable struct Cg <: CC
  func::Func
  inputs::Vector{Int64}
  numoutputs::Int64
  output::Goal
  key::String
end

mutable struct Need
  goal::Goal    # dictionary key
  inputs::Vector{Int64}
  circuit_key::String   # If not "", the key of the all_circuits dictionary that meets the need
  best_fitness::Float64
end


text:
module Ops21_gpu

include("activation.jl")
include("mmul.jl")

end


text:
using RadioTools
using PyPlot

data = readdlm("hydrostatic.bak")#("../../Tesis/data/VALC/hydrostatic.bak")

n2  = data[:,5]
ν2  = 12e9
T2  = data[:,3]
I02 = 5e-13

Iν  = intensity(n2, ν2, T2, I02)
loglog(Iν)
title("Intensity at 12GHz")
xlabel("Step")
ylabel("Intensity")

show()

Tb = Iν * c^2 / (2 * k_b * ν2^2)

loglog(data[:,2],Tb)
title("Emisión total")
ylabel("Tb (K)")
xlabel("Height over photosphere (Km)")
#yscale("log")
ylim(1e3,1e5)

show()


text:
include("algebra.jl")
include("CachedOperator.jl")
include("FiniteOperator.jl")
include("OperatorLayout.jl")
include("PartialInverseOperator.jl")
include("InterlaceOperator.jl")
include("OperatorFunction.jl")


text:
# Arithmetic

@simple_unary(Base.abs)
@simple_unary(Base.exp)
@simple_unary(Base.log)
@simple_unary(Base.log10)
@simple_unary(Base.log2)
@simple_unary(Base.sqrt)
@simple_unary(Base.cbrt)
@simple_unary(Base.sign)
@simple_unary_self_inverse(Base.:-)
@simple_unary_self_inverse(Base.inv)

@simple_binary(Base.:+)
@simple_binary(Base.:-)
@simple_binary(Base.:*)
@simple_binary(Base.:/)
@simple_binary(Base.:^)
@simple_binary(Base.min)
@simple_binary(Base.max)

# Logical

@simple_unary_self_inverse(Base.:!)

@simple_binary(Base.:>)
@simple_binary(Base.:<)
@simple_binary(Base.:>=)
@simple_binary(Base.:<=)

# Linear algebra

@simple_binary(LinearAlgebra.dot)


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

__precompile__(true)

module BAT

using Base.Threads

# import Base.Math.JuliaLibm

using Dates
using Distributed
using LinearAlgebra
using Markdown
using Printf
using Random
using Statistics

using ArgCheck
using ArraysOfArrays
using Clustering
using Colors
using Distributions
using DocStringExtensions
using DoubleFloats
using ElasticArrays
using FFTW
using FillArrays
using ForwardDiffPullbacks
using IntervalSets
using KernelDensity
using LaTeXStrings
using ParallelProcessingTools
using Parameters
using PDMats
using PositiveFactorizations
using ProgressMeter
using RecipesBase
using Requires
using StaticArrays
using StatsBase
using StructArrays
using Tables
using ValueShapes

import AdvancedHMC
import ChainRulesCore
import DiffResults
import DistributionsAD
import EmpiricalDistributions
import FiniteDiff
import ForwardDiff
import HypothesisTests
import Measurements
import NamedArrays
import NLSolversBase
import Optim
import Random123
import Sobol
import TypedTables
import Zygote
import ZygoteRules

using ChainRulesCore: AbstractTangent, Tangent, NoTangent, ZeroTangent

#for AHMI
using DataStructures
using QuadGK

#for Space Partitioning
import CPUTime

include("utils/utils.jl")
include("rngs/rngs.jl")
include("distributions/distributions.jl")
include("variates/variates.jl")
include("transforms/transforms.jl")
include("densities/densities.jl")
include("algotypes/algotypes.jl")
include("initvals/initvals.jl")
include("statistics/statistics.jl")
include("differentiation/differentiation.jl")
include("optimization/optimization.jl")
include("samplers/samplers.jl")
include("integration/integration.jl")
include("algodefaults/algodefaults.jl")
include("io/io.jl")
include("plotting/plotting.jl")

include("precompile.jl")

const _PLOTS_MODULE = Ref{Union{Module,Nothing}}(nothing)
_plots_module() = _PLOTS_MODULE[]

function __init__()
    @require Cuba = "8a292aeb-7a57-582c-b821-06e4c11590b1" include("integration/cuba_integration.jl")
    @require HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" include("io/hdf5_specific.jl")
    @require NestedSamplers = "41ceaf6f-1696-4a54-9b49-2e7a9ec3782e" include("samplers/nested_sampling/ellipsoidal_nested_sampling/ellipsoidal_nested_sampling.jl")
    @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" _PLOTS_MODULE[] = Plots
    @require UltraNest = "6822f173-b0be-4018-9ee2-28bf56348d09" include("samplers/nested_sampling/ultranest.jl")
end

end # module


text:
ENV["GKSwstype"] = "nul"  # For StatsPlots

import Pkg; Pkg.activate(joinpath(@__DIR__, "../../../"))
using CytofDensityEstimation
const CDE = CytofDensityEstimation
using CSV, DataFrames
using LaTeXStrings
using Distributions
using StatsPlots
import Random
using RCall
using BSON
using StatsFuns
@rimport stats as rstats
VD = Vector{Dict{Symbol, Any}}

ygrid = collect(range(-10, 10, length=200))  # TODO: length=200
default_plotsize = (400, 400)
flatten = Iterators.flatten

# plotsize = (400, 400)
# Plots.scalefontsizes()
# Plots.scalefontsizes(1.5)

# Plot observed data histogram.
function plot_observed_hist(yC, yT, imgdir; bins=:auto, binsT=:auto,
                            alpha=0.3, digits=5, plotsize=default_plotsize,
                            legendpos=:topleft)
  binsC = bins
  binsT == nothing && (binsT = bins)
  histogram(yC[isfinite.(yC)], bins=binsC, label=L"\tilde y_C", alpha=alpha,
            color=:blue, linealpha=0, normalize=true)
  histogram!(yT[isfinite.(yT)], bins=binsT, label=L"\tilde y_T", alpha=alpha,
             color=:red, linealpha=0, normalize=true)
  p0C = round(mean(isinf.(yC)), digits=digits)
  p0T = round(mean(isinf.(yT)), digits=digits)
  # title!("prop. 0 in C: $(p0C) | prop. -∞ in T: $(p0T)", titlefont=font(10))
  # xlabel!(L"\tilde y_i")
  ylabel!("density")
  plot!(size=plotsize, legend=:topleft)
  savefig(joinpath(imgdir, "data-hist.pdf"))
  closeall()
end

function postprocess(chain, laststate, summarystats, yC, yT, imgdir;
                     bw_postpred=0.2, density_legend_pos=:best,
                     simdata=nothing, ygrid=ygrid)
  # Print summary statistics.
  CDE.Model.printsummary(chain, summarystats)
  flush(stdout)

  # Print KS Statistic.
  ks_fit = rstats.ks_test(yC, yT)
  println("KS-test p-value: ", ks_fit["p.value"][1])
  flush(stdout)

  # Plot results.
  @time CDE.Model.plotpostsummary(chain, summarystats, yC, yT, imgdir,
                                  bw_postpred=bw_postpred, ygrid=ygrid,
                                  density_legend_pos=density_legend_pos,
                                  simdata=simdata)
  flush(stdout)
end

function postprocess(chain0, chain1, data, imgdir, awsbucket;
                     simdata=nothing, bw_postpred=0.20,
                     density_legend_pos=:best, binsC=:auto, binsT=:auto,
                     ygrid=ygrid, p=nothing)
  mkpath(imgdir)
  CDE.Util.redirect_stdout_to_file(joinpath(imgdir, "bf.txt")) do
    println("Start merging results...")

    # Plot observed data.
    plot_observed_hist(data.yC, data.yT, imgdir,
                       alpha=0.3, digits=5, legendpos=:topleft)

    # Compute Bayes factor in favor of M1.
    p == nothing && (p = CDE.Model.group(:p, chain0)[1])
    @time lbf = CDE.Model.log_bayes_factor(data, VD(chain0[1]), VD(chain1[1]))
    pm1 = CDE.Model.posterior_prob1(lbf, logpriorodds=logit(p))
    println("Log Bayes Factor: $(lbf) | P(M1|y): $(pm1)"); flush(stdout)

    # Print KS Statistic.
    ks_fit = rstats.ks_test(data.yC, data.yT)
    println("KS-test p-value: ", ks_fit["p.value"][1])
    flush(stdout)

    # Plot marginal posterior density
    println("Plot average posterior predictive....")
    @time CDE.Model.plot_posterior_predictive(data.yC, data.yT, chain0, chain1, pm1,
                                              imgdir, bw_postpred=bw_postpred,
                                              density_legend_pos=density_legend_pos,
                                              ygrid=ygrid, plotsize=(300, 300),
                                              simdata=simdata, lw=.5, ls=:solid,
                                              binsC=binsC, binsT=binsT,
                                              digits=5, fontsize=7)
    flush(stdout)

    CDE.Model.plot_gamma(data.yC, data.yT, chain0, chain1, pm1, imgdir, xrotation=nothing)

    # DIC
    @time dic0, dic1 = CDE.dic(chain0[1], data), CDE.dic(chain1[1], data)
    println("(DIC0, DIC1): ($(round(dic0, digits=3)), $(round(dic1, digits=3)))")

    @time dic_average = CDE.dic(chain0[1], chain1[1], pm1, data)
    println("DIC average: $(round(dic_average, digits=3))")

    flush(stdout)
  end

  # Send results
  if awsbucket != nothing
    CDE.Util.s3sync(from=imgdir, to=awsbucket, tags=`--exclude '*.nfs'`)
  end
end

function defaults(yC, yT, K; seed=nothing)
  seed == nothing || Random.seed!(seed)

  data = CDE.Model.Data(yC, yT)
  y = [data.yC_finite; data.yT_finite]
  prior = CDE.Model.Prior(K,
                          mu=Normal(mean(y), (maximum(y) - minimum(y))/2),
                          nu=LogNormal(3, .5),
                          p=Beta(100, 100),
                          psi=Normal(-1, .5),
                          a_omega=2.5,
                          tau=Gamma(.5, 1.),
                          eta=Dirichlet(K, 1/K))
  state = CDE.Model.State(data, prior)
  tuners = CDE.Model.Tuners(K, 1.0)
  return (state=state, data=data, prior=prior, tuners=tuners)
end

function run(config)
  resultsdir = config[:resultsdir]
  mkpath(resultsdir)
  path_to_log = joinpath(resultsdir, "log.txt")
  CDE.Util.redirect_stdout_to_file(path_to_log) do
    _run(config)
  end
end

function _run(config)
  println("pid: $(getpid())"); flush(stdout)

  resultsdir = config[:resultsdir]
  awsbucket = config[:awsbucket]
  yC = config[:yC]
  yT = config[:yT]
  K = config[:K]
  nsamps = config[:nsamps]
  nburn = config[:nburn]
  thin = config[:thin]

  # Make image path id needed.
  imgdir = joinpath(resultsdir, "img")
  mkpath(imgdir)

  # Create initial state, data, prior, and tuners objects.
  state, data, prior, tuners = defaults(yC, yT, K, seed=0)
  state.p = config[:p]
  state.beta = config[:beta]
  println("Priors:")
  foreach(fn -> println(fn, " => ", getfield(prior, fn)),
          fieldnames(CDE.Model.Prior))

  # Plot data.
  plot_observed_hist(yC, yT, imgdir,
                     alpha=0.3, digits=5, legendpos=:topleft)

  # Run analysis.
  println("Run Chain ..."); flush(stdout)
  @time chain, laststate, summarystats = CDE.fit(
      state, data, prior, tuners, nsamps=[nsamps], fix=[:p, :beta],
      nburn=nburn, thin=thin, rep_aux=1, flags=CDE.Model.Flag[])
  flush(stdout)

  # Save results
  BSON.bson("$(resultsdir)/results.bson",
            Dict(:chain=>chain, :laststate=>laststate, :data=>data,
                 :prior=>prior, :summarystats=>summarystats))

  # Load data:
  out = BSON.load("$(resultsdir)/results.bson")

  # Post process
  postprocess(out[:chain], out[:laststate], out[:summarystats], out[:data].yC,
              out[:data].yT, imgdir, bw_postpred=0.3,
              ygrid=ygrid,
              density_legend_pos=:topleft)

  # Send results
  if awsbucket != nothing
    CDE.Util.s3sync(from=resultsdir, to=awsbucket, tags=`--exclude '*.nfs'`)
  end

  # Print done.
  println("Done!"); flush(stdout)
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
using Revise
using Test
using VulkanCore
using LinearAlgebra
using GeometryTypes
using Lava
using Lava: @class, @scope, @autodestroy
include("TestAPIs.jl")

# include("class_method.jl")
# include("class_method_benchmark.jl")
#
# println("Testing...")
# @test testClassMethod() == 0
# @test benchmarkClassMethod() == 0

# include("AutoDestroyTest.jl")


text:
@testset "Expr conversion" begin
    @testset "Quote nodes" begin
        @test parseall(Expr, ":(a)", rule=:atom) == QuoteNode(:a)
        @test parseall(Expr, ":(:a)", rule=:atom) == Expr(:quote, QuoteNode(:a))
        @test parseall(Expr, ":(1+2)", rule=:atom) == Expr(:quote, Expr(:call, :+, 1, 2))
        # Compatibility hack for VERSION >= v"1.4"
        # https://github.com/JuliaLang/julia/pull/34077
        @test parseall(Expr, ":true", rule=:atom) == Expr(:quote, true)
    end

    @testset "Short form function line numbers" begin
        # A block is added to hold the line number node
        @test parseall(Expr, "f() = xs", rule=:statement) ==
            Expr(:(=),
                 Expr(:call, :f),
                 Expr(:block,
                      LineNumberNode(1, :none),
                      :xs))
        # flisp parser quirk: In a for loop the block is not added, despite
        # this defining a short-form function.
        @test parseall(Expr, "for f() = xs\nend", rule=:statement) ==
            Expr(:for,
                 Expr(:(=), Expr(:call, :f), :xs),
                 Expr(:block))
    end
end


text:
# Optical constants of water
#  wavelen.  ext. fact.   index of refr.
#  λ(nm)      k           n
const WATERCONST =
[   200     1.10e-7     1.396;
    225     4.90e-8     1.373;
    250     3.35e-8     1.362;
    275     2.35e-8     1.354;
    300     1.60e-8     1.349;
    325     1.08e-8     1.346;
    350     6.50e-9     1.343;
    375     3.50e-9     1.341;
    400     1.86e-9     1.339;
    425     1.30e-9     1.338;
    450     1.02e-9     1.337;
    475     9.35e-10    1.336;
    500     1.00e-9     1.335;
    525     1.32e-9     1.334;
    550     1.96e-9     1.333;
    575     3.60e-9     1.333;
    600     1.09e-8     1.332;
    625     1.39e-8     1.332;
    650     1.64e-8     1.331;
    675     2.23e-8     1.331;
    700     3.35e-8     1.331;
    725     9.15e-8     1.330;
    750     1.56e-7     1.330;
    775     1.48e-7     1.330;
    800     1.25e-7     1.329;
    825     1.82e-7     1.329;
    850     2.93e-7     1.329;
    875     3.91e-7     1.328;
    900     4.86e-7     1.328;
    925     1.06e-6     1.328;
    950     2.93e-6     1.327;
    975     3.48e-6     1.327;
   1000     2.89e-6     1.327;
   1200     9.89e-6     1.324;
   1400     1.38e-4     1.321;
   1600     8.55e-5     1.317;
   1800     1.15e-4     1.312;
   2000     1.10e-3     1.306;
   2200     2.89e-4     1.296;
   2400     9.56e-4     1.279;
   2600     3.17e-3     1.242;
   2650     6.70e-3     1.219;
   2700     0.01900     1.188;
   2750     0.05900     1.157;
   2800     0.11500     1.142;
   2850     0.18500     1.149;
   2900     0.26800     1.201;
   2950     0.29800     1.292;
   3000     0.27200     1.371;
   3050     0.24000     1.426;
   3100     0.19200     1.467;
   3150     0.13500     1.483;
   3200     0.09240     1.478;
   3250     0.06100     1.467;
   3300     0.03680     1.450;
   3350     0.02610     1.432;
   3400     0.01950     1.420;
   3450     0.01320     1.410;
   3500     0.00940     1.400;
   3600     0.00515     1.385;
   3700     0.00360     1.374;
   3800     0.00340     1.364;
   3900     0.00380     1.357;
   4000     0.00460     1.351;
   4100     0.00562     1.346;
   4200     0.00688     1.342;
   4300     0.00845     1.338;
   4400     0.01030     1.334;
   4500     0.01340     1.332;
   4600     0.01470     1.330;
   4700     0.01570     1.330;
   4800     0.01500     1.330;
   4900     0.01370     1.328;
   5000     0.01240     1.325;
   5100     0.01110     1.322;
   5200     0.01010     1.317;
   5300     0.00980     1.312;
   5400     0.01030     1.305;
   5500     0.01160     1.298;
   5600     0.01420     1.289;
   5700     0.02030     1.277;
   5800     0.03300     1.262;
   5900     0.06220     1.248;
   6000     0.10700     1.265;
   6100     0.13100     1.319;
   6200     0.08800     1.363;
   6300     0.05700     1.357;
   6400     0.04490     1.347;
   6500     0.03920     1.339;
   6600     0.03560     1.334;
   6700     0.03370     1.329;
   6800     0.03270     1.324;
   6900     0.03220     1.321;
   7000     0.03200     1.317;
   7100     0.03200     1.314;
   7200     0.03210     1.312;
   7300     0.03220     1.309;
   7400     0.03240     1.307;
   7500     0.03260     1.304;
   7600     0.03280     1.302;
   7700     0.03310     1.299;
   7800     0.03350     1.297;
   7900     0.03390     1.294;
   8000     0.03430     1.291;
   8200     0.03510     1.286;
   8400     0.03610     1.281;
   8600     0.03720     1.275;
   8800     0.03850     1.269;
   9000     0.03990     1.262;
   9200     0.04150     1.255;
   9400     0.04330     1.247;
   9600     0.04540     1.239;
   9800     0.04790     1.229;
  10000     0.05080     1.218;
  10500     0.06620     1.185;
  11000     0.09680     1.153;
  11500     0.14200     1.126;
  12000     0.19900     1.111;
  12500     0.25900     1.123;
  13000     0.30500     1.146;
  13500     0.34300     1.177;
  14000     0.37000     1.210;
  14500     0.38800     1.241;
  15000     0.40200     1.270;
  15500     0.41400     1.297;
  16000     0.42200     1.325;
  16500     0.42800     1.351;
  17000     0.42900     1.376;
  17500     0.42900     1.401;
  18000     0.42600     1.423;
  18500     0.42100     1.443;
  19000     0.41400     1.461;
  19500     0.40400     1.476;
  20000     0.39300     1.480;
  21000     0.38200     1.487;
  22000     0.37300     1.500;
  23000     0.36700     1.511;
  24000     0.36100     1.521;
  25000     0.35600     1.531;
  26000     0.35000     1.539;
  27000     0.34400     1.545;
  28000     0.33800     1.549;
  29000     0.33300     1.551;
  30000     0.32800     1.551;
  32000     0.32400     1.546;
  34000     0.32900     1.536;
  36000     0.34300     1.527;
  38000     0.36100     1.522;
  40000     0.38500     1.519;
  42000     0.40900     1.522;
  44000     0.43600     1.530;
  46000     0.46200     1.541;
  48000     0.48800     1.555;
  50000     0.51400     1.587;
  60000     0.58700     1.703;
  70000     0.57600     1.821;
  80000     0.54700     1.886;
  90000     0.53600     1.924;
 100000     0.53200     1.957;
 110000     0.53100     1.966;
 120000     0.52600     2.004;
 130000     0.51400     2.036;
 140000     0.50000     2.056;
 150000     0.49500     2.069;
 160000     0.49600     2.081;
 170000     0.49700     2.094;
 180000     0.49900     2.107;
 190000     0.50100     2.119;
 200000     0.50400     2.130]


text:
using HTTP

escape(s::String) = replace(s, "\\" => "\\\\")

function download_schema(schema_git_tag)
    schema_path = joinpath(@__DIR__, "$(schema_git_tag).json")
    r = HTTP.request("GET", "https://raw.githubusercontent.com/odow/MathOptFormat/$(schema_git_tag)/mof.schema.json")
    if r.status == 200
        open(schema_path, "w") do io
            write(io, String(r.body))
        end
        open(joinpath(@__DIR__, "deps.jl"), "w") do io
            write(io, "const SCHEMA_PATH = \"$(escape(schema_path))\"\n")
        end
    else
        error("Unable to download the latest MathOptFormat schema.\n" *
              "HTTP status: $(r.status)\n" *
              "HTTP body: $(String(r.body))")
    end
end

# Update this tag whenever github.com/odow/MathOptFormat releases a new update
# to the schema (providing that any code in this package is also updated).
const SCHEMA_GIT_TAG = "v0.2.1"

download_schema(SCHEMA_GIT_TAG)


text:
import Dates

@testset "T5.jl" begin
    filename = "data\\test.t5"
    df = DDR2import.T5.read(filename)
    @test df[2,:FLIGHT_ID] == 213765625
    @test df[2,:SECTOR_NAME] == "EHAAFIR"
    @test Dates.year(df[2,:ENTRY_DATETIME]) == 2017
    @test Dates.day(df[2,:ENTRY_DATETIME]) == 14
    @test Dates.minute(df[2,:ENTRY_DATETIME]) == 49
    @test Dates.month(df[2,:EXIT_DATETIME]) == 12
    @test Dates.hour(df[2,:EXIT_DATETIME]) == 22
    @test Dates.second(df[2,:EXIT_DATETIME]) == 0
    @test df[2,:ENTRY_FL] == 0.0
    @test df[2,:EXIT_FL] == 255.0
    @test df[2,:ENTRY_SEGMENT_NAME] == "EHAM_*AM80"
    @test df[2,:EXIT_SEGMENT_NAME] == "!AAEb_SONEB"
    @test df[2,:TOT_DISTANCE_IN_SEGMENT_M] == 84.060844 * 1852.0
    @test df[2,:TOT_TIME_IN_SEGMENT_S] == 840.0
    #LAST
    @test df[4,:SECTOR_NAME] == "EH_DN"
end


text:
## SCENARIO PARTITION
function GAPM(Ind_old, duals; α_val = 0.)
    dualRange = maximum(duals)-minimum(duals)
    MaxDiff = α_val*dualRange # Maximum difference between scenario duals

    NP = maximum(Ind_old)

    Scenarios = [duals Ind_old collect(1:N)]

    Ind,NsubParts = ones(N), 1

    for p in 1:NP
        DualPart = Scenarios[Scenarios[:,2].==p,:]
        DualPart=DualPart[sortperm(DualPart[:, 1]), :]
        sizePart = size(DualPart, 1)
        Ind[Integer(DualPart[1,3])] = NsubParts
        ref = DualPart[1,1]
        for n in 2:sizePart
            # Evaluating distance between duals
            if abs(DualPart[n,1] - ref) > MaxDiff
                NsubParts += 1
                ref = DualPart[n,1] # Reference (1st) value for the partition
            end
            Ind[Integer(DualPart[n,3])] = NsubParts
        end
        NsubParts += 1
    end
        return Integer.(Ind)
end

max_stars_repo_name: gonzalez-alvaro/MultidimensionalGAPM

max_issues_repo_name: gonzalez-alvaro/MultidimensionalGAPM

max_forks_repo_name: gonzalez-alvaro/MultidimensionalGAPM


text:
using Test, DynamicEnergyBudgets, Microclimate, Unitful, Photosynthesis

import DynamicEnergyBudgets: get_environment, temp, apply_environment!

rad = [1000.0, 900.0]u"W*m^-2"
snow = [1.0, 0.9]u"m"
airt = [25.0 24.0
        30.0 29.0]u"°C" .|> u"K"
ws = [1.0 2.0
      3.0 4.0]u"m*s^-1"
rh = [0.7 0.5
      0.8 0.6]
soilt = [20.0 19.0 18.0 17.0 16.0 15.0 14.0 13.0
         12.0 11.0 10.0 9.0  8.0  7.0  6.0  5.0]u"°C" .|> u"K"
soilwp = [-100.0 -90.0 -120.0 -80.0 -110.0 -70.0 -130.0 -60.0
          -200.0 -190.0 -220.0 -280.0 -110.0 -170.0 -230.0 -160.0]u"kPa"
soilwc = [0.3 0.3 0.2 0.3 0.2 0.3 0.2 0.3
          0.4 0.5 0.4 0.5 0.4 0.5 0.4 0.5]

environment = MicroclimPoint{0,Microclimate.LAYERINCREMENTS,Microclimate.NEXTLAYER,Microclimate.LAYERRANGE}(
    rad, snow, airt, rh, ws, soilt, soilwp, soilwc
)

u = [9.0,8.0,7.0]u"mol"

@testset "apply environment to deb vars" begin
    o = DynamicEnergyBudgets.Plant(
        environment=environment,
        environment_start=0u"hr",
       )
    o1, o2 = o.organs
    v1 = o1.vars; v2 = o2.vars
    p1 = o1.params; p2 = o2.params
    av1 = p1.assimilation_pars.vars; av2 = p2.assimilation_pars.vars

    av1.soilwaterpotential = 0.0u"kPa"
    av2.soilwaterpotential = 0.0u"kPa"
    av1.J_L_F == 0.0u"W*m^-2"

    t = 0u"hr"
    apply_environment!(o, o.organs, u, t)
    @test temp(v1) |> u"°C" == 25.0u"°C"
    @test temp(v2) |> u"°C"  == 20.0u"°C"
    @test av1.J_L_F == 1000.0u"W*m^-2" * 4.57u"mol*W^-1*s^-1"
    @test av1.soilwaterpotential == -100.0u"kPa"
    # @test av2.soilwaterpotential == -100.0u"kPa"

    t = 1u"hr"
    apply_environment!(o, o.organs, u, t)
    @test temp(v1) |> u"°C"  == 30.0u"°C"
    @test temp(v2) |> u"°C"  == 12.0u"°C"
    @test av1.J_L_F == 900.0u"W*m^-2" * 4.57u"mol*W^-1*s^-1"
    @test av1.soilwaterpotential == -200.0u"kPa"
    # @test av2.soilwaterpotential == -200.0u"kPa"
end

@testset "apply environment to deb vars" begin
    o = DynamicEnergyBudgets.Plant(
        params=(Params(assimilation_pars=BallBerryPotentialCAssim()), Params()),
        environment=environment,
        environment_start=0u"hr",
    )

    o1, o2 = o.organs;
    v1 = o1.vars;
    p1 = o1.params;
    av = p1.assimilation_pars.vars

    av.tair = 0.0u"°C"
    av.windspeed = 0.0u"m*s^-1"
    av.rh = 0.0
    av.rnet = 0.0u"W*m^-2"
    av.par = 0.0u"W*m^-2" * 0.0u"mol*W^-1*s^-1"
    av.soilmoist = 0.0
    av.swp = 0.0u"kPa"

    t = 0u"hr"
    apply_environment!(o, o.organs, u, t)
    @test av.tair == 25.0u"°C"
    @test av.tleaf > av.tair
    @test temp(v1) == av.tleaf
    @test av.windspeed == 1.0u"m*s^-1"
    @test av.rh == 0.7
    @test av.rnet == 1000.0u"W*m^-2"
    @test av.par == 1000.0u"W*m^-2" * 4.57u"mol*W^-1*s^-1"
    # @test av.soilmoist == 0.2 not used currently
    @test av.swp == -100.0u"kPa"

    t = 1u"hr"
    apply_environment!(o, o.organs, u, t)
    @test av.tair == 30.0u"°C"
    @test av.tleaf > av.tair
    @test temp(v1) == av.tleaf
    @test av.windspeed == 3.0u"m*s^-1"
    @test av.rh == 0.8
    @test av.rnet == 900.0u"W*m^-2"
    @test av.par == 900.0u"W*m^-2" * 4.57u"mol*W^-1*s^-1"
    # @test av.soilmoist == 0.4 not used currently
    @test av.swp == -200.0u"kPa"
end

nothing


text:
@testset "MolecularSpec" begin
    @testset "Atom" begin
        a = Atom(:H, (1, 2, 3))
        b = Atom(:H, (1, 2, 3))
        @test a == b
        @test a === b
        @test_throws MethodError copy(a)
        @test_throws ErrorException Atom(:NonExistentElement, (1, 2, 3))
    end

    @testset "Geometry" begin
        a = Atom(:H, (0.0, 0.0, 0.0))
        b = Atom(:H, (0.0, 0.0, 0.74))
        g = Geometry(a, b)
        g1 = copy(g)
        @test g == g1
        @test length(g) == 2
        @test g[1] == a
        @test g[2] == b
    end
end


text:
# module for testing each model modules

module TestModel
  using Test

  # target modules
  include(joinpath(split(@__FILE__, "test")[1], "src/model/world/world.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/model/world/puddle_world.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/model/world/world_sample.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/model/agent/agent.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/model/agent/puddle_ignore_agent.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/model/agent/dp_policy_agent.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/model/agent/q_agent.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/model/map/map.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/model/object/object.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/model/sensor/sensor.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/model/sensor/robot_sensing_sample.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/model/robot/differential_wheeled_robot/differential_wheeled_robot.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/model/robot/warp_robot/warp_robot.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/model/robot/robot_move_sample.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/model/puddle/puddle.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/model/goal/goal.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/model/uncertainty/movement/kidnap/anime_move_kidnap.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/model/uncertainty/movement/random_noise/anime_move_random_noise.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/model/uncertainty/movement/speed_bias/anime_move_speed_bias.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/model/uncertainty/movement/stuck/anime_move_stuck.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/model/uncertainty/observation/bias/anime_observe_bias.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/model/uncertainty/observation/noise/anime_observe_noise.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/model/uncertainty/observation/occlusion/anime_observe_occlusion.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/model/uncertainty/observation/oversight/anime_observe_oversight.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/model/uncertainty/observation/phantom/anime_observe_phantom.jl"))

  function main()
    @testset "Model" begin
      @testset "World" begin
        w = World(-5.0, 5.0, -5.0, 5.0)
        @test w.x_min == -5.0
        @test w.x_max == 5.0
        @test w.y_min == -5.0
        @test w.y_max == 5.0
        @test_nowarn draw(w)
      end
      @testset "PuddleWorld" begin
        pw = PuddleWorld(-5.0, 5.0, -5.0, 5.0, is_test=true)
        @test pw.x_min == -5.0
        @test pw.x_max == 5.0
        @test pw.y_min == -5.0
        @test pw.y_max == 5.0
        @test length(pw.objects) == 0
        @test pw.delta_time == 0.1
        @test pw.end_time == 30
        @test pw.is_test == true
        @test pw.save_path === nothing
        @test length(pw.puddles) == 0
        @test length(pw.robots) == 0
        @test length(pw.goals) == 0
        append(pw, DifferentialWheeledRobot([1, 1, 1], 0.5, "red",
                                            PuddleIgnoreAgent(), 0.1))
        append(pw, Puddle([0.0, 0.0], [6.0, 6.0], 0.1))
        append(pw, Goal(1.0, 2.0, radius=0.5, value=1.0))
        @test length(pw.puddles) == 1
        @test length(pw.robots) == 1
        @test length(pw.goals) == 1
        @test puddle_depth(pw, [1.0, 1.0, 1.0]) == 0.1
        @test_nowarn draw(pw)
      end
      @testset "WorldSample" begin
        @test_nowarn WorldSample.main(is_test=true)
      end
      @testset "Agent" begin
        a = Agent(0.1, 1.0)
        obsrv = []
        speed, yaw_rate = draw_decision!(a, obsrv)
        @test speed == 0.1
        @test yaw_rate == 1.0
        @test a.time_interval == 0.1
        @test a.estimator === nothing
        @test a.prev_spd == 0.0
        @test a.prev_yr == 0.0
      end
      @testset "PuddleIgnoreAgent" begin
        pia = PuddleIgnoreAgent()
        @test pia.speed == 0.0
        @test pia.yaw_rate == 0.0
        @test pia.delta_time == 0.1
        @test pia.estimator === nothing
        @test pia.prev_spd == 0.0
        @test pia.prev_yr == 0.0
        @test pia.puddle_coef == 100
        @test pia.puddle_depth == 0.0
        @test pia.total_reward == 0.0
        @test pia.in_goal == false
        @test pia.final_value == 0.0
        @test pia.goal === nothing
        @test reward_per_sec(pia) == -1.0
      end
      @testset "DpPolicyAgent" begin
        dpa = DpPolicyAgent()
        @test dpa.speed == 0.0
        @test dpa.yaw_rate == 0.0
        @test dpa.delta_time == 0.1
        @test dpa.estimator === nothing
        @test dpa.prev_spd == 0.0
        @test dpa.prev_yr == 0.0
        @test dpa.puddle_coef == 100
        @test dpa.puddle_depth == 0.0
        @test dpa.total_reward == 0.0
        @test dpa.in_goal == false
        @test dpa.final_value == 0.0
        @test dpa.goal === nothing
        @test dpa.pose_min == [-4.0, -4.0, 0.0]
        @test dpa.pose_max == [4.0, 4.0, 2*pi]
        @test dpa.widths == [0.2, 0.2, pi/18]
        @test dpa.index_nums == round.(Int64, (dpa.pose_max - dpa.pose_min)./dpa.widths)
        @test reward_per_sec(dpa) == -1.0
        @test to_index(dpa, [3.0, 3.0, 0.0]) == [35, 35, 0]
        @test to_index(dpa, [5.0, 5.0, 0.0]) == [39, 39, 0]
        @test_nowarn draw_decision!(dpa, [])
      end
      @testset "QAgent" begin
        qa = QAgent()
        @test qa.speed == 0.0
        @test qa.yaw_rate == 0.0
        @test qa.delta_time == 0.1
        @test qa.estimator === nothing
        @test qa.prev_spd == 0.0
        @test qa.prev_yr == 0.0
        @test qa.puddle_coef == 100
        @test qa.puddle_depth == 0.0
        @test qa.total_reward == 0.0
        @test qa.in_goal == false
        @test qa.final_value == 0.0
        @test qa.goal === nothing
        @test qa.pose_min == [-4.0, -4.0, 0.0]
        @test qa.pose_max == [4.0, 4.0, 2*pi]
        @test qa.widths == [0.2, 0.2, pi/18]
        @test qa.index_nums == round.(Int64, (qa.pose_max - qa.pose_min)./qa.widths)
        @test qa.alpha == 0.5
        @test qa.s === nothing
        @test qa.a === nothing
        @test qa.update_end == false
        @test qa.stuck_timer == 0.0
        @test reward_per_sec(qa) == -1.0
        @test to_index(qa, [3.0, 3.0, 0.0]) == [35, 35, 0]
        @test to_index(qa, [5.0, 5.0, 0.0]) == [39, 39, 0]
        @test_nowarn draw_decision!(qa, [])
      end
      @testset "Object" begin
        o = Object(1.0, 2.0, shape=:circle, color=:red, size=20, id=4)
        @test o.pose == [1.0, 2.0]
        @test o.id == 4
        @test o.size == 20
        @test o.color == :red
        @test o.shape == :circle
      end
      @testset "Map" begin
        m = Map()
        @test m.objects == []
        add_object(m, Object(1.0, 2.0))
        @test m.objects[1].pose == [1.0, 2.0]
      end
      @testset "Sensor" begin
        m = Map()
        add_object(m, Object(2.0, -2.0))
        s = Sensor(m)
        @test s.last_data == []
        @test s.dist_rng == (0.5, 6.0)
        @test s.dir_rng == (-pi/3, pi/3)
        @test visible(s, [0.5, -pi/3]) == true
        @test visible(s, [6.0, pi/3]) == true
        @test visible(s, [0.2, -pi/3]) == false
        @test visible(s, [0.5, pi]) == false
      end
      @testset "RobotSensingSample" begin
        @test_nowarn RobotSensingSample.main(is_test=true)
      end
      @testset "DifferentialWheeledRobot" begin
        a = Agent(0.1, 1.0)
        r = DifferentialWheeledRobot([1, 1, 1], 0.5, "red", a, 0.1)
        @test r.pose == [1, 1, 1]
        @test r.radius == 0.5
        @test r.color == "red"
        @test r.delta_time == 0.1
        @test r.agent.speed == 0.1
        @test r.agent.yaw_rate == 1.0
      end
      @testset "WarpRobot" begin
        a = Agent(0.1, 1.0)
        wr = WarpRobot([1, 1, 1], 0.5, "red", a, 0.1)
        @test wr.pose == [1, 1, 1]
        @test wr.radius == 0.5
        @test wr.color == "red"
        @test wr.delta_time == 0.1
        @test wr.agent.speed == 0.1
        @test wr.agent.yaw_rate == 1.0
        @test wr.is_stuck == false
        @test wr.traj_x[1] == 1
        @test wr.traj_y[1] == 1
      end
      @testset "RobotMoveSample" begin
        @test_nowarn RobotMoveSample.main(is_test=true)
      end
      @testset "Goal" begin
        g = Goal(1.0, 2.0, radius=0.5, value=1.0)
        @test g.pose[1] == 1.0
        @test g.pose[2] == 2.0
        @test g.radius == 0.5
        @test g.value == 1.0
        @test inside(g, [1.2, 2.0]) == true
        @test inside(g, [1.0, 3.0]) == false
        @test_nowarn draw!(g)
      end
      @testset "Puddle" begin
        p = Puddle([0.0, 0.0], [6.0, 6.0], 0.1)
        @test p.lower_left == [0.0, 0.0]
        @test p.upper_right == [6.0, 6.0]
        @test p.depth == 0.1
        @test_nowarn draw!(p)
        @test inside(p, [2.0, 2.0]) == true
        @test inside(p, [-1.0, -2.0]) == false
      end
      @testset "uncertainty" begin
        @testset "Movement" begin
          @testset "Kidnap" begin
            @test_nowarn AnimeMoveKidnap.main(is_test=true)
          end
          @testset "RandomNoise" begin
            @test_nowarn AnimeMoveRandomNoise.main(is_test=true)
          end
          @testset "SpeedBias" begin
            @test_nowarn AnimeMoveSpeedBias.main(is_test=true)
          end
          @testset "Stuck" begin
            @test_nowarn AnimeMoveStuck.main(is_test=true)
          end
        end
        @testset "Observation" begin
          @testset "Bias" begin
            @test_nowarn AnimeObserveBias.main(is_test=true)
          end
          @testset "Noise" begin
            @test_nowarn AnimeObserveNoise.main(is_test=true)
          end
          @testset "Occlusion" begin
            @test_nowarn AnimeObserveOcclusion.main(is_test=true)
          end
          @testset "Oversight" begin
            @test_nowarn AnimeObserveOversight.main(is_test=true)
          end
          @testset "Phantom" begin
            @test_nowarn AnimeObservePhantom.main(is_test=true)
          end
        end
      end
    end
  end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
struct FdHullWhiteSwaptionEngine{Y <: TermStructure, F <: FdmSchemeDescType} <: PricingEngine
  model::HullWhite{AffineModelType, Y}
  tGrid::Int
  xGrid::Int
  dampingSteps::Int
  invEps::Float64
  schemeDesc::FdmSchemeDesc{F}
  ts::Y
end

FdHullWhiteSwaptionEngine(model::HullWhite{AffineModelType, Y},
                          tGrid::Int = 100,
                          xGrid::Int = 100,
                          dampingSteps::Int = 0,
                          invEps::Float64 = 1e-5,
                          schemeDesc::FdmSchemeDesc{F} = FdmSchemeDesc(Douglas())) where {F <: FdmSchemeDescType, Y <: TermStructure} =
                  FdHullWhiteSwaptionEngine{Y, F}(model, tGrid, xGrid, dampingSteps, invEps, schemeDesc, model.ts)
# methods #
function _calculate!(pe::FdHullWhiteSwaptionEngine, swaption::Swaption)
  # 1. Term structure
  ts = pe.ts

  # 2. Mesher
  dc = ts.dc
  refDate = reference_date(ts)
  maturity = year_fraction(dc, refDate, swaption.exercise.dates[end])

  process = OrnsteinUhlenbeckProcess(get_a(pe.model), get_sigma(pe.model))
  shortRateMesher = FdmSimpleProcess1dMesher(pe.xGrid, process, maturity, 1, pe.invEps)
  mesher = FdmMesherComposite(shortRateMesher)

  # 3. Inner Value calculator
  exerciseDates = swaption.exercise.dates
  t2d = Dict{Float64, Date}()
  @simd for i = 1:length(exerciseDates)
    @inbounds t = year_fraction(dc, refDate, exerciseDates[i])
    @inbounds t2d[t] = exerciseDates[i]
  end

  disTs = pe.model.ts
  fwdTs = swaption.swap.iborIndex.ts

  # TODO check that day counts and ref dates match btwn the two term structures
  fwdModel = HullWhite(fwdTs, get_a(pe.model), get_sigma(pe.model))
  calculator = FdmAffineModelSwapInnerValue(pe.model, fwdModel, swaption.swap, t2d, mesher, 1)

  # 4. Step Conditions
  conditions = vanilla_FdmStepConditionComposite(DividendSchedule(), swaption.exercise, mesher, calculator, refDate, dc)

  # 5. Boundary conditions
  boundaries = FdmBoundaryConditionSet()

  # 6. Solver
  solverDesc = FdmSolverDesc(mesher, boundaries, conditions, calculator, maturity, pe.tGrid, pe.dampingSteps)

  solver = FdmHullWhiteSolver(pe.model, solverDesc, pe.schemeDesc)
  swaption.results.value = value_at(solver, 0.0)
end


text:
# adds a triangle to mesh where a b c d are the relative index coordinates
macro addtri(a,b,c,d)
    return esc(quote
               P2 = center;
               P1 = zeros(3); P1[fixdim[l]]=fixdim_val[l]; P3 = copy(P1);
               P1[freedim[1]]=grid[ii+$a]; P1[freedim[2]]=grid[jj+$b];
               P3[freedim[1]]=grid[ii+$c]; P3[freedim[2]]=grid[jj+$d];
               push!(mesh,Triangle(P1,P2,P3,area=area))
    end)
end
# generates a mesh of unit cube withe a fichera corner. N determines the size.
function gen_ficheramesh(N)

    if (mod(N, 2) != 0)   # Enforce N to be even
	N = N + 1;
    end

    nn = round(Int,N / 2); # Half the size

    # Prepare the mesh: Vector of Triangle
    mesh = Vector{Triangle}(undef,0);

    # compute the area of each triangle
    area = 0.25/N/N;

    # Uniform grid on [0,1]
    grid = Vector(0:N) / N;

    # Metadata for the discretization
    fixdim=[1;2;3;1;1;1;1;2;2;2;2;3;3;3;3]
    fixdim_val=[0;0;0;1;1;1;0.5;1;1;1;0.5;1;1;1;0.5]
    freedims=[2 3; 3 1; 1 2; 2 3; 2 3; 2 3 ; 2 3;3 1;3 1; 3 1; 3 1; 1 2; 1 2; 1 2; 1 2];
    Nvals=[1 N 1 N;
           1 N 1 N;
           1 N 1 N;
           1 nn 1 nn;
           nn+1 N 1 nn;
           1 nn nn+1 N;
           nn+1 N nn+1 N
           1 nn 1 nn
           nn+1 N 1 nn
           1 nn nn+1 N
           nn+1 N nn+1 N
           1 nn 1 nn
           nn+1 N 1 nn
           1 nn nn+1 N
           nn+1 N nn+1 N]
    # Ordering:
    # left face
    # front
    # bottom
    # right
    # (1/2, 1) x (0, 1/2)
    # (0, 1/2) x (1/2, 1)
    # (1/2, 1) x (1/2, 1)
    # back face: (0, 1/2) x (0, 1/2)
    # (1/2, 1) x (0, 1/2)
    # (0, 1/2) x (1/2, 1)
    # (1/2, 1) x (1/2, 1)
    # top face: (0, 1/2) x (0, 1/2)
    # (1/2, 1) x (0, 1/2)
    # (0, 1/2) x (1/2, 1)
    # (1/2, 1) x (1/2, 1)

    for l=1:size(fixdim,1);
        for ii = Nvals[l,1]:Nvals[l,2]
	    for jj = Nvals[l,3]:Nvals[l,4]

                center = zeros(3);
                center[fixdim[l]]=fixdim_val[l];         # Fix dimension
                freedim=freedims[l,:];
                # freedim now contains freedim[1] and freedim[2] which
                # correspond to the non-fixed dimensions
                center[freedim[1]]=(grid[ii]+grid[ii+1])/2;
                center[freedim[2]]=(grid[jj]+grid[jj+1])/2;

                # Now create the four triangles associated with this "center"
                if (l<4)
                    @addtri(0,0,1,0);
                    @addtri(1,0,1,1);
                    @addtri(1,1,0,1);
                    @addtri(0,1,0,0);
                else # For constency with Effenberger
                    @addtri(0,0,0,1);
                    @addtri(0,1,1,1);
                    @addtri(1,1,1,0);
                    @addtri(1,0,0,0);
	        end
            end
        end
    end

    # Auxiliary mesh data
    for tri = mesh
        # determine midpoint of triangle
        tri.aux.midpoint = (tri.P1 + tri.P2 + tri.P3) / 3;

        # normalized tangent vectors
        tri.aux.tau1 = normalize(tri.P2 - tri.P1);
        tri.aux.tau2 = normalize(tri.P3 - tri.P2);
        tri.aux.tau3 = normalize(tri.P1 - tri.P3);

        # exterior (normalized) normal of triangle
        tri.aux.normal = normalize(cross( tri.aux.tau1, tri.aux.tau2 ));

        # exterior (normalized) normals of the edges
        tri.aux.nu1 = normalize(cross(tri.aux.tau1, tri.aux.normal));
        tri.aux.nu2 = normalize(cross(tri.aux.tau2, tri.aux.normal));
        tri.aux.nu3 = normalize(cross(tri.aux.tau3, tri.aux.normal));
    end
    return mesh
end


text:
test_LP() = CoreModel(
    zeros(4, 3),
    zeros(4),
    ones(3),
    ones(3),
    ones(3),
    ["r$x" for x = 1:3],
    ["m$x" for x = 1:4],
)

test_simpleLP() = CoreModel(
    [
        1.0 1.0
        -1.0 1.0
    ],
    [3.0, 1.0],
    [-0.25, 1.0],
    -ones(2),
    2.0 * ones(2),
    ["r$x" for x = 1:2],
    ["m$x" for x = 1:2],
)

test_simpleLP2() = CoreModel(
    zeros(2, 2),
    [0.0, 0.0],
    [-0.25, 1.0],
    -ones(2),
    2.0 * ones(2),
    ["r$x" for x = 1:2],
    ["m$x" for x = 1:2],
)

test_sparseLP() = CoreModel(
    sprand(4000, 3000, 0.5),
    sprand(4000, 0.5),
    sprand(3000, 0.5),
    sprand(3000, 0.5),
    sprand(3000, 0.5),
    ["r$x" for x = 1:3000],
    ["m$x" for x = 1:4000],
)

test_coupledLP() = CoreModelCoupled(
    CoreModel(
        sprand(4000, 3000, 0.5),
        sprand(4000, 0.5),
        sprand(3000, 0.5),
        sprand(3000, 0.5),
        sprand(3000, 0.5),
        ["r$x" for x = 1:3000],
        ["m$x" for x = 1:4000],
    ),
    sprand(2000, 3000, 0.5),
    sprand(2000, 0.5),
    sprand(2000, 0.5),
)

test_toyModel() = CoreModel(
    [
        -1.0 1.0 0.0 0.0 0.0 0.0 0.0
        -2.0 0.0 1.0 0.0 0.0 0.0 0.0
        1.0 0.0 0.0 0.0 0.0 0.0 -1.0
        0.0 -1.0 0.0 -1.0 0.0 0.0 0.0
        0.0 0.0 -1.0 0.0 -1.0 0.0 0.0
        0.0 0.0 0.0 0.0 0.0 -1.0 1.0
    ],
    zeros(6),
    [0, 0, 0, 0, 0, 0, 1.0],
    fill(-1000.0, 7),
    fill(1000.0, 7),
    ["r1", "m1t", "m3t", "EX_m1(e)", "EX_m3(e)", "EX_biomass(e)", "biomass1"],
    ["m1[c]", "m3[c]", "m2[c]", "m1[e]", "m3[e]", "biomass[c]"],
)

const reaction_standard_gibbs_free_energies = Dict(
    #=
    ΔᵣG⁰ data from Equilibrator using the E. coli core model's reactions
    To generate this data manually, go to https://equilibrator.weizmann.ac.il/ and
    enter each reaction into the search bar, fill in the ΔG⁰ below from there. To generate
    automatically, use the eQuilibrator.jl package.
    =#
    "ACALD" => -21.268198981756314,
    "PTAr" => 8.651025243027988,
    "ALCD2x" => 17.47260052762408,
    "PDH" => -34.246780702043225,
    "PYK" => -24.48733600711958,
    "CO2t" => 0.0,
    "MALt2_2" => -6.839209921974724,
    "CS" => -39.330940148207475,
    "PGM" => -4.470553692565886,
    "TKT1" => -1.4976299544215408,
    "ACONTa" => 8.46962176350985,
    "GLNS" => -15.771242654033706,
    "ICL" => 9.53738025507404,
    "FBA" => 23.376920310319235,
    "SUCCt3" => -43.97082007726285,
    "FORt2" => -3.4211767267069124,
    "G6PDH2r" => -7.3971867270165035,
    "AKGDH" => -28.235442362320782,
    "TKT2" => -10.318436817276165,
    "FRD7" => 73.61589524884772,
    "SUCOAS" => -1.1586968559555615,
    "FBP" => -11.606887851480572,
    "ICDHyr" => 5.398885342794983,
    "AKGt2r" => 10.085299238039797,
    "GLUSy" => -47.21690395283849,
    "TPI" => 5.621932460512994,
    "FORt" => 13.509690780237293,
    "ACONTb" => -1.622946931741609,
    "GLNabc" => -30.194559792842753,
    "RPE" => -3.3881719029747615,
    "ACKr" => 14.027197131450492,
    "THD2" => -33.84686533309243,
    "PFL" => -19.814421615735,
    "RPI" => 4.477649590945703,
    "D_LACt2" => -3.4223903975852004,
    "TALA" => -0.949571985515206,
    "PPCK" => 10.659841402564751,
    "ACt2r" => -3.4196348363397995,
    "NH4t" => -13.606633927039097,
    "PGL" => -25.94931748161696,
    "NADTRHD" => -0.014869680795754903,
    "PGK" => 19.57192102020454,
    "LDH_D" => 20.04059765689044,
    "ME1" => 12.084968268076864,
    "PIt2r" => 10.415108493818785,
    "ATPS4r" => -37.570267233299816,
    "PYRt2" => -3.422891289768689,
    "GLCpts" => -45.42430981510088,
    "GLUDy" => 32.834943812395665,
    "CYTBD" => -59.700410815493775,
    "FUMt2_2" => -6.845105500839001,
    "FRUpts2" => -42.67529760694199,
    "GAPD" => 0.5307809794271634,
    "H2Ot" => -1.5987211554602254e-14,
    "PPC" => -40.81304419704113,
    "NADH16" => -80.37770501380615,
    "PFK" => -18.546314942995934,
    "MDH" => 25.912462872631522,
    "PGI" => 2.6307087407442395,
    "O2t" => 0.0,
    "ME2" => 12.099837948872533,
    "GND" => 10.312275879236381,
    "SUCCt2_2" => -6.82178244356977,
    "GLUN" => -14.381960140443113,
    "ETOHt2r" => -16.930867506944217,
    "ADK1" => 0.3893321583896068,
    "ACALDt" => -3.197442310920451e-14,
    "SUCDi" => -73.61589524884772,
    "ENO" => -3.8108376097261782,
    "MALS" => -39.22150045995042,
    "GLUt2r" => -3.499043558772904,
    "PPS" => -6.0551989457468665,
    "FUM" => -3.424133018702122,
)


text:
const rex_resrange = r"/(\d+)-(\d+)$"
function residue_range(name::AbstractString)
    m = match(rex_resrange, name)
    m === nothing && error("no residue range found")
    return parse(Int, m.captures[1]):parse(Int, m.captures[2])
end
function strip_residue_range(name::AbstractString)
    m = match(rex_resrange, name)
    m === nothing && error("no residue range found")
    return name[begin:prevind(name, m.offset)]
end

const rex_uniprot_species = r"_([A-Z0-9]{1,5})($|/)"

"""
    species(name)

Extract the species identifier from a [UniProt "X_Y" entry](https://www.uniprot.org/help/entry_name)
or elaborated variant (e.g., PFAM sequence name).

See also [`uniprotX`](@ref).

# Examples

```jldoctest
julia> species("Q8VGW6_MOUSE/31-308")
"MOUSE"
```
"""
function species(name::AbstractString)
    m = match(rex_uniprot_species, name)
    m === nothing && error("no species found")
    return m.captures[1]::AbstractString
end

const rex_uniprotX_Swiss = r"^([A-Z0-9]{1,5})($|_)"   # Swiss
const rex_uniprot_accession = r"^([OPQ][0-9][A-Z0-9]{3}[0-9]|[A-NR-Z][0-9]([A-Z][A-Z0-9]{2}[0-9]){1,2})($|_)"

"""
    uniprotX(name)

Extract the [UniProt "X" entry](https://www.uniprot.org/help/entry_name) from an X_Y entry or elaborated variant (e.g., PFAM sequence name).

See also [`species`](@ref).

# Examples

```jldoctest
julia> uniprotX("Q8VGW6_MOUSE/31-308")
"Q8VGW6"
```
"""
function uniprotX(name::AbstractString)::AbstractString
    m = match(rex_uniprotX_Swiss, name)
    m !== nothing && return m.captures[1]
    m = match(rex_uniprot_accession, name)
    m !== nothing && return m.captures[1]
    error("no UniProt X identifier found")
end


text:
"""
This test is the using the same procedure as identificationRobustInference.jmd
with using new funciton which have been added to scr/HypothesisTests.
For compating, the same methods in identificationRobustInference.jmd also
has been presented.
"""

using Optim, ForwardDiff, LinearAlgebra, Distributions, LaTeXStrings
"""
    AR_gmm_obj(gi:Function)

Returns a function `ar(θ)` where:
``ar(θ) = n [1/n \\sum g_i(θ)]' \\widehat{Var}(g_i(θ))^{-1}[1/n \\sum g_i(θ)] ``
"""
function AR_gmm_obj(gi::Function)
  function(θ)
    g = gi(θ)
    m = mean(g,dims=1)
    W = inv(cov(g))
    (size(g)[1]*( m*W*m')[1])
  end
end
"""
    AR_test(n::Integer, β::AbstractVector,
                      π::AbstractMatrix, ρ)

Returns a function `ar(θ)` where:
``ar(θ) = n [1/n \\sum g_i(θ)]' \\widehat{Var}(g_i(θ))^{-1}[1/n \\sum g_i(θ)] ``

Or equivalent to:
``AR_gmm_obj(get_gi(IVLogitShare(n,β0,π0,ρ)))``
"""
function AR_test(n::Integer, β::AbstractVector,
                      π::AbstractMatrix, ρ)
    ivlogitshare = IVLogitShare(n,β,π,ρ)
    gi = get_gi(ivlogitshare)
    return AR_gmm_obj(gi)
end

"""
There are also identification robust versions of likelihood ratio and
lagrange multiplier test. Moreire (2003)[@moreira2003] proposed a conditional
likelihood ratio (CLR) test for weakly identified linear IV
models. Kleibergen (2005)[@kleibergen2005] developed a Lagrange multiplier (often called
the KLM) test and extended Moreira's CLR test to weakly identified GMM
models.  More recently, Andrews and Guggenberge (2015) [@andrews2015]
and Andrews and Guggenberge (2017) [@andrews2017] showed the
validity of these tests under more general conditions. These tests are
somewhat more complicated than the AR test, but they have the
advantage that they are often more powerful. The AR test statistic has
a $\chi^2_{m}$ distribution, where $m$ is the number of moments. The
CLR and KLM statistics under strong identification have $\chi^2_k$
distributions (as does the Wald statistic), where $k$ is the number of
parameters. Consequently, when the model is overidentified, the CLR
and LM tests are more powerful than the AR test.

Here is an implementation of the KLM and CLR statistics. The names of
variables roughly follows the notation of Andrews and
Guggenberge(2017

"""

"""
    statparts(gi::Function)

Return a function.

compute common components of klm, rk, & clr stats
follows notation of Andrews & Guggenberger 2017, section 3.1
"""
function statparts(gi::Function)
  # compute common components of klm, rk, & clr stats
  # follows notation of Andrews & Guggenberger 2017, section 3.1
    function P(A::AbstractMatrix) # projection matrix
    A*pinv(A'*A)*A'
    end
    function(θ)
        giθ = gi(θ)
        p = length(θ)
        (n, k) = size(giθ)
        Ω = cov(giθ)
        gn=mean(gi(θ), dims=1)'
        #G = ForwardDiff.jacobian(θ->mean(gi(θ),dims=1),θ)
        Gi= ForwardDiff.jacobian(gi,θ)
        Gi = reshape(Gi, n , k, p)
        G = mean(Gi, dims=1)
        Γ = zeros(eltype(G),p,k,k)
        D = zeros(eltype(G),k, p)
        for j in 1:p
        for i in 1:n
          Γ[j,:,:] += (Gi[i,:,j] .- G[1,:,j]) * giθ[i,:]'
        end
        Γ[j,:,:] ./= n
        D[:,j] = G[1,:,j] - Γ[j,:,:]*inv(Ω)*gn
        end
        return(n,k,p,gn, Ω, D, P)
    end
end

"""
    klm(gi::Function)

Return test statistic for KLM method.
"""

function klm(gi::Function)
    SP = statparts(gi)
    function(θ)
        (n,k,p,gn, Ω, D, P) = SP(θ)
        return n*(gn'*Ω^(-1/2)*P(Ω^(-1/2)*D)*Ω^(-1/2)*gn)[1]
    end
end

"""
    clr(gi::Function)

Return test statistic for CLR method.
"""
function clr(gi::Function)
  SP = statparts(gi)
  function(θ)
        (n,k,p,gn, Ω, D, P) = SP(θ)
        rk = eigmin(n*D'*inv(Ω)*D)
        AR  = (n*gn'*inv(Ω)*gn)[1]
        lm = (n*(gn'*Ω^(-1/2)*P(Ω^(-1/2)*D)*Ω^(-1/2)*gn))[1]
        lr = 1/2*(AR - rk + sqrt( (AR-rk)^2 + 4*lm*rk))

        # simulate to find p-value
        S = 5000
        function randc(k,p,r,S)
        χp = rand(Chisq(p),S)
        χkp = rand(Chisq(k-p),S)
        0.5.*(χp .+ χkp .- r .+
              sqrt.((χp .+ χkp .- r).^2 .+ 4 .* χp.*r))
        end
        csim = randc(k,p,rk,S)
        return mean(csim.<=lr)
    end
end

function plot_cr(β,V, tests::AbstractArray{Function}, labels; ngrid=30)
  lb = β - sqrt.(diag(V))*5
  ub = β + sqrt.(diag(V))*5
  fig=scatter([β0[1]], [β0[2]], markersize=8, legend=false,
              xlabel=L"\beta_1", ylabel=L"\beta_2")
  ntest = 1000
  βtest = [rand(2).*(ub-lb) .+ lb for i in 1:ntest]
  pval = tests[1].(βtest)
  βtest = vcat(βtest'...)
  crit = 0.9
  fig=scatter!(βtest[:,1],βtest[:,2], group=(pval.<crit), legend=false,
               markersize=4, markerstrokewidth=0.0, seriesalpha=0.5,
               palette=:heat)
  b1 = lb[1]:(ub[1]-lb[1])/ngrid:ub[1]
  b2 = lb[2]:(ub[2]-lb[2])/ngrid:ub[2]
  colors = [:black, :red, :blue, :green]
  for t in 1:length(tests)
    fig=contour!(b1,b2,(a,b)->tests[t]([a,b]),
             levels = [0.9, 0.95],
             contour_labels=false, legend=false,
             label = labels[t],
             c=cgrad([colors[t],colors[t]],[0.0,1.0]))
  end
  fig
end

function gmmVar(θ,gi,W)
  g = gi(θ)
  n = size(g)[1]
  D = ForwardDiff.jacobian(θ->mean(gi(θ),dims=1),θ)
  Σ = cov(gi(θ))
  1/n*inv(D'*W*D)*(D'*W*Σ*W*D)*inv(D'*W*D)
end
function gmmObj(θ,gi,W)
  g = gi(θ)
  m = mean(g,dims=1)
  (size(g)[1]*( m*W*m')[1]) # return scalar, not 1x1 array
end

n = 50
k = 2
iv =3
π0 = vcat(0.1*diagm(0=>ones(k)),0.2*ones(iv-k,k))
ρ = 0.5
data = IVLogitShare(n,β0,π0,ρ)
opt1 = optimize(θ->gmmObj(θ, get_gi(data) ,I),
                β0, BFGS(), autodiff =:forward)
β1 = opt1.minimizer
V1 = gmmVar(β1,get_gi(data),I)

KLM_f = klm(get_gi(data))
pklm = θ -> cdf(Chisq(length(β1)), KLM_f(θ) )

AR_f = AR_gmm_obj(get_gi(data))
par  = θ ->cdf(Chisq(size(data.z)[2]), AR_f(θ) )
pclr  = clr(get_gi(data))
pwald = θ -> cdf(Chisq(length(β1)),(θ-β1)'*inv(V1)*(θ-β1))
plot_cr(β1,V1, [pclr, pklm, par, pwald],
        ["CLR","KLM","AR","Wald"], ngrid=40)

"""
Now the equivalent codes from part of identificationRobustInference.jmd
For comparison
"""
function simulate_ivshare(n,β,γ,ρ)
  z = randn(n, size(γ)[1])
  endo = randn(n, length(β))
  x = z*γ .+ endo
  ξ = rand(Normal(0,sqrt((1.0-ρ^2))),n).+endo[:,1]*ρ
  y = cdf.(Logistic(), x*β .+ ξ)
  return((y=y,x=x,z=z))
end
n = 100
k = 2
iv = 3
β0 = ones(k)
π0 = vcat(5*I,ones(iv-k,k))
ρ = 0.5
(y,x,z) = simulate_ivshare(n,β0,π0,ρ)

function gi_ivshare(β,y,x,z)
  ξ = quantile.(Logistic(),y) .- x*β
  ξ.*z
end

function gmmObj(θ,gi,W)
  g = gi(θ)
  m = mean(g,dims=1)
  (size(g)[1]*( m*W*m')[1]) # return scalar, not 1x1 array
end

function gmmVar(θ,gi,W)
  g = gi(θ)
  n = size(g)[1]
  D = ForwardDiff.jacobian(θ->mean(gi(θ),dims=1),θ)
  Σ = cov(gi(θ))
  1/n*inv(D'*W*D)*(D'*W*Σ*W*D)*inv(D'*W*D)
end

function ar(θ,gi)
  gmmObj(θ,gi,inv(cov(gi(θ))))
end

using ForwardDiff, Plots, Optim
Plots.gr()
function statparts(θ,gi)
  # compute common components of klm, rk, & clr stats
  # follows notation of Andrews & Guggenberger 2017, section 3.1
  function P(A::AbstractMatrix) # projection matrix
    A*pinv(A'*A)*A'
  end
  giθ = gi(θ)
  p = length(θ)
  (n, k) = size(giθ)
  Ω = cov(giθ)
  gn=mean(gi(θ), dims=1)'
  #G = ForwardDiff.jacobian(θ->mean(gi(θ),dims=1),θ)
  Gi= ForwardDiff.jacobian(gi,θ)
  Gi = reshape(Gi, n , k, p)
  G = mean(Gi, dims=1)
  Γ = zeros(eltype(G),p,k,k)
  D = zeros(eltype(G),k, p)
  for j in 1:p
    for i in 1:n
      Γ[j,:,:] += (Gi[i,:,j] .- G[1,:,j]) * giθ[i,:]'
    end
    Γ[j,:,:] ./= n
    D[:,j] = G[1,:,j] - Γ[j,:,:]*inv(Ω)*gn
  end
  return(n,k,p,gn, Ω, D, P)
end

function klm(θ,gi)
  (n,k,p,gn, Ω, D, P) = statparts(θ,gi)
  lm = n*(gn'*Ω^(-1/2)*P(Ω^(-1/2)*D)*Ω^(-1/2)*gn)[1]
end

function clr(θ,gi)
  (n,k,p,gn, Ω, D, P) = statparts(θ,gi)

  rk = eigmin(n*D'*inv(Ω)*D)
  AR  = (n*gn'*inv(Ω)*gn)[1]
  lm = (n*(gn'*Ω^(-1/2)*P(Ω^(-1/2)*D)*Ω^(-1/2)*gn))[1]
  lr = 1/2*(AR - rk + sqrt( (AR-rk)^2 + 4*lm*rk))

  # simulate to find p-value
  S = 5000
  function randc(k,p,r,S)
    χp = rand(Chisq(p),S)
    χkp = rand(Chisq(k-p),S)
    0.5.*(χp .+ χkp .- r .+
          sqrt.((χp .+ χkp .- r).^2 .+ 4 .* χp.*r))
  end
  csim = randc(k,p,rk,S)
  pval = mean(csim.<=lr)
end

function plot_cr(β,V, tests::AbstractArray{Function}, labels; ngrid=30)
  lb = β - sqrt.(diag(V))*5
  ub = β + sqrt.(diag(V))*5
  fig=scatter([β0[1]], [β0[2]], markersize=8, legend=false,
              xlabel=L"\beta_1", ylabel=L"\beta_2")
  ntest = 1000
  βtest = [rand(2).*(ub-lb) .+ lb for i in 1:ntest]
  pval = tests[1].(βtest)
  βtest = vcat(βtest'...)
  crit = 0.9
  fig=scatter!(βtest[:,1],βtest[:,2], group=(pval.<crit), legend=false,
               markersize=4, markerstrokewidth=0.0, seriesalpha=0.5,
               palette=:heat)
  b1 = lb[1]:(ub[1]-lb[1])/ngrid:ub[1]
  b2 = lb[2]:(ub[2]-lb[2])/ngrid:ub[2]
  colors = [:black, :red, :blue, :green]
  for t in 1:length(tests)
    fig=contour!(b1,b2,(a,b)->tests[t]([a,b]),
             levels = [0.9, 0.95],
             contour_labels=false, legend=false,
             label = labels[t],
             c=cgrad([colors[t],colors[t]],[0.0,1.0]))
  end
  fig
end

n = 50
k = 2
iv =3
π0 = vcat(0.1*diagm(0=>ones(k)),0.2*ones(iv-k,k))
ρ = 0.5
(y,x,z) = simulate_ivshare(n,β0,π0,ρ)
opt1 = optimize(θ->gmmObj(θ, β->gi_ivshare(β,y,x,z) ,I),
                β0, BFGS(), autodiff =:forward)
β1 = opt1.minimizer
V1 = gmmVar(β1,β->gi_ivshare(β,y,x,z),I)

pklm = θ->cdf(Chisq(length(βcue)),klm(θ, β->gi_ivshare(β,y,x,z)))
par  = θ->cdf(Chisq(size(z)[2]), ar(θ, β->gi_ivshare(β,y,x,z)))
pclr  = θ->clr(θ, β->gi_ivshare(β,y,x,z))
pwald = θ -> cdf(Chisq(length(β1)),(θ-β1)'*inv(V1)*(θ-β1))
plot_cr(β1,V1, [pclr, pklm, par, pwald],
        ["CLR","KLM","AR","Wald"], ngrid=40)


text:
"""
An ellipsoidal representation of the Earth, for converting between LLA and
other co-ordinate systems such as ECEF.
"""
struct Ellipsoid
    a::Float64        # Semi-major axis
    b::Float64        # Semi-minor axis
    f::Float64        # Flattening
    e2::Float64       # Eccentricity squared
    name::Symbol      # Conventional name - for clarity, should match the name
                      # of the const instance in the package!
end

function Ellipsoid(; a::String="", b::String="", f_inv::String="", name=:UNKNOWN)
    if isempty(a) || isempty(b) == isempty(f_inv)
        throw(ArgumentError("Specify parameter 'a' and either 'b' or 'f_inv'"))
    end
    if isempty(b)
        _ellipsoid_af(parse(BigFloat, a), parse(BigFloat, f_inv), name)
    else
        _ellipsoid_ab(parse(BigFloat, a), parse(BigFloat, b), name)
    end
end

function _ellipsoid_ab(a::BigFloat, b::BigFloat, name)
    f = 1 - b/a
    e2 = f*(2-f)
    Ellipsoid(a, b, f, e2, name)
end
function _ellipsoid_af(a::BigFloat, f_inv::BigFloat, name)
    b = a * (1 - inv(f_inv))

    _ellipsoid_ab(a, b, name)
end

function Base.show(io::IO, el::Ellipsoid)
    if el.name != :UNKNOWN
        # To clarify that these are Ellipsoids, we wrap the name in
        # 'Ellipsoid', even though the name itself should resolve to the
        # correct ellipsoid instance.
        print(io, "Ellipsoid($(el.name))")
    else
        print(io, "Ellipsoid(a=$(el.a), b=$(el.b))")
    end
end

#-------------------------------------------------------------------------------
# Standard ellipsoids

# Worldwide
const wgs84_ellipsoid   = Ellipsoid(a = "6378137.0", f_inv = "298.257223563", name=:wgs84_ellipsoid)
const grs80   = Ellipsoid(a = "6378137.0", f_inv = "298.2572221008827112431628366", name=:grs80) #NB: not the definition - f_inv is derived.

# Regional
const airy1830   = Ellipsoid(a = "6377563.396", b = "6356256.909", name=:airy1830)
const clarke1866 = Ellipsoid(a = "6378206.4",   b = "6356583.8", name=:clarke1866)

ellipsoid(x) = error("No ellipsoid defined for datum $x")


text:
using Genie.Router

route("/") do
  serve_static_file("/welcome.html")
end


text:
"""
    initialize_dual_potential!(::AbstractDivergence, a::DiscreteMeasure) -> Nothing

Apply an initialization for the dual potential, for use in
[`unbalanced_sinkhorn!`](@ref); falls back to zeroing out the dual potential.
Specialized implementations can improve performance, but should not affect
correctness.
"""
function initialize_dual_potential!(::AbstractDivergence, a::DiscreteMeasure)
    a.dual_potential .= 0
end

"""
    function unbalanced_sinkhorn!(
        D::AbstractDivergence,
        C,
        a::DiscreteMeasure,
        b::DiscreteMeasure,
        ϵ = 1e-1;
        tol = 1e-5,
        max_iters = 10^5,
        warn::Bool = true,
    ) -> NamedTuple

Implements algorithm 1 of [[SFVTP19](@ref)]. The `dual_potential` fields of `a`
and `b` are updated to hold the optimal dual potentials. The `density`,
`log_density`, and `set` fields are not modified. The parameters are

* `D`: the [`AbstractDivergence`](@ref) used for measuring the cost of
  creating or destroying mass
* `ϵ`: the regularization parameter
* `C`: either a function from `a.set` × `b.set` to real numbers; should satisfy
  `C(x,y) = C(y,x)` and `C(x,x)=0` when applicable, or a precomputed cost matrix, see [`cost_matrix`](@ref)
* `tol`: the convergence tolerance
* `max_iters`: the maximum number of iterations to perform.
* `warn`: whether or not to warn when the maximum number of iterations is reached.

Returns a NamedTuple of the number of iterations performed (`iters`), and the
maximum residual (`max_residual`), which is the maximum infinity norm difference
between consecutive iterates of the dual potentials, at the end of the process.
If `max_iters` is not reached, iteration stops when the `max_residual` falls
below `tol`.

"""
function unbalanced_sinkhorn!(
    D::AbstractDivergence,
    C::AbstractMatrix,
    a::DiscreteMeasure,
    b::DiscreteMeasure,
    ϵ = 1e-1;
    tol = 1e-5,
    max_iters = 10^5,
    warn::Bool = true,
)
    if D isa Balanced && warn && sum(a.density) ≉ sum(b.density)
        @warn "Should have `sum(a.density) ≈ sum(b.density)` for `D==Balanced()`"
    end

    initialize_dual_potential!(D, a)
    initialize_dual_potential!(D, b)

    f = a.dual_potential
    tmp_f = a.cache

    g = b.dual_potential
    tmp_g = b.cache

    max_residual = Inf
    iters = 0
    while iters < max_iters && max_residual > tol
        iters += 1
        max_residual = 0.0
        @inbounds for j in eachindex(g)
            for i in eachindex(a.log_density, f, tmp_f)
                tmp_f[i] = a.log_density[i] + (f[i] - C[i, j]) / ϵ
            end
            new_g = -ϵ * logsumexp!(tmp_f)
            new_g = -aprox(D, ϵ, -new_g)
            diff = abs(g[j] - new_g)
            if diff > max_residual
                max_residual = diff
            end
            g[j] = new_g
        end
        @inbounds for j in eachindex(f)
            for i in eachindex(b.log_density, g, tmp_g)
                tmp_g[i] = b.log_density[i] + (g[i] - C[j, i]) / ϵ
            end
            new_f = -ϵ * logsumexp!(tmp_g)
            new_f = -aprox(D, ϵ, -new_f)
            diff = abs(f[j] - new_f)
            if diff > max_residual
                max_residual = diff
            end
            f[j] = new_f
        end
    end

    if warn && iters == max_iters
        @warn "Maximum iterations ($max_iters) reached" max_residual
    end

    return (iters = iters, max_residual = max_residual)
end

"""
    function OT!(
        D::AbstractDivergence,
        C,
        a::DiscreteMeasure,
        b::DiscreteMeasure,
        ϵ = 1e-1;
        C = (x, y) -> norm(x - y),
        kwargs...,
    ) -> Number

Computes the optimal transport cost between `a` and `b`, using
[`unbalanced_sinkhorn!`](@ref); see that function for the meaning of the
parameters and the keyword arguments. Implements Equation (15) of
[[SFVTP19](@ref)].
"""
function OT!(
    D::AbstractDivergence,
    C::AbstractMatrix,
    a::DiscreteMeasure,
    b::DiscreteMeasure,
    ϵ = 1e-1;
    kwargs...,
)

    unbalanced_sinkhorn!(D, C, a, b, ϵ; kwargs...)

    f = a.dual_potential
    g = b.dual_potential

    T = promote_type(eltype(a), eltype(b))
    _nφstar = q -> -φstar(D, -q)

    t1 = fdot(_nφstar, a.density, f)
    t2 = fdot(_nφstar, b.density, g)
    t3 = zero(T)
    for i in eachindex(f), j in eachindex(g)
        t3 -=
            ϵ *
            a.density[i] *
            b.density[j] *
            (exp((f[i] + g[j] - C[i, j]) / ϵ) - one(T))
    end
    return t1 + t2 + t3
end

# Generic method
function _sinkhorn_divergence!(
    D::AbstractDivergence,
    C,
    a::DiscreteMeasure,
    b::DiscreteMeasure,
    ϵ;
    kwargs...,
)
    OT_aa = OT!(D, C, a, a, ϵ; kwargs...)
    OT_bb = OT!(D, C, b, b, ϵ; kwargs...)
    OT_ab = OT!(D, C, a, b, ϵ; kwargs...)

    OT_ab + (-OT_aa - OT_bb + ϵ * (sum(a.density) - sum(b.density))^2) / 2
end

"""
    sinkhorn_divergence!(
        D::AbstractDivergence,
        C,
        a::DiscreteMeasure,
        b::DiscreteMeasure,
        ϵ = 1e-1;
        kwargs...,
    ) -> Number

Computes the unbalanced sinkhorn divergence between `a` and `b` as defined in
Def. 6 of [[SFVTP19](@ref)], using [`unbalanced_sinkhorn!`](@ref); see that
function for the meaning of the parameters and the keyword arguments. Sets the
optimal `dual_potential`'s of `a` and `b`.
"""
sinkhorn_divergence!(
    D::AbstractDivergence,
    C,
    a::DiscreteMeasure,
    b::DiscreteMeasure,
    ϵ = 1e-1;
    kwargs...,
) = _sinkhorn_divergence!(D, C, a, b, ϵ; kwargs...)

sinkhorn_divergence!(
    D::AbstractDivergence,
    C::AbstractMatrix,
    a::DiscreteMeasure,
    b::DiscreteMeasure,
    ϵ = 1e-1;
    kwargs...,
) = throw(ArgumentError("Must pass a cost function `C`, not a cost matrix."))

"""
    function optimal_coupling!(
        D::AbstractDivergence,
        C,
        a::DiscreteMeasure,
        b::DiscreteMeasure,
        ϵ = 1e-1;
        dual_potentials_populated::Bool = false,
        kwargs...) -> Matrix

Computes the optimal coupling between `a` and `b` using the optimal dual
potentials, the regularization parameter `ϵ`, and the cost function `C`.

If `dual_potentials_populated = false`, [`unbalanced_sinkhorn!`](@ref) is
called to populate the dual potentials of `a` and `b`, using the divergence `D`.
If `dual_potentials_populated = true`, one of [`unbalanced_sinkhorn!`](@ref)
or [`OT!`](@ref) or [`sinkhorn_divergence!`](@ref) must be called first to
set the optimal dual potentials, with the same choice of `ϵ` and `C`. In this case,
`a` and `b` are not mutated.

This function implements Prop.
6 of [[SFVTP19](@ref)].
"""
function optimal_coupling!(
    D::AbstractDivergence,
    C::AbstractMatrix,
    a::DiscreteMeasure,
    b::DiscreteMeasure,
    ϵ = 1e-1;
    dual_potentials_populated::Bool = false,
    kwargs...,
)
    if !dual_potentials_populated
        unbalanced_sinkhorn!(D, C, a, b, ϵ; kwargs...)
    end

    f = a.dual_potential
    g = b.dual_potential
    return [
        exp((f[i] + g[j] - C[i, j]) / ϵ) * a.density[i] * b.density[j]
        for i in eachindex(a.density), j in eachindex(b.density)
    ]
end

# Provide default choice of norm as cost function
for fun in [:unbalanced_sinkhorn!, :OT!, :sinkhorn_divergence!, :optimal_coupling!]
    @eval begin
        $fun(
            D::AbstractDivergence,
            a::DiscreteMeasure,
            b::DiscreteMeasure,
            args...;
            kwargs...
        ) = $fun(D, (x,y)->norm(x-y), a, b, args...; kwargs...)
    end
end

# If a function is provided, precompute the cost matrix. sinkhorn_divergence is omitted since it requires three different cost matrices.
for fun in [:unbalanced_sinkhorn!, :OT!, :optimal_coupling!]
    @eval begin
        $fun(
            D::AbstractDivergence,
            C, # Not restricting to Function to allow callable structs etc.
            a::DiscreteMeasure,
            b::DiscreteMeasure,
            args...;
            kwargs...
        ) = $fun(D, cost_matrix(C, a, b), a, b, args...; kwargs...)
    end
end


text:
__precompile__()

module LogProbs

import Base: show, log, float, rand, isapprox
import Base: ==, hash
import Base: zero, one, +, -, *, /, isless

using StatsFuns: logaddexp, log1mexp

export LogProb, information

"""
    LogProb(x)

Create a container of a positive real number `x`
for convenient calculations with logspace probabilities

see https://en.wikipedia.org/wiki/Log_probability
"""
struct LogProb <: Number
    log :: Float64
    LogProb(number; islog::Bool=false) = islog ? new(number) : new(log(number))
end

show(io::IO, x::LogProb)  = print(io, "LogProb($(float(x)))")

==(x::LogProb,y::LogProb) = x.log == y.log
hash(x::LogProb)          = hash(LogProb, hash(x.log))

rand(::Type{LogProb})     = LogProb(rand())

float(x::LogProb)         = exp(x.log)
log(x::LogProb)           = x.log

"""
    information(p::LogProb)

Calculate the Shannon information content of `p` in bits
"""
information(x::LogProb)   = - log(x) / log(2)

one( ::Type{LogProb})     = LogProb(1)
zero(::Type{LogProb})     = LogProb(0)

*(x::LogProb, y::LogProb) = LogProb(x.log + y.log, islog=true)
/(x::LogProb, y::LogProb) = LogProb(x.log - y.log, islog=true)
+(x::LogProb, y::LogProb) = LogProb(logaddexp(x.log, y.log), islog=true)
-(x::LogProb, y::LogProb) = LogProb(x.log + log1mexp(y.log-x.log), islog=true)

isless(  x::LogProb, y::LogProb) = isless(  x.log, y.log)
isapprox(x::LogProb, y::LogProb) = isapprox(x.log, y.log)

end # Module