text:
include("../src/tenedyn_module.jl")
include("../../__lib__/math/ode/src/ode_module.jl")

module ttest
using ..tenedyn
using ..datastructs
using ..ode
using PyPlot
PyPlot.pygui(true)

struct OdeParams_PiPt1

    setpoint_amplitude::Float64
    setpoint_time_delay::Float64
    ctrl_kp::Float64
    ctrl_ki::Float64
    ctrl_integral_init_value::Float64
    plant_T::Float64
    plant_integral_init_value::Float64

    function OdeParams_PiPt1()
        setpoint_amplitude  = 2.0
        setpoint_time_delay = 0.5

        ctrl_kp = 10.0
        ctrl_ki = 5.0
        ctrl_integral_init_value = 0.0

        plant_T = 0.2
        plant_integral_init_value = 0.5

        return new(
            setpoint_amplitude,
            setpoint_time_delay,
            ctrl_kp,
            ctrl_ki,
            ctrl_integral_init_value,
            plant_T,
            plant_integral_init_value)
    end
end

function step(time, a, td)
    return time >= td ? a : zero(typeof(time))
end

function manual_ode_fcn(der_q, time, q, params::OdeParams_PiPt1)

    u2 = q[1]
    x  = q[2]

    e  = step(time, params.setpoint_amplitude, params.setpoint_time_delay) - x
    u1 = params.ctrl_kp * e
    u  = u1 + u2

    der_q[1] = params.ctrl_ki * e
    der_q[2] = (u - x) / params.plant_T

    return der_q
end

function create_model_pi_controller_pt1(
    setpoint_amplitude,
    setpoint_time_delay,
    ctrl_kp,
    ctrl_ki,
    ctrl_integral_init_value,
    plant_T,
    plant_integral_init_value)

    # model definition
    setpoint = tenedyn.Step(setpoint_amplitude, setpoint_time_delay)

    ctrl_minus      = tenedyn.Minus()
    ctrl_prop_gain  = tenedyn.Gain(ctrl_kp)
    ctrl_int_gain   = tenedyn.Gain(ctrl_ki)
    ctrl_plus       = tenedyn.Plus()
    ctrl_integrator = tenedyn.Integrator(ctrl_integral_init_value)

    plant_pt1_minus      = tenedyn.Minus()
    plant_pt1_mult1T     = tenedyn.Gain(1/plant_T)
    plant_pt1_integrator = tenedyn.Integrator(plant_integral_init_value)

    sensor = tenedyn.Identity()

    model = tenedyn.Model()

    # add elements to the model
    push!(model, setpoint,              "ref")                  # 1

    push!(model, ctrl_minus,            "ctrl_minus")           # 2
    push!(model, ctrl_prop_gain,        "ctrl_prop_gain")       # 3
    push!(model, ctrl_int_gain,         "ctrl_int_gain")        # 4
    push!(model, ctrl_plus,             "ctrl_plus")            # 5
    push!(model, ctrl_integrator,       "ctrl_integrator")      # 6

    push!(model, plant_pt1_minus,       "plant_pt1_minus")      # 7
    push!(model, plant_pt1_mult1T,      "plant_pt1_mult1T")     # 8
    push!(model, plant_pt1_integrator,  "plant_pt1_integrator") # 9

    push!(model, sensor,                "sensor")               # 10

    # topology
    tenedyn.connect!(model, setpoint.y,             ctrl_minus.x1)

    tenedyn.connect!(model, ctrl_minus.y,           ctrl_prop_gain.x)
    tenedyn.connect!(model, ctrl_prop_gain.y,       ctrl_plus.x1)

    tenedyn.connect!(model, ctrl_minus.y,           ctrl_int_gain.x)
    tenedyn.connect!(model, ctrl_int_gain.y,        ctrl_integrator.x)
    tenedyn.connect!(model, ctrl_integrator.y,      ctrl_plus.x2)

    tenedyn.connect!(model, ctrl_plus.y,            plant_pt1_minus.x1)
    tenedyn.connect!(model, plant_pt1_integrator.y, plant_pt1_minus.x2)

    tenedyn.connect!(model, plant_pt1_minus.y,      plant_pt1_mult1T.x)
    tenedyn.connect!(model, plant_pt1_mult1T.y,     plant_pt1_integrator.x)

    tenedyn.connect!(model, plant_pt1_integrator.y, sensor.x)
    tenedyn.connect!(model, sensor.y,               ctrl_minus.x2)

    return model
end

function create_parametrized_model()

    # model parameters
    setpoint_amplitude  = 2.0
    setpoint_time_delay = 0.5

    ctrl_kp = 10.0
    ctrl_ki = 5.0
    ctrl_integral_init_value = 0.0

    plant_T = 0.2
    plant_integral_init_value = 0.5

    @time model = create_model_pi_controller_pt1(
        setpoint_amplitude,
        setpoint_time_delay,
        ctrl_kp,
        ctrl_ki,
        ctrl_integral_init_value,
        plant_T,
        plant_integral_init_value)

    return model
end

function create_odefcn_pi_controller_pt1()

    # model parameters
    setpoint_amplitude  = 2.0
    setpoint_time_delay = 0.5

    ctrl_kp = 10.0
    ctrl_ki = 5.0
    ctrl_integral_init_value = 0.0

    plant_T = 0.2
    plant_integral_init_value = 0.5

    @time model = create_model_pi_controller_pt1(
        setpoint_amplitude,
        setpoint_time_delay,
        ctrl_kp,
        ctrl_ki,
        ctrl_integral_init_value,
        plant_T,
        plant_integral_init_value)

    @time (forest, integrator_indices, tree_node_vecvec, evalulation_indices_vecvec) = tenedyn.compile(model)

    ode_fcn = tenedyn.create_ode_function(
        model,
        forest,
        integrator_indices,
        tree_node_vecvec,
        evalulation_indices_vecvec)

    return ode_fcn
end

function run_simulation(ode_fcn)

    time_interval = ode.TimeInterval(50.0)
    init_cond     = [0.0, 0.5]

    ode_problem = ode.OdeProblem(
        ode_fcn,
        time_interval,
        init_cond)

    ode_solver = ode.Tsitouras45()

    basic_options = ode.BasicOdeOptions(
        ode_problem;
        abs_tol = 1e-5,
        rel_tol = 1e-5,
        max_iter = 1000
    )

    step_size_controller = ode.ClassicalStepSizeControl()

    stepper_options = ode.StepperOptions(;
        max_step_size = 5.0)

    @time (ode_res, step_stat) = ode.solve(
        ode_problem,
        ode_solver;
        basic_options        = basic_options,
        step_size_controller = step_size_controller,
        stepper_options      = stepper_options
        ) #

    # t = ode.time(ode_res)
    # q = ode.variables(ode_res)

    return ode_res
end

function plot_ode_res(ode_res)

    t = ode.time(ode_res)
    q = ode.variables(ode_res)

    figh = PyPlot.figure()
    a1 = PyPlot.subplot(2,1,1)

    PyPlot.grid()
    PyPlot.plot( t, q[1, :] )

    PyPlot.subplot(2,1,2, sharex = a1)
    PyPlot.grid()
    PyPlot.plot( t, q[2, :] )

    return figh
end

# after creation of a function string from the tree, only this is needed for more speed:
# fcn_str         = "fcn(t) = sin(t) + 3"
# parsed_expr     = Meta.parse(fcn_str)
# callable_object = eval(parsed_expr)

#=
ode_fcn_str_prefix  = "ode_fcn(der_q, time, q) = "
str_begin           = "begin "
str_end             = " end"
body                = "der_q[1] = - 0.2 * q[1] + sin(time) ;"

ode_fcn_str = ode_fcn_str_prefix * str_begin * body * str_end
=#

function create_callable(ode_fcn_str)

    parsed_expr = Meta.parse(ode_fcn_str)
    ode_fcn     = eval(parsed_expr)

    return ode_fcn
end

function sim2(ode_fcn)

    time_interval = ode.TimeInterval(5.0)
    init_cond     = [0.0, 0.5]

    ode_problem = ode.OdeProblem(
        ode_fcn,
        time_interval,
        init_cond)

    ode_solver = ode.Tsitouras45()

    basic_options = ode.BasicOdeOptions(
        ode_problem;
        abs_tol = 1e-5,
        rel_tol = 1e-5,
        max_iter = 1000
    )

    step_size_controller = ode.ClassicalStepSizeControl()

    stepper_options = ode.StepperOptions(;
        max_step_size = 0.1)

    @time (ode_res, step_stat) = ode.solve(
        ode_problem,
        ode_solver;
        basic_options        = basic_options,
        step_size_controller = step_size_controller,
        stepper_options      = stepper_options
        ) #

    return ode_res
end

function plot_sim2(ode_res)

    t = ode.time(ode_res)
    q = ode.variables(ode_res)

    figh = PyPlot.figure()
    PyPlot.grid()
    PyPlot.plot( t, q[1, :] )

    return figh
end

model = create_parametrized_model()

(forest,
integrator_indices,
tree_node_vecvec,
evaluation_indices_vecvec) = tenedyn.compile(model)

@time fcn_body = tenedyn.model_to_string(
    model,
    forest,
    integrator_indices,
    tree_node_vecvec,
    evaluation_indices_vecvec)

function create_callable_as_closure(ode_fcn_str, model)
    parsed_expr             = Meta.parse(ode_fcn_str)
    gen_ode_fcn             = eval(parsed_expr)
    ode_fcn(der_q, time, q) = gen_ode_fcn(der_q, time, q, model)

    return ode_fcn
end

ode_fcn = create_callable_as_closure(fcn_body, model)
# ode_res = sim2(ode_fcn);

end # module

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
include("helm_mg.jl")
include("jo_multigrid.jl")

ext: jl


text:
using MatrixCompletion

@testset "$(format("Sampling: BernoulliModel[VecOrMat]"))" begin
    sample_bernoulli_model0 = Sampler(BernoulliModel(0.8))
    tc1 = sample_bernoulli_model0.draw(ones(5,5))
    @test isa(tc1, Array{MaybeMissing{Float64},2}) || Array{Float64,2}
    tc2 = sample_bernoulli_model0.draw([1,2,3,4,5])
    @test isa(tc2, Array{MaybeMissing{Int64},1}) || isa(tc2,Array{Int64})
    sample_bernoulli_model1= provide(Sampler{BernoulliModel}(),rate = 0.8)
    tc3 = sample_bernoulli_model1.draw(ones(5,5))
    @test isa(tc3, Array{MaybeMissing{Float64},2}) || Array{Float64,2}
    tc4 = sample_bernoulli_model1.draw([1,2,3,4,5])
    @test isa(tc4, Array{MaybeMissing{Int64},1}) || isa(tc4,Array{Int64})
end

@testset "$(format("Sampling: UniformModel[VecOrMat]"))" begin
    #==================== Vector Case ====================#
    let
        tc1 = rand(10000)
        output = Sampler(UniformModel(0.5)).draw(tc1)
        @test count(x->ismissing(x),output) > 1000
        @test count(x->!ismissing(x),output) > 1000
        # stronger
        tc2 = rand(10000)
        output = Sampler(UniformModel(0.1)).draw(tc2)
        @test count(x->ismissing(x),output) >1000
        @test count(x->!ismissing(x),output) <1000
    end

    #==================== Matrix Case ====================#
    let
        @test 4 <= count(x -> !ismissing(x),Sampler(UniformModel(0.1)).draw(ones(10,10))) <= 10
        @test 10 < count(x -> ismissing(x), Sampler(UniformModel(0.1)).draw(ones(10,10)))
    end
    #==================== Factory Mode ====================#
    let
        sampler = provide(Sampler{UniformModel}(),rate = 0.5)
        tc1 = rand(10000)
        output = sampler.draw(tc1)
        @test count(x->ismissing(x),output) > 1000
        @test count(x->!ismissing(x),output) > 1000
        # stronger
        tc2 = rand(10000)
        sampler2 = provide(Sampler{UniformModel}(),rate = 0.1)
        output = sampler2.draw(tc2)
        @test count(x->ismissing(x),output) >1000
        @test count(x->!ismissing(x),output) <=1000
        @test 4 <= count(x -> !ismissing(x),sampler2.draw(ones(10,10))) <= 10
        @test 10 < count(x -> ismissing(x),sampler2.draw(ones(10,10)))
    end
end


text:
# Functions related to an individual's genome:
# - meiosis
# - mutation
# - genome creation
# - auxiliary functions

"""
    meiosis(genome, maternal)

Carry out meiosis on a genome (marked as maternal or not). Returns a haploid
gamete genome. (genome => array of chromosomes)
"""
function meiosis(genome::Array{Chromosome,1}, maternal::Bool, heterozygosity::Bool)
    #NOTE This function assumes that a genome is sorted, with all paternal chromosomes
    # in the first half and all maternals in the second (or vice versa).
    gametelength = Int(length(genome)/2)
    gamete = Array{Chromosome,1}(undef, gametelength)
    i = 1
    while i <= gametelength
        rand(Bool) ? g = i : g = i+gametelength
        # We do not deepcopy parental genes for performance reasons. However, that
        # means that if these genes are to be mutated, they must be deepcopied first.
        gamete[i] = heterozygosity ? LineageChromosome(genome[g].genes, maternal, genome[g].lineage) :
            DefaultChromosome(genome[g].genes, maternal)
        i += 1
    end
    gamete
end

"""
    gettraitdict(chromosomes, traitnames)

Convert a genome (an array of chromosomes) into a dict of traits and their values.
"""
function gettraitdict(chrms::Array{Chromosome, 1}, traitnames::Array{String, 1})
    #TODO can this be made more efficient? It's called really often...
    traitdict = Dict{String, Float64}()
    traits = Array{Trait,1}()
    nchrms = 0
    ngenes = 0
    for chrm in chrms
        nchrms += 1
        for gene in chrm.genes
            ngenes += 1
            append!(traits, gene.codes)
        end
    end
    for traitidx in eachindex(traitnames)
        wantedtraits = skipmissing(map(x -> x.value, filter(x -> x.nameindex == traitidx, traits)))
        traitdict[traitnames[traitidx]] = mean(wantedtraits)
        traitdict[traitnames[traitidx] * "sd"] = std(wantedtraits)
    end
    traitdict["ngenes"] = ngenes
    traitdict["nlnkgunits"] = nchrms
    traitdict
end

"""
    gettraitdictfast(chromosomes, traitnames)

Convert a genome (an array of chromosomes) into a dict of traits and their values.
This is an optimised version that can be run if `degpleiotropy` is 0 and `linkage` is "none".
"""
function gettraitdictfast(chrms::Array{Chromosome, 1}, traitnames::Array{String, 1})
    # Makes use of the fact that with `degpleiotropy == 0` and `linkage == "none"`,
    # there is exactly one trait per chromosome (one gene per chromosome and one trait per gene),
    # and the chromosomes are arranged in trait-order.
    genomesize = length(chrms)
    traitdict = Dict{String, Float64}()
    haploidlength = Int(genomesize/2)
    values = Array{Float64,}(undef,genomesize)
    for c in eachindex(chrms)
        @inbounds values[c] = chrms[c].genes[1].codes[1].value
    end
    for traitidx in eachindex(traitnames)
        wantedtraits = (values[traitidx], values[traitidx+haploidlength])
        traitdict[traitnames[traitidx]] = mean(wantedtraits)
        traitdict[traitnames[traitidx] * "sd"] = std(wantedtraits)
    end
    traitdict["ngenes"] = genomesize
    traitdict["nlnkgunits"] = genomesize
    traitdict
end

"""
    getseqsimilarity(seqone, seqtwo)

Compare two strings and return similarity.
"""
function getseqsimilarity(indgene::AbstractString, mategene::AbstractString)
    basediffs = 0
    #FIXME this should find out whether indgene or mategene is shorter and use that,
    # then count the length difference towards the base diff
    for i in eachindex(indgene)
        try
            indgene[i] != mategene[i] && (basediffs += 1) # alternatively use bioinformatics tools
        catch # e.g., in case of differently long genes
            basediffs += 1
        end
    end
    seqidentity = 1 - (basediffs / length(indgene))
end

"""
    getseqs(genome, traitname)

Find and return the sequences of genes that code for the given trait.
"""
function getseqs(genome::Array{Chromosome, 1}, traitname::String)
    seqs = String[]
    traitidx = findfirst(x -> x == traitname, setting("traitnames"))
    for chrm in genome
        for gene in chrm.genes
            if any(x -> x.nameindex == traitidx, gene.codes)
                setting("compressgenes") ? push!(seqs, num2seq(gene.sequence)) : push!(seqs, gene.sequence)
            end
        end
    end
    seqs
end

"""
    iscompatible(mate, individual)

Check to see whether two individual organisms are reproductively compatible by
comparing gene sequences (which sequence depends on setting("speciation")).
"""
function iscompatible(mate::Individual, ind::Individual)
    mate.lineage != ind.lineage && return false
    indgenes, mategenes = "", ""
    if setting("speciation") == "off"
        return true # all conspecifics are compatible, no speciation can happen
    elseif setting("speciation") == "neutral"
        # default, use a non-coding sequence for mutation-order speciation
        indgenes = join(getseqs(ind.genome, "compat"))
        mategenes = join(getseqs(mate.genome, "compat"))
    elseif setting("speciation") == "ecological"
        # use all coding sequences for ecological speciation
        for tr in setting("traitsforecospec")
            (tr == "compat") && continue
            indgenes *= join(getseqs(ind.genome, tr))
            mategenes *= join(getseqs(mate.genome, tr))
        end
    end
    seqidentity = getseqsimilarity(indgenes, mategenes)
    seqidentity >= ind.traits["seqsimilarity"]
end

"""
    seq2num(sequence)

Convert a DNA base sequence (a string) into binary and then into an integer.
This saves memory.
"""
#XXX Actually, this hardly seems to make a difference :-(
function seq2num(sequence::String)
    # This function effectively uses the same technique as seq2bignum, but it skips the
    # intermediate allocations and should therefore be more efficient.
    num::Int64 = 0  # Int64 allows for max length of 21bp
    for b in eachindex(sequence)
        if sequence[end+1-b] == 'a'
            num += 2^(3*(b-1)) * 4 # b'100'
        elseif sequence[end+1-b] == 'c'
            num += 2^(3*(b-1)) * 5 # b'101'
        elseif sequence[end+1-b] == 'g'
            num += 2^(3*(b-1)) * 6 # b'110'
        elseif sequence[end+1-b] == 't'
            num += 2^(3*(b-1)) * 7 # b'111'
        end
    end
    num
end

function intseq(n::Int)
    #XXX doesn't seem to work either? I'm starting to think I'm optimising in the wrong place :-(
    num::Int64 = 0  # Int64 allows for max length of 21bp
    for b in 0:(n-1)
        # basically like seq2bignum(), but skips all intermediate allocations
        num += 2^(3*b) * rand(4:7)
    end
    num
end

"""
    seq2bignum(sequence)

Convert a DNA base sequence (a string) into binary and then into a BigInt (for
larger genes). This saves memory.
"""
function seq2bignum(sequence::String)
    bases = "acgt"
    binary = ""
    for base in sequence
        binary *= string(findfirst(x -> x == base, bases) + 3, base = 2)
    end
    parse(BigInt, binary, base = 2) # BigInt allows arbitrary-length sequences
end

"""
    num2seq(n)

Convert an integer into binary and then into a DNA base sequence string.
"""
function num2seq(n::Integer)
    bases = "acgt"
    binary = string(n, base = 2)
    sequence = ""
    for i in 1:3:(length(binary) - 2)
        sequence *= string(bases[parse(Int, binary[i:(i + 2)], base = 2) - 3])
    end
    sequence
end

"""
    createtraits()

Create an array of trait objects generated from the default trait values (with a
random offset).
"""
function createtraits()
    #XXX this is all very ugly. (case/switch w/ v. 2.0+?)
    traitnames = setting("traitnames")
    traits = Array{Trait,1}(undef,length(traitnames))
    # exponential distributions of body sizes:
    repoffset = setting("maxrepsize") - setting("minrepsize")
    seedoffset = setting("maxseedsize") - setting("minseedsize")
    tempoffset = setting("maxtemp") - setting("mintemp")
    precoffset = setting("maxprec") - setting("minprec")
    #FIXME this needs to be changed to account for allometric relationships
    repsize, seedsize = 0, 0
    while true
        repsize = exp(setting("minrepsize") + repoffset * rand())
        seedsize = exp(setting("minseedsize") + seedoffset * rand())
        repsize > seedsize && break
    end
    for idx in eachindex(traitnames)
        if occursin("dispshape", traitnames[idx])
            traits[idx] = Trait(idx, rand() * setting("dispshape"))
        elseif occursin("dispmean", traitnames[idx])
            traits[idx] = Trait(idx, rand() * setting("dispmean"))
        elseif occursin("precopt", traitnames[idx])
            traits[idx] = Trait(idx, setting("minprec") + rand() * precoffset)
        elseif occursin("prectol", traitnames[idx])
            traits[idx] = Trait(idx, rand() * setting("maxbreadth"))
        elseif occursin("repsize", traitnames[idx])
            traits[idx] = Trait(idx, repsize)
        elseif occursin("seqsimilarity", traitnames[idx]) && setting("fixtol")
            traits[idx] = Trait(idx, setting("tolerance"))
        elseif occursin("seqsimilarity", traitnames[idx]) && !setting("fixtol")
            traits[idx] = Trait(idx, rand()) # assortative mating might evolve
        elseif occursin("seedsize", traitnames[idx])
            traits[idx] = Trait(idx, seedsize)
        elseif occursin("tempopt", traitnames[idx])
            traits[idx] = Trait(idx, setting("mintemp") + rand() * tempoffset)
        elseif occursin("temptol", traitnames[idx])
            traits[idx] = Trait(idx, rand() * setting("maxbreadth"))
        else
            traits[idx] = Trait(idx, rand())
        end
    end
    traits
end

"""
    creategenes(ngenes, traits)

Randomly create a given number of gene objects, with their base sequence and
associated traits. Returns the result as an array of AbstractGenes.
"""
function creategenes(ngenes::Int, traits::Array{Trait,1})
    genes = Array{AbstractGene,1}(undef, ngenes)
    bases = ['a','c','g','t']
    compatidx = findfirst(x -> x == "compat", setting("traitnames"))
    # initialise each gene with an arbitrary sequence
    for i in 1:ngenes
        if setting("compressgenes") #default
            genes[i] = Gene(intseq(setting("smallgenelength")), Trait[])
        else
            sequence = String(rand(bases, setting("smallgenelength")))
            genes[i] = StringGene(sequence, Trait[])
        end
    end
    # assign traits to the genes
    if setting("degpleiotropy") == 0
        # Disable polygenic inheritance and pleiotropy: make sure one gene codes for one trait.
        for tr in eachindex(traits)
            (tr == compatidx) && (genes[tr] = createcompatgene(bases, compatidx))
            genes[tr].codes = [traits[tr]]
        end
    else # allow for pleiotropy as well as polygenic inheritance
        for trait in traits
            (trait.nameindex == compatidx) && continue
            # We're actually setting polygenic inheritance here, rather than pleiotropy.
            # Pleiotropy is introduced indirectly, because a higher number of coding genes
            # increases the likelihood of picking a gene that already codes for other traits.
            ncodinggenes = rand(Geometric(1 - setting("degpleiotropy"))) + 1
            codinggenes = rand(genes,ncodinggenes)
            for gene in codinggenes
                push!(gene.codes,deepcopy(trait))
            end
        end
        push!(genes, createcompatgene(bases, compatidx))
    end
    genes
end

"""
    createcompatgene(bases, compatidx)

Create the gene that is used to determine reproductive compatibility.
"""
function createcompatgene(bases::Array{Char,1}, compatidx::Int)::AbstractGene
    # Note: we only need big genes for the compatibility gene, because we need a longer base
    # sequence than offered by `smallgenelength` if we want to do phylogenetic analyses.
    setting("usebiggenes") ? seql = setting("biggenelength") : seql = setting("smallgenelength")
    cseq = String(rand(bases, seql))
    ctrait = [Trait(compatidx, 0.5)]
    if setting("compressgenes")
        if setting("usebiggenes")
            return BigGene(seq2bignum(cseq), ctrait)
        else
            return Gene(seq2num(cseq), ctrait)
        end
    else
        return StringGene(cseq, ctrait)
    end
end

"""
    createchrms(nchrms, genes)

Randomly distribute the passed genes into the given number of haploid chromosomes.
Returns a diploid genome (array of chromosome objects).
"""
function createchrms(nchrms::Int,genes::Array{AbstractGene,1},lineage::String)::Array{Chromosome,1}
    chromosomes = Array{Chromosome,1}(undef,nchrms*2)
    if nchrms == 1 # linkage == "full"
        chromosomes[1] = @Chromosome(genes, true, lineage)
        chromosomes[2] = @Chromosome(deepcopy(genes),false, lineage)
    elseif nchrms == length(genes) # linkage == "none"
        for g in eachindex(genes)
            chromosomes[g] = @Chromosome([genes[g]], true, lineage)
            chromosomes[nchrms+g] = @Chromosome([deepcopy(genes[g])], false, lineage)
        end
    else # linkage == "random"
        chrmsplits = sort(rand(1:length(genes), nchrms-1))
        for chr in 1:nchrms
            if chr==1 # first chromosome
                cgenes = genes[1:chrmsplits[chr]]
            elseif chr==nchrms # last chromosome
                cgenes = genes[(chrmsplits[chr-1]+1):end]
            else
                cgenes = genes[(chrmsplits[chr-1]+1):chrmsplits[chr]]
            end
            chromosomes[chr] = @Chromosome(cgenes, true, lineage)
            chromosomes[nchrms+chr] = @Chromosome(deepcopy(cgenes), false, lineage)
        end
    end
    chromosomes
end

"""
    varyalleles!(genes, locivar)

Mutate gene traits in the passed array of genes.
"""
function varyalleles!(genes::Array{AbstractGene, 1}, locivar::Float64)
    locivar == 0 && return
    for gene in genes
        mutate!(gene.codes, locivar)
    end
end

"""
    varyalleles!(chromosomes, locivar)

Mutate gene traits in the passed array of chromosomes.
"""
function varyalleles!(chrms::Array{Chromosome, 1}, locivar::Float64)
    locivar == 0 && return
    for chrm in chrms
        varyalleles!(chrm.genes, locivar)
    end
end

"""
    mutate!(traits, locivar)

Loop over an array of traits, mutating each value in place along a normal distribution.
`locivar` can be used to scale the variance of the normal distribution used to draw new
trait values (together with `setting(phylconstr]`).
"""
function mutate!(traits::Array{Trait, 1}, locivar::Float64 = 1.0)
    setting("phylconstr") * locivar == 0 && return
    for trait in traits
        traitname = setting("traitnames")[trait.nameindex]
        occursin("seqsimilarity", traitname) && setting("fixtol") && continue
        oldvalue = trait.value
        occursin("tempopt", traitname) && (oldvalue -= 273)
        while oldvalue <= 0 # make sure sd of Normal dist != 0
            oldvalue = abs(rand(Normal(0,0.01)))
        end
        newvalue = rand(Normal(oldvalue, oldvalue * setting("phylconstr") * locivar))
        #FIXME there are no traits with "prob" in their name?!
        # Is this meant to cap compat/seqsimilarity/selfing?
        (newvalue > 1 && occursin("prob", traitname)) && (newvalue=1.0)
        while newvalue <= 0
            newvalue = rand(Normal(oldvalue, oldvalue * setting("phylconstr") * locivar))
        end
        occursin("tempopt", traitname) && (newvalue += 273)
        trait.value = newvalue
    end
end

"""
    mutate!(individual, temp)

Mutate an individual's genome (sequence and traits) in place.
"""
function mutate!(ind::Individual, temp::Float64)
    muts = setting("mutationrate") * exp(-act/(boltz*temp))
    nmuts = rand(Poisson(muts))
    nmuts == 0 && return
    chrmidcs = rand(eachindex(ind.genome), nmuts)
    for c in chrmidcs
        # We need to deepcopy genes here, because they are not automatically deepcopied during
        # reproduction for performance reasons
        ind.genome[c] = deepcopy(ind.genome[c])
        length(ind.genome[c].genes) == 0 && continue
        g = rand(eachindex(ind.genome[c].genes))
        gseq = ind.genome[c].genes[g].sequence
        setting("compressgenes") ? charseq = collect(num2seq(gseq)) : charseq = collect(gseq)
        i = rand(eachindex(charseq))
        newbase = rand(collect("acgt"),1)[1]
        while newbase == charseq[i]
            newbase = rand(collect("acgt"),1)[1]
        end
        charseq[i] = newbase
        mutate!(ind.genome[c].genes[g].codes)
        if setting("compressgenes")
            if length(charseq) > 21
                ind.genome[c].genes[g].sequence = seq2bignum(String(charseq))
            else
                ind.genome[c].genes[g].sequence = seq2num(String(charseq))
            end
        else
            ind.genome[c].genes[g].sequence = String(charseq)
        end
    end
    if setting("degpleiotropy") == 0 && setting("linkage") == "none"
        ind.traits = gettraitdictfast(ind.genome, setting("traitnames"))
    else
        ind.traits = gettraitdict(ind.genome, setting("traitnames"))
    end
end

"""
    mutate!(patch)

Mutate all seed individuals in a patch.
"""
function mutate!(patch::Patch)
    setting("mode") == "zosterops" ? temp = setting("bodytemp") : temp = patch.temp
    for ind in patch.seedbank
        mutate!(ind, temp)
    end
end

"""
    mutate!(world)

Mutate the world. (That sounds scary!)
"""
function mutate!(world::Array{Patch, 1})
    for patch in world
        (patch.isisland || !setting("static")) && mutate!(patch)
    end
end

max_stars_repo_name: CCTB-Ecomods/gemm

max_issues_repo_name: CCTB-Ecomods/gemm

max_forks_repo_name: CCTB-Ecomods/gemm


text:
export MtlError

const MTLError = Ptr{NsError}

mutable struct MtlError <: Exception
	code::NsInteger
	domain::String
	userinfo::String

	ptr::MTLError
end

Base.convert(::Type{MTLError}, err::MtlError) = err.ptr
Base.unsafe_convert(::Type{MTLError}, err::MtlError) = convert(MTLError, err.ptr)

function MtlError(err::MTLError)
	code = mtErrorCode(err)
	domain = mtErrorDomain(err) |> unsafe_string
	userinfo = mtErrorUserInfo(err) |> unsafe_string

	obj = MtlError(code, domain, userinfo, err)
	finalizer(unsafe_destroy!, obj)
	return obj
end

function unsafe_destroy!(err::MtlError)
	mtErrorRelease(err)
end

description(err::MtlError) = unsafe_string_maybe(mtErrorLocalizedDescription(err))
recoverySuggestion(err::MtlError) = unsafe_string_maybe(mtErrorLocalizedRecoverySuggestion(err))
failureReason(err::MtlError) = unsafe_string_maybe(mtErrorLocalizedFailureReason(err))

function recoveryOptions(err::MtlError)
	options = mtErrorLocalizedRecoveryOptions(err)
	opts = Vector{String}();
	for i=0:typemax(Int)
		_opt = Base.unsafe_load(options + i * sizeof(Ptr{NsError}))
		_opt == C_NULL && break
		push!(opts, unsafe_string(_opt))
	end
	Base.Libc.free(options)
	return opts
end

Base.show(io::IO, ::MIME"text/plain", err::MtlError) = _show(io, err)
Base.show(io::IO, err::MtlError) = _show(io, err)

function _show(io::IO,  err::MtlError)
	println(io, "MtlError (Error in Metal Runtime):")
	println(io, " code     : ", err.code)
	println(io, " domain   : ", err.domain)
	println(io, " userinfo : ", replace(err.userinfo, "\\n"=>"\n"))

	isempty(description(err)) || println(io, "Description:", description(err))
	isempty(recoverySuggestion(err)) || println(io, "Suggestion:", recoverySuggestion(err))
	isempty(failureReason(err)) || println(io, "Failure Reason:", failureReason(err))
	isempty(recoveryOptions(err)) || println(io, "Recovery Options:", recoveryOptions(err))
end


text:
# Define default values for keyword arguments.
mutable struct SLDefs
  port::Int64
  gap::Int64
  kai::Int64
  mode::String
  refresh::Real
  x_on_err::Bool

  SLDefs( port::Int64,
          gap::Int64,
          kai::Int64,
          mode::String,
          refresh::Real,
          x_on_err::Bool
          ) = new(port, gap, kai, mode, refresh, x_on_err)
end

mutable struct FiltDefs
  fl::Float64
  fh::Float64
  np::Int64
  rp::Int64
  rs::Int64
  rt::String
  dm::String

  FiltDefs( fl::Float64,
            fh::Float64,
            np::Int64,
            rp::Int64,
            rs::Int64,
            rt::String,
            dm::String) = new(fl, fh, np, rp, rs, rt, dm)
end

mutable struct KWDefs
  SL::SLDefs
  Filt::FiltDefs
  comp::UInt8
  evw::Array{Real,1}
  fmt::String
  full::Bool
  mag::Array{Float64,1}
  nd::Real
  nev::Int64
  n_zip::Int64
  nx_add::Int64
  nx_new::Int64
  opts::String
  pha::String
  prune::Bool
  rad::Array{Float64,1}
  reg::Array{Float64,1}
  si::Bool
  src::String
  to::Int64
  v::Int64
  w::Bool
  y::Bool
end

"""
    SeisIO.KW

An immutable structure containing default keyword argument values in SeisIO.
Arguments that accept keywords in SeisIO.KW use the default values when a
keyword isn't specified.

### Main Module Keywords

| KW       | Default    | Allowed Data Types | Meaning                        |
|----------|:-----------|:-------------------|:-------------------------------|
| comp     | 0x00       | UInt8              | compress data on write?[^1]    |
| fmt      | "miniseed" | String             | request data format            |
| full     | false      | Bool               | read full headers?             |
| nd       | 1          | Real               | number of days per subrequest  |
| n_zip    | 100000     | Int64              | compress if length(x) > n_zip  |
| nx_add   | 360000     | Int64              | minimum length increase of an  |
|          |            |                    |  undersized data array         |
| nx_new   | 8640000    | Int64              | number of samples allocated    |
|          |            |                    |   for a new data channel       |
| opts     | ""         | String             | user-specified options[^2]     |
| prune    | true       | Bool               | call prune! after get_data?    |
| rad      | []         | Array{Float64,1}   | radius search: `[center_lat,`  |
|          |            |                    |   `center_lon, r_min, r_max]`  |
|          |            |                    |   in decimal degrees (°)       |
| reg      | []         | Array{Float64,1}   | geographic search region:      |
|          |            |                    |   `[min_lat, max_lat,`         |
|          |            |                    |    `min_lon, max_lon,`         |
|          |            |                    |   `min_dep, max_dep]`          |
|          |            |                    |   lat, lon in degrees (°)      |
|          |            |                    |   dep in km with down = +      |
| si       | true       | Bool               | autofill request station info? |
| to       | 30         | Int64              | timeout (s) for web requests   |
| v        | 0          | Int64              | verbosity                      |
| w        | false      | Bool               | write requests to disc?        |
| y        | false      | Bool               | sync after web requests?       |

[^1]: If `comp == 0x00`, never compress data; if `comp == 0x01`, only compress channel `i` if `length(S.x[i]) > KW.n_zip`; if `comp == 0x02`, always compress data.
[^2]: Format as for an http request request URL, e.g. "szsrecs=true&repo=realtime" for FDSN (note: string should not begin with an ampersand).

### Quake Submodule Keywords

| KW       | Default    | Allowed Data Types | Meaning                        |
|----------|:-----------|:-------------------|:-------------------------------|
| evw      | [600,      | Array{Real,1}      | search for events in window    |
|          |  600]      |                    |   [ot-abs(t1), ot+abs(t2)]     |
| mag      | [6.0, 9.9] | Array{Float64,1}   | search magitude range          |
| nev      | 0          | Int64              | number of events per query     |
|          |            |                    | (if nev=0, return all matches) |
| pha      | "P"        | String             | phases to get (comma-separated |
|          |            |                    |   list; use "ttall" for all)   |

### Substructures

    SeisIO.KW.SL: Seedlink-specific keyword defaults. SeedLink also uses some
general keywords.

| Name        | Default | Type            | Description                       |
|:------------|:--------|:----------------|:----------------------------------|
| gap         | 3600    | Real            | max. gap since last packet [s]    |
| kai         | 600     | Real            | keepalive interval [s]            |
| mode        | "DATA"  | String          | TIME, DATA, or FETCH              |
| port        | 18000   | Int64           | port number                       |
| refresh     | 20      | Real            | base refresh interval [s]         |
| xonerr      | true    | Bool            | exit on error?                    |

SeisIO.KW.Filt: Defaults parameters for time-series filtering.

| Name  | Default       | Type    | Description                         |
|:------|:--------------|:--------|:------------------------------------|
| fl    | 1.0           | Float64 | lower corner frequency [Hz] [^1]    |
| fh    | 15.0          | Float64 | upper corner frequency [Hz] [^1]    |
| np    | 4             | Int64   | number of poles                     |
| rp    | 8             | Int64   | pass-band ripple (dB)               |
| rs    | 30            | Int64   | stop-band ripple (dB)               |
| rt    | "Bandpass"    | String  | response type (type of filter)      |
| dm    | "Butterworth" | String  | design mode (name of filter)        |

[^1]: Remember the (counter-intuitive) convention that the lower corner frequency (fl) is used in a Highpass filter, and fh is used in a Lowpass filter. This convention is preserved in SeisIO.

"""
const KW = KWDefs(
           SLDefs(18000,    # port::Int64
                   3600,    # gap::Int64
                    600,    # kai::Int64
                 "DATA",    # mode::String
                   20.0,    # refresh::Real
                   true ),  # x_on_err::Bool

           FiltDefs(1.0,    # fl::Float64
                   15.0,    # fh::Float64
                      4,    # np::Int64
                      8,    # rp::Int64
                     30,    # rs::Int64
             "Bandpass",    # rt::String
          "Butterworth" ),  # dm::String

                     0x00,  # comp::Bool
               [600, 600],  # evw::Real
               "miniseed",  # fmt::String
                    false,  # full::Bool
        Float64[6.0, 9.9],  # mag::Array{Float64,1}
                        1,  # nd::Real
                        0,  # nev::Int64
                   100000,  # n_zip::Int64
                   360000,  # nx_add::Int64
                  8640000,  # nx_new::Int64
                        "", # opts::String
                       "P", # pha::String
                      true, # prune::Bool
                 Float64[], # rad: Array{Float64,1}
                 Float64[], # reg: Array{Float64,1}
                      true, # si::Bool
                    "IRIS", # src::String
                        30, # to::Int64
                         0, # v::Int (verbosity)
                     false, # w::Bool (write to disk)
                     false) # y::Bool (syc)


text:
### generic DataSet

"""
    mix(f, ϕ,    ds::DataSet)
    mix(f, ϕ, θ, ds::DataSet)

Compute the mixed `(f°, ϕ°)` from the unlensed field `f` and lensing potential
`ϕ`, given the definition of the mixing matrices in `ds` evaluated at parameters
`θ` (or at fiducial values if no `θ` provided).
"""
mix(f, ϕ, ds::DataSet) = mix(f,ϕ,(;),ds)
function mix(f, ϕ, θ, ds::DataSet)
    @unpack D,G,L = ds(θ)
    L(ϕ)*D*f, G*ϕ
end

"""
    unmix(f°, ϕ°,    ds::DataSet)
    unmix(f°, ϕ°, θ, ds::DataSet)

Compute the unmixed/unlensed `(f, ϕ)` from the mixed field `f°` and mixed
lensing potential `ϕ°`, given the definition of the mixing matrices in `ds`
evaluated at parameters `θ` (or at fiducial values if no `θ` provided).
"""
unmix(f°, ϕ°, ds::DataSet) = unmix(f°,ϕ°,(;),ds)
function unmix(f°, ϕ°, θ, ds::DataSet)
    @unpack D,G,L = ds(θ)
    ϕ = G\ϕ°
    D\(L(ϕ)\f°), ϕ
end

@doc doc"""
    lnP(t, fₜ, ϕₜ,    ds::DataSet)
    lnP(t, fₜ, ϕₜ, θ, ds::DataSet)

Compute the log posterior probability in the joint parameterization as a
function of the field, $f_t$, the lensing potential, $\phi_t$, and possibly some
cosmological parameters, $\theta$. The subscript $t$ can refer to either a
"time", e.g. passing `t=0` corresponds to the unlensed parametrization and `t=1`
to the lensed one, or can be `:mix` correpsonding to the mixed parametrization.
In all cases, the arguments `fₜ` and `ϕₜ` should then be $f$ and $\phi$ in
that particular parametrization.

If any parameters $\theta$ are provided, we also include the determinant terms for
covariances which depend on $\theta$. In the mixed parametrization, we
also include any Jacobian determinant terms that depend on $\theta$.

The argument `ds` should be a `DataSet` and stores the masks, data, etc...
needed to construct the posterior.
"""
lnP(t, fₜ, ϕₜ,    ds::DataSet) = lnP(Val(t), fₜ, ϕₜ, (;), ds)
lnP(t, fₜ, ϕₜ, θ, ds::DataSet) = lnP(Val(t), fₜ, ϕₜ,  θ,  ds)

function lnP(::Val{t}, fₜ, ϕ, θ, ds::DataSet) where {t}

    @unpack Cn,Cf,Cϕ,L,M,B,d = ds

    f,f̃ = t==0 ? (fₜ, L(ϕ)*fₜ) : (L(ϕ)\fₜ, fₜ)
    Δ = d - M(θ)*B(θ)*f̃ - nonCMB_data_components(θ,ds)
    (
        -(1//2) * (
            Δ'*pinv(Cn(θ))*Δ + logdet(Cn,θ) +
            f'*pinv(Cf(θ))*f + logdet(Cf,θ) +
            ϕ'*pinv(Cϕ(θ))*ϕ + logdet(Cϕ,θ)
        )
        + lnPriorθ(θ,ds)
    )

end

function lnP(::Val{:mix}, f°, ϕ°, θ, ds::DataSet)
    lnP(Val(0), unmix(f°,ϕ°,θ,ds)..., θ, ds) - logdet(ds.D,θ) - logdet(ds.G,θ)
end

# can be specialized for specific DataSet types:
nonCMB_data_components(θ, ds::DataSet) = 0
lnPriorθ(θ, ds::DataSet) = 0

## NoLensingDataSet

lnP(   f,    ds::NoLensingDataSet) = lnP(f, (;), ds)
lnP(_, f, θ, ds::NoLensingDataSet) = lnP(f, θ,   ds)

function lnP(f, θ, ds::NoLensingDataSet)

    @unpack Cn,Cf,M,B,d = ds

    Δ = d - M(θ)*B(θ)*f - nonCMB_data_components(θ,ds)
    (
        -(1//2) * (
            Δ'*pinv(Cn(θ))*Δ + logdet(Cn,θ) +
            f'*pinv(Cf(θ))*f + logdet(Cf,θ)
        )
        + lnPriorθ(θ,ds)
    )

end


text:
function mean!(μ::VecIO, r::VecI, X::MatI)
    N = sum(r)
    if iszero(N)
        @simd for m in eachindex(μ)
            @inbounds μ[m] = 0.0
        end
    else
        BLAS.gemv!('T', 1.0, X, r, 0.0, μ)
        @simd for m in eachindex(μ)
            @inbounds μ[m] /= N
        end
    end
    return nothing
end

function mean!(μ::MatIO, R::MatI, X::MatI, N::VecI; init::Bool=false)
    if init
        @inbounds for i in eachindex(N)
            N[i] = sum(view(R,i,:))
        end
    end
    BLAS.gemm!('N', 'N', 1.0, R, X, 0.0, μ)
    for k in axes(μ, 1)
        @inbounds Nk = N[k]
        if iszero(Nk)
            @simd for m in axes(μ, 2)
                @inbounds μ[k,m] = 0.0
            end
        else
            @simd for m in axes(μ, 2)
                @inbounds μ[k,m] /= Nk
            end
        end
    end
    return nothing
end

function deviation!(Δ::MatIO, X::MatI, μ::VecI)
    n, m = size(X)
    for j in eachindex(1:m)
        @inbounds μj = μ[j]
        @simd for i in eachindex(1:n)
            @inbounds Δ[i,j] = X[i,j] - μj
        end
    end
    return nothing
end

function covariance!(C::MatIO, X::MatI, μ::VecI, r::VecI, Δ::MatB)
    deviation!(Δ, X, μ)
    for n in eachindex(r)
        @inbounds rn = sqrt(r[n])
        @simd for m in eachindex(μ)
            @inbounds Δ[n,m] *= rn
        end
    end
    BLAS.syrk!('U', 'T', 1.0, Δ, 0.0, C)
    @inbounds for j in eachindex(μ), i in 1:j-1
        C[j,i] = C[i,j]
    end
    return nothing
end


text:
# Generating solarwind wave fluctuations.

export BackgroundVariable, BoundaryVariable
export generate_signal, save_signal

using DelimitedFiles, UnPack

abstract type Variable end

struct BackgroundVariable{T} <: Variable
   n::T
   T::T
   Vx::T
   Vy::T
   Vz::T
   Bx::T
   By::T
   Bz::T
end

struct BoundaryVariable{T} <: Variable
   t::AbstractVector{T}
   n::Vector{T}
   T::Vector{T}
   Vx::Vector{T}
   Vy::Vector{T}
   Vz::Vector{T}
   Bx::Vector{T}
   By::Vector{T}
   Bz::Vector{T}
end

"""
    generate_signal(varBG; fsignal=1.0, dv=1e3, fsample=5.0, tstart=0.0, tend=10.0,
       signal=:alfven, dir="xyz", model=:vlasiator)

Generate a fluctuation signal as an upstream input file to plasma models.
Currently only support input perturbation in velocity `dv`.

# Optional Arguments
- `fsignal=1.0`: fluctuation signal frequency in [Hz].
- `fsample=5.0`: sampling frequency in [Hz].
- `tstart=0.0`: starting time in [s].
- `tend=10.0`: end time in [s].
- `dv=1e3`: velocity perturbation magnitude in [m/s].
- `dn=1`: mass density perturbation magnitude in [amu/cc].
- `signal=:alfven`: [:alfven, :fast].
- `dir="xyz"`: direction(s) for the applied Alfvénic perturbations in Cartesian coordinates,
which  can be any combination of "x", "y" and "z".
- `model=:vlasiator`: model specific settings. In Vlasiator, the background magnetic field
includes both a dipole and a constant. Therefore in the solarwind input files magnetic field
values should not contain the constant field.
"""
function generate_signal(varBG::BackgroundVariable; fsignal=1.0, fsample=5.0,
   tstart=0.0, tend=10.0, signal=:alfven, dir="xyz", dv=1e3, dn=1, model=:vlasiator)

   n0, T0, Vx0, Vy0, Vz0, Bx0, By0, Bz0 = varBG.n, varBG.T, varBG.Vx, varBG.Vy, varBG.Vz,
      varBG.Bx, varBG.By, varBG.Bz

   B0 = hypot(Bx0, By0, Bz0)
   μ₀ = 4π*1e-7
   mp = 1.6726219*1e-27
   VA = B0 / √(μ₀ * n0 * mp) # background Alfven speed

   len = floor(Int, (tend - tstart) * fsample) + 1
   t = range(tstart, tend, length=len)

   n = fill(n0, len)
   T = fill(T0, len)
   Vx, Vy, Vz = fill(Vx0, len), fill(Vy0, len), fill(Vz0, len)
   if model == :vlasiator
      Bx, By, Bz = zeros(len), zeros(len), zeros(len)
   else
      Bx, By, Bz = fill(Bx0, len), fill(By0, len), fill(Bz0, len)
   end

   if signal == :alfven
      δV = @. dv * sinpi(2 * fsignal * (t - tstart))
      δB = @. B0 / VA * δV
      if occursin("x", dir)
         Vx .+= δV
         Bx .+= δB
      end
      if occursin("y", dir)
         Vy .+= δV
         By .+= δB
      end
      if occursin("z", dir)
         Vz .+= δV
         Bz .+= δB
      end

   elseif signal == :fast
      δn = @. dn * 1e6 * sinpi(2 * fsignal * (t - tstart))
      n .+= δn
   end

   v = BoundaryVariable(t, n, T, Vx, Vy, Vz, Bx, By, Bz)
end

"Save generated signal data `var` to `file`."
function save_signal(var::BoundaryVariable; file="sw.dat", sigdigits=5)
   @unpack t, n, T, Vx, Vy, Vz, Bx, By, Bz = var
   open(file, "w") do io
      write(io, "# time density temperature Vx Vy Vz Bx By Bz\n")
      writedlm(io, round.([t n T Vx Vy Vz Bx By Bz]; sigdigits))
   end
end


text:
module Mimikkyu

using Printf, LaTeXStrings

import PyPlot as plt

export @L_str, init_rc!,
       save_pdf, save_png,
       xticks!, yticks!, zticks!

function init_rc!()
    rcParams = plt.PyDict(plt.matplotlib."rcParams")
    pgfpreamble = string("\\usepackage{siunitx}",
                         "\\usepackage{mhchem}",
                         "\\usepackage{unicode-math}",
                         "\\setmainfont{XITS}",
                         "\\setmathfont{XITS Math}")
    # Latex
    rcParams["pgf.rcfonts"] = false
    rcParams["pgf.preamble"] = pgfpreamble
    rcParams["text.usetex"] = true
    rcParams["text.latex.preamble"] = "\\usepackage{siunitx}\\usepackage{mhchem}"
    # Font
    rcParams["font.size"] = 14
    rcParams["font.family"] = "serif"
    # Axes
    rcParams["axes.linewidth"] = 1.5
    rcParams["xtick.major.width"] = 1.5
    rcParams["ytick.major.width"] = 1.5
    # Grid
    rcParams["grid.alpha"] = 0.5
    rcParams["grid.linewidth"] = 1.5
    # Lines
    rcParams["lines.linewidth"] = 2.0
    # Legend
    rcParams["legend.framealpha"] = 0.0
    return nothing
end

include("./ticks.jl")

save_pdf(f::String, fig) = fig.savefig(f, format="pdf", dpi=600, backend="pgf")
save_png(f::String, fig) = fig.savefig(f, format="png", dpi=600, backend="pgf")

end # module


text:
function check_syntax(line)
    char_pair = Dict(
        '(' => ')',
        '{' => '}',
        '[' => ']',
        '<' => '>'
    )

    open_chars = Char[]
    close_chars = Char[]
    for (i, char) in enumerate(line)
        if char in keys(char_pair)
            push!(open_chars, char)
            push!(close_chars, char_pair[char])
        else
            if char == close_chars[end]
                pop!(open_chars)
                pop!(close_chars)
            else
                return "corrupt", line[i:end]
            end
        end
    end

    if !isempty(close_chars)
        return "incomplete", join(reverse(close_chars))
    end

    "complete", nothing
end

function score_error(error)
    points = Dict(
        ')' => 3,
        ']' => 57,
        '}' => 1197,
        '>' => 25137
    )
    points[error[1]]
end

function score_completion(completion)
    points = Dict(
        ')' => 1,
        ']' => 2,
        '}' => 3,
        '>' => 4
    )
    score = 0
    for char in completion
        score *= 5
        score += points[char]
    end
    score
end

lines = readlines("day10/lines.txt")
results = check_syntax.(lines)

errors = [error for (status, error) in results if status == "corrupt"]
p1 = sum(score_error.(errors))

@show p1

completions = [completion for (status, completion) in results if status == "incomplete"]
p2 = sort(score_completion.(completions))[convert(Int, ceil(length(completions) / 2))]

@show p2


text:
using WeaveUtil

function bld()
	buildJMD("main.jmd")
end


text:
using SparseGrids
using Test

@testset "N_q^D" begin
    @testset "Only one element with D == q" begin
        ndq = listNdq(3, 3)

        @test ndq == [ones(Int, 3)]
    end

    @testset "Split a dice" begin
        ndq = listNdq(2, 6)

        @test ndq == [[5, 1], [4, 2], [3, 3], [2, 4], [1, 5]]
    end
end


text:
facts("KDTrees.inball") do

context("KDTrees.inball.point") do

for rn in [true, false]
    data = [0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0;
            0.0 0.0 1.0 1.0 0.0 0.0 1.0 1.0;
            0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0] # 8 node cube

    tree = KDTree(data, leafsize = 2, reorder = rn)

    idxs = inball(tree, [1.1, 1.1, 1.1], 0.2)
    @fact idxs --> [8] # Only corner 8 at least 0.2 distance away from [1.1, 1.1, 1.1]

    idxs = inball(tree, [0.0, 0.0, 0.5], 0.6)
    @fact allin(idxs, [1, 2]) --> true# Corner 1 and 2 at least 0.6 distance away from [0.0, 0.0, 0.5]

    idxs = inball(tree, [0, 0, 0], 0.6)
    @fact idxs --> [1]

    idxs = inball(tree, [1//3, 1//3, 1//3], 1)
    @fact allin(idxs, [1, 2, 3, 5]) --> true

    idxs = inball(tree, [0.5, 0.5, 0.5], 0.2)
    @fact idxs --> [] #

    idxs = inball(tree, [0.5, 0.5, 0.5], 1.0)
    @fact allin(idxs, [1, 2, 3, 4, 5, 6, 7, 8]) --> true #

    @fact_throws inball(tree, [0.1], 1.0) # n_dim != trees dim
end
end # context

context("KDTrees.inball.tree") do

    data = [0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0;
            0.0 0.0 1.0 1.0 0.0 0.0 1.0 1.0;
            0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0] # 8 node cube

    tree1 = KDTree(data, leafsize = 2, reorder = false)
    tree2 = KDTree([[0.0, 0.0, 0.5] [0.0, 0.0, 0.0]], leafsize = 1, reorder = false)

    idxs = inball(tree1, tree2, 0.6)
    @fact allin(idxs[1], [1, 2]) --> true
    @fact allin(idxs[2], [1]) --> true
end # context

end # facts


text:
using Documenter, HetaSimulator

makedocs(
    sitename = "HetaSimulator Docs",
    authors = "Ivan Borisov, Evgeny Metelkin",
    modules = [HetaSimulator],
    format = Documenter.HTML(
        analytics = "UA-149749027-1",
        #assets = ["assets/favicon.ico"],
        canonical="https://hetalang.github.io/HetaSimulator.jl/dev/"
    ),
    pages = [
        "Home" => "index.md", # readme
        "Basics" => [ # background, discussion
            "basics/overview.md",
            "table-formats/scenario.md",
            "table-formats/measurement.md",
            "table-formats/parameters.md"
        ],
        "Tutorial" => [ # methods
            "tutorial/intro.md",
            "tutorial/sim.md",
            "tutorial/mc.md",
            "tutorial/fit.md",
        ],
        "API" => "api.md"
    ],
)

deploydocs(
    repo = "github.com/hetalang/HetaSimulator.jl.git",
    target = "build",
)


text:
# rpca.jl - Robust PCA solvers using ProximalOperators
#
#   minimize 0.5*||A - S - L||^2_F + lam1*||S||_1 + lam2*||L||_*
#
# See Parikh, Boyd "Proximal Algorithms", §7.2

using LinearAlgebra
using Random
using SparseArrays
using ProximalOperators

Random.seed!(0)

# Define solvers

function rpca_fista(A, lam1, lam2, S, L; tol=1e-3, maxit=50000)
  S_prev = copy(S)
  L_prev = copy(L)
  g = SeparableSum(NormL1(lam1), NuclearNorm(lam2))
  gam = 0.5
  for it = 1:maxit
    # extrapolation step
    S_extr = S + (it-2)/(it+1)*(S - S_prev)
    L_extr = L + (it-2)/(it+1)*(L - L_prev)
    # compute residual
    res = A - S - L
    # compute gradient (forward) step
    y_S = S_extr + gam*res
    y_L = L_extr + gam*res
    # store current iterates
    S_prev .= S
    L_prev .= L
    # compute proximal (backward) step
    prox!((S, L), g, (y_S, y_L), gam)
    # stopping criterion
    fix_point_res = max(norm(S_extr-S, Inf), norm(L_extr-L, Inf))/gam
    rel_fix_point_res = fix_point_res/(1+max(norm(S,Inf), norm(L,Inf)))
    if rel_fix_point_res <= tol
      break
    end
  end
  return S, L
end

# Generate random problem

println("Generating random robust PCA problem")

m, n, r, p, sig = 200, 500, 4, 0.05, 1e-3
L1 = randn(m, r)
L2 = randn(r, n)
L = L1*L2
S = sprand(m, n, p)
V = sig*randn(m, n)
A = L + S + V
lam1 = 0.15*norm(A, Inf)
lam2 = 0.15*opnorm(A)

# Call solvers

println("Calling solvers")

S_fista, L_fista = rpca_fista(A, lam1, lam2, zeros(m, n), zeros(m, n))
println("FISTA")
println("      nnz(S)    = $(count(!isequal(0), S_fista))")
println("      rank(L)   = $(rank(L_fista))")
println("      ||A||     = $(norm(A, Inf))")
println("      ||A-S-L|| = $(norm(A - S_fista - L_fista, Inf))")


text:
using HDF5,JLD,KUparser,Base.Test
include("pxyequal.jl")

# All interfaces give the same result
info("Load data")
@date @load "conll07.tst.jld4"
ft = Flist.acl11eager
pt = ArcEager13
ncpu = 10
ndeps = length(deprel)
s1 = corpus[1]

info("Single sentence")
@date p1=oparse(pt, [s1], ndeps)
info("with features")
@date pxy2=(p2,x2,y2)=oparse(pt, [s1], ndeps, ft)
@test @show isequal(p2,p1)

info("Multiple sentences")
@date p3=oparse(pt, corpus, ndeps)
@test @show isequal(p3[1],p2[1])
info("with features")
@date pxy4=(p4,x4,y4)=oparse(pt, corpus, ndeps, ft)
@test @show isequal(p4, p3)
I2 = 1:size(x2[1],2)
@test @show isequal(x4[1][:,I2], x2[1])
@test @show isequal(y4[1][:,I2], y2[1])

info("Multi-cpu processing")
@date p5=oparse(pt, corpus, ndeps; usepmap=true)
@test @show isequal(p5,p4)
info("with features")
@date pxy6=(p6,x6,y6)=oparse(pt, corpus, ndeps, ft; usepmap=true)
@test @show pxyequal(pxy6,pxy4)


text:
mutable struct BlockMPSState <: QState
    sites::Vector{ITensor}
    dict::Dict{Int, Int}
    llim::Int
    rlim::Int

    BlockMPSState(sites::Vector{ITensor},
                  dict::Dict{Int, Int},
                  llim::Int,
                  rlim::Int) = new(sites, dict, llim, rlim)
end # struct

function BlockMPSState(N::Int, config::Dict{Int, Int})

end


text:
if abspath(PROGRAM_FILE) == @__FILE__
    x::Int64 = 0
    y::Int64 = 0
    @assert(!(x > 2))
    @assert(y < 10)
    @assert((x + 2 * y) == 0)
end


text:
using myJuliaUtils
using Test
using LinearAlgebra
using Statistics

@testset "testing genLPM function - Monte Carlo test" begin
    T = 10000
    K = 400000
    X = [ones(T) randn(T)];
    lammbda = [rand(),randn() / 10.]

    lambbdaMean = 0. .* lammbda
    y = BitArray(undef, T)
    @time for i = 1:K
        y .= genLPM(lammbda, T, X)
        lambbdaMean .+= X \ y

        if mod(i, 100000) == 0
            println("keep")
        end

    end
    lambbdaMean ./= K

    @test_broken isapprox(lambbdaMean, lammbda, atol = 0.01)

end

@testset "testing genMC function" begin
    noStates = round(Int, 2 + rand() * 50)
    S0 = floor(Int, 1 + rand() * noStates)
    PI = eye(noStates)
    T = 100

    @test all(simMC(PI, T, S0 = S0) .== S0)

    noStates = round(Int, 2 + rand() * 20)
    S0 = floor(Int, 1 + rand() * noStates)

    A = randn(noStates, noStates)
    A = abs.(A)
    A = A ./ sum(A;dims = 2)
    B = copy(A)
    A[end,:] = [zeros(1, noStates - 1) 1]
    B[1,:] = [1 zeros(1, noStates - 1)]

    F = simMC(A, 1000, S0 = S0)
    @test F[end] == noStates
    @test F[1] == S0

    F = simMC(B, 1000, S0 = S0)
    @test F[end] == 1
    @test F[1] == S0
end

max_stars_repo_name: zymbuzz/myJuliaUtils

max_issues_repo_name: zymbuzz/myJuliaUtils

max_forks_repo_name: zymbuzz/myJuliaUtils


text:
#=== device apis ===#
@ocl_func(clCreateSubDevices, CL_int,
              (CL_device_id, CL_device_partition_property, CL_uint, Ptr{CL_device_id}, Ptr{CL_uint}))

@ocl_func(clRetainDevice, CL_int, (CL_device_id,))

@ocl_func(clReleaseDevice, CL_int, (CL_device_id,))

#=== memory object apis ===#
@ocl_func(clCreateImage, CL_mem,
              (CL_context, CL_mem_flags, CL_image_format, CL_image_desc, Ptr{Void}, Ptr{CL_int}))

#=== program object apis ===#
@ocl_func(clCreateProgramWithBuiltInKernels, CL_program,
              (CL_context, CL_uint, Ptr{CL_device_id}, Ptr{Cchar}, Ptr{CL_int}))

@ocl_func(clCompileProgram, CL_int,
              (CL_program, CL_uint, Ptr{CL_device_id}, Ptr{CL_device_id}, Ptr{Cchar},
               CL_uint, Ptr{CL_program}, Ptr{Ptr{Char}}, CL_callback, Ptr{Void}))

@ocl_func(clLinkProgram, CL_program,
              (CL_context, CL_uint, Ptr{CL_device_id}, Ptr{Cchar}, CL_uint,
               CL_callback, Ptr{Void}, Ptr{CL_int}))

@ocl_func(clUnloadPlatformCompiler, CL_int, (CL_platform_id,))

#=== kernel object apis ===#
@ocl_func(clGetKernelArgInfo, CL_int,
              (CL_kernel, CL_uint, CL_kernel_arg_info, Csize_t, Ptr{Void}, Ptr{Csize_t}))

#=== enqueued commands apis ===#
@ocl_func(clEnqueueFillBuffer, CL_int,
              (CL_command_queue, CL_mem, Ptr{Void}, Csize_t, Csize_t, Csize_t,
               CL_uint, Ptr{CL_event}, Ptr{CL_event}))

@ocl_func(clEnqueueFillImage, CL_int,
              (CL_command_queue, CL_mem, Ptr{Void}, Ptr{Csize_t}, Ptr{Csize_t},
               CL_uint, Ptr{CL_event}, Ptr{CL_event}))

@ocl_func(clEnqueueMigrateMemObjects, CL_int,
              (CL_command_queue, CL_uint, Ptr{CL_mem}, CL_mem_migration_flags,
               CL_uint, Ptr{CL_event}, Ptr{CL_event}))

@ocl_func(clEnqueueMarkerWithWaitList, CL_int,
              (CL_command_queue, CL_uint, Ptr{CL_event}, Ptr{CL_event}))

@ocl_func(clEnqueueBarrierWithWaitList, CL_int,
              (CL_command_queue, CL_uint, Ptr{CL_event}, Ptr{CL_event}))

#=== extension function access ===#
@ocl_func(clGetExtensionFunctionAddressForPlatform, Ptr{Void},
              (CL_platform_id, Ptr{Cchar}))

#=== opengl interop functions ===#

@ocl_func(clCreateFromGLTexture, CL_mem,
              (CL_context, CL_mem_flags, GL_enum, GL_int, GL_uint, Ptr{CL_int}))

#=== deprecation ===#

# @deprecate clGetExtensionFunctionAddress clGetExtensionFunctionAddressForPlatform

# @deprecate clCreateImage2D clCreateImage
# @deprecate clCreateImage3D clCreateImage

# @deprecate clEnqueueMarker clEnqueueMarkerWithWaitList
# @deprecate clEnqueueBarrier clEnqueueMarkerWithWaitList
# @deprecate clEnqueueWaitForEvents clEnqueueMarkerWithWaitList
# @deprecate clUnloadCompiler Void()

# @deprecate clCreateFromGLTexture2D clCreateFromGLTexture
# @deprecate clCreateFromGLTexture3D clCreateFromGLTexture


text:
using BayesLaplace
using BayesLaplace.Curvature
using BayesLaplace.Data
using Flux

@testset "Jacobian" begin

    # One layer:
    nn = Chain(Dense(2,1))
    la = Laplace(nn)

end

@testset "Hessian" begin

end


text:
defaultOutputFormat(plt::Plot) = "png"

function png(plt::Plot, fn::AbstractString)
    fn = addExtension(fn, "png")
    io = open(fn, "w")
    show(io, MIME("image/png"), plt)
    close(io)
end
png(fn::AbstractString) = png(current(), fn)

function svg(plt::Plot, fn::AbstractString)
    fn = addExtension(fn, "svg")
    io = open(fn, "w")
    show(io, MIME("image/svg+xml"), plt)
    close(io)
end
svg(fn::AbstractString) = svg(current(), fn)

function pdf(plt::Plot, fn::AbstractString)
    fn = addExtension(fn, "pdf")
    io = open(fn, "w")
    show(io, MIME("application/pdf"), plt)
    close(io)
end
pdf(fn::AbstractString) = pdf(current(), fn)

function ps(plt::Plot, fn::AbstractString)
    fn = addExtension(fn, "ps")
    io = open(fn, "w")
    show(io, MIME("application/postscript"), plt)
    close(io)
end
ps(fn::AbstractString) = ps(current(), fn)

function eps(plt::Plot, fn::AbstractString)
    fn = addExtension(fn, "eps")
    io = open(fn, "w")
    show(io, MIME("image/eps"), plt)
    close(io)
end
eps(fn::AbstractString) = eps(current(), fn)

function tex(plt::Plot, fn::AbstractString)
    fn = addExtension(fn, "tex")
    io = open(fn, "w")
    show(io, MIME("application/x-tex"), plt)
    close(io)
end
tex(fn::AbstractString) = tex(current(), fn)

function json(plt::Plot, fn::AbstractString)
    fn = addExtension(fn, "json")
    io = open(fn, "w")
    show(io, MIME("application/vnd.plotly.v1+json"), plt)
    close(io)
end
json(fn::AbstractString) = json(current(), fn)

function html(plt::Plot, fn::AbstractString)
    fn = addExtension(fn, "html")
    io = open(fn, "w")
    show(io, MIME("text/html"), plt)
    close(io)
end
html(fn::AbstractString) = html(current(), fn)

function txt(plt::Plot, fn::AbstractString)
    fn = addExtension(fn, "txt")
    io = open(fn, "w")
    show(io, MIME("text/plain"), plt)
    close(io)
end
txt(fn::AbstractString) = txt(current(), fn)

# ----------------------------------------------------------------

const _savemap = Dict(
    "png" => png,
    "svg" => svg,
    "pdf" => pdf,
    "ps" => ps,
    "eps" => eps,
    "tex" => tex,
    "json" => json,
    "html" => html,
    "tikz" => tex,
    "txt" => txt,
)

const _extension_map = Dict("tikz" => "tex")

function addExtension(fn::AbstractString, ext::AbstractString)
    oldfn, oldext = splitext(fn)
    oldext = chop(oldext, head = 1, tail = 0)
    if get(_extension_map, oldext, oldext) == ext
        return fn
    else
        return string(fn, ".", ext)
    end
end

"""
    savefig([plot,] filename)

Save a Plot (the current plot if `plot` is not passed) to file. The file
type is inferred from the file extension. All backends support png and pdf
file types, some also support svg, ps, eps, html and tex.
"""
function savefig(plt::Plot, fn::AbstractString)
    fn = abspath(expanduser(fn))

    # get the extension
    _, ext = splitext(fn)
    ext = chop(ext, head = 1, tail = 0)
    if isempty(ext)
        ext = defaultOutputFormat(plt)
    end

    # save it
    if haskey(_savemap, ext)
        func = _savemap[ext]
        return func(plt, fn)
    else
        error("Invalid file extension: ", fn)
    end
end
savefig(fn::AbstractString) = savefig(current(), fn)

# ---------------------------------------------------------

"""
    gui([plot])

Display a plot using the backends' gui window
"""
gui(plt::Plot = current()) = display(PlotsDisplay(), plt)

# IJulia only... inline display
function inline(plt::Plot = current())
    isijulia() || error("inline() is IJulia-only")
    Main.IJulia.clear_output(true)
    display(Main.IJulia.InlineDisplay(), plt)
end

function Base.display(::PlotsDisplay, plt::Plot)
    prepare_output(plt)
    _display(plt)
end

_do_plot_show(plt, showval::Bool) = showval && gui(plt)
function _do_plot_show(plt, showval::Symbol)
    showval == :gui && gui(plt)
    showval in (:inline, :ijulia) && inline(plt)
end

# ---------------------------------------------------------

const _best_html_output_type =
    KW(:pyplot => :png, :unicodeplots => :txt, :plotlyjs => :html, :plotly => :html)

# a backup for html... passes to svg or png depending on the html_output_format arg
function _show(io::IO, ::MIME"text/html", plt::Plot)
    output_type = Symbol(plt.attr[:html_output_format])
    if output_type == :auto
        output_type = get(_best_html_output_type, backend_name(plt.backend), :svg)
    end
    if output_type == :png
        # @info("writing png to html output")
        print(
            io,
            "",
        )
    elseif output_type == :svg
        # @info("writing svg to html output")
        show(io, MIME("image/svg+xml"), plt)
    elseif output_type == :txt
        show(io, MIME("text/plain"), plt)
    else
        error("only png or svg allowed. got: $(repr(output_type))")
    end
end

# delegate showable to _show instead
function Base.showable(m::M, plt::P) where {M<:MIME,P<:Plot}
    return hasmethod(_show, Tuple{IO,M,P})
end

function _display(plt::Plot)
    @warn("_display is not defined for this backend.")
end

Base.show(io::IO, m::MIME"text/plain", plt::Plot) = show(io, plt)
# for writing to io streams... first prepare, then callback
for mime in (
    "text/html",
    "image/png",
    "image/eps",
    "image/svg+xml",
    "application/eps",
    "application/pdf",
    "application/postscript",
    "application/x-tex",
    "application/vnd.plotly.v1+json",
)
    @eval function Base.show(io::IO, m::MIME{Symbol($mime)}, plt::Plot)
        if haskey(io, :juno_plotsize)
            showjuno(io, m, plt)
        else
            prepare_output(plt)
            _show(io, m, plt)
        end
        return nothing
    end
end

Base.show(io::IO, m::MIME"application/prs.juno.plotpane+html", plt::Plot) =
    showjuno(io, MIME("text/html"), plt)

"Close all open gui windows of the current backend"
closeall() = closeall(backend())

# function html_output_format(fmt)
#     if fmt == "png"
#         @eval function Base.show(io::IO, ::MIME"text/html", plt::Plot)
#             print(io, "")
#         end
#     elseif fmt == "svg"
#         @eval function Base.show(io::IO, ::MIME"text/html", plt::Plot)
#             show(io, MIME("image/svg+xml"), plt)
#         end
#     else
#         error("only png or svg allowed. got: $fmt")
#     end
# end
#
# html_output_format("svg")

# ---------------------------------------------------------
# Atom PlotPane
# ---------------------------------------------------------
function showjuno(io::IO, m, plt)
    dpi = plt[:dpi]

    jratio = get(io, :juno_dpi_ratio, 1)

    plt[:dpi] = jratio * Plots.DPI

    prepare_output(plt)
    try
        _showjuno(io, m, plt)
    finally
        plt[:dpi] = dpi
    end
end

function _showjuno(io::IO, m::MIME"image/svg+xml", plt)
    if Symbol(plt.attr[:html_output_format]) ≠ :svg
        throw(MethodError(show, (typeof(m), typeof(plt))))
    else
        _show(io, m, plt)
    end
end

Base.showable(::MIME"application/prs.juno.plotpane+html", plt::Plot) = false

_showjuno(io::IO, m, plt) = _show(io, m, plt)


text:
using Test

@testset "TypicalTerm" begin
    include("typical.jl")
end

@testset "linear regression" begin
    include("linear_regression.jl")
end


text:
"""
    ExtensionDifferentialOperator(x̄, method::DifferenceMethod)

Returns a discretized differential operator of `length(x̄)` by `length(x̄) + 2` matrix
whose first and last columns are applied to the ghost nodes just before `x̄[1]` and `x̄[end]` respectively
under no boundary condition using finite difference method specified by `method`.

# Examples
```jldoctest; setup = :(using SimpleDifferentialOperators)
julia> x̄ = 0:5
0:5

julia> ExtensionDifferentialOperator(x̄, BackwardFirstDifference())
4×6 SparseArrays.SparseMatrixCSC{Float64,Int64} with 11 stored entries:
  [1, 1]  =  -1.0
  [1, 2]  =  1.0
  [2, 2]  =  -1.0
  [1, 3]  =  0.0
  [2, 3]  =  1.0
  [3, 3]  =  -1.0
  [2, 4]  =  0.0
  [3, 4]  =  1.0
  [4, 4]  =  -1.0
  [3, 5]  =  0.0
  [4, 5]  =  1.0

julia> ExtensionDifferentialOperator(x̄, ForwardFirstDifference())
4×6 SparseArrays.SparseMatrixCSC{Float64,Int64} with 11 stored entries:
  [1, 2]  =  -1.0
  [2, 2]  =  0.0
  [1, 3]  =  1.0
  [2, 3]  =  -1.0
  [3, 3]  =  0.0
  [2, 4]  =  1.0
  [3, 4]  =  -1.0
  [4, 4]  =  0.0
  [3, 5]  =  1.0
  [4, 5]  =  -1.0
  [4, 6]  =  1.0

julia> ExtensionDifferentialOperator(x̄, CentralSecondDifference())
4×6 SparseArrays.SparseMatrixCSC{Float64,Int64} with 12 stored entries:
  [1, 1]  =  1.0
  [1, 2]  =  -2.0
  [2, 2]  =  1.0
  [1, 3]  =  1.0
  [2, 3]  =  -2.0
  [3, 3]  =  1.0
  [2, 4]  =  1.0
  [3, 4]  =  -2.0
  [4, 4]  =  1.0
  [3, 5]  =  1.0
  [4, 5]  =  -2.0
  [4, 6]  =  1.0
```
"""
function ExtensionDifferentialOperator(x̄::AbstractRange, method::BackwardFirstDifference)
    T = eltype(x̄)
    M = length(x̄) - 2
    Δ = step(x̄)
    return spdiagm(0 => -ones(T, M), 1 => ones(T, M), 2 => zeros(T, M))[1:M,:] / Δ
end

function ExtensionDifferentialOperator(x̄::AbstractRange, method::ForwardFirstDifference)
    T = eltype(x̄)
    M = length(x̄) - 2
    Δ = step(x̄)
    return spdiagm(0 => zeros(T, M), 1 => -ones(T, M), 2 => ones(T, M))[1:M,:] / Δ
end

function ExtensionDifferentialOperator(x̄::AbstractRange, method::CentralSecondDifference)
    T = eltype(x̄)
    M = length(x̄) - 2
    Δ = step(x̄)
    return spdiagm(0 => ones(T, M), 1 => -2*ones(T, M), 2 => ones(T, M))[1:M,:] / Δ^2
end

function ExtensionDifferentialOperator(x̄::AbstractArray, method::BackwardFirstDifference)
    T = eltype(x̄)
    M = length(x̄) - 2
    d = diff(x̄)
    Δ₋⁻¹ = 1 ./ d[1:end-1]
    return spdiagm(0 => -Δ₋⁻¹, 1 => Δ₋⁻¹, 2 => zeros(T, M))[1:M,:]
end

function ExtensionDifferentialOperator(x̄::AbstractArray, method::ForwardFirstDifference)
    T = eltype(x̄)
    M = length(x̄) - 2
    d = diff(x̄)
    Δ₊⁻¹ = 1 ./ d[2:end]
    return spdiagm(0 => zeros(T, M), 1 => -Δ₊⁻¹, 2 => Δ₊⁻¹)[1:M,:]
end

function ExtensionDifferentialOperator(x̄::AbstractArray, method::CentralSecondDifference)
    T = eltype(x̄)
    M = length(x̄) - 2
    d = diff(x̄)
    Δ₋⁻¹ = 1 ./ d[1:end-1] # 1 ./ Δ₋
    Δ₊⁻¹ = 1 ./ d[2:end] # 1 ./ Δ₊
    Δ⁻¹ = 1 ./ (d[1:end-1] + d[2:end]) # 1 ./ (Δ₋ + Δ₊)
    return 2*spdiagm(0 => Δ⁻¹.*Δ₋⁻¹, 1 => -Δ₋⁻¹ .* Δ₊⁻¹, 2 => Δ⁻¹.*Δ₊⁻¹)[1:M,:]
end


text:
using Test

function example_test()
    # Test stuff here.
    ret = (1==1)
    return ret
end

function tests()

    @testset "Tests" begin
        @testset "Example Test" begin
            @test example_test() == true
        end
    end
end

tests()


text:
load_met_readable_ids() = ldat(HumanGEM,
    "met_readable_ids", ".jls";
    verbose = false
)


text:
include("opt_berpaths_mc_serial.jl")

n_sim =131072
tic()
optval = cal_optvalue(32, 5, 0.3, 10, 0.06, 1, 1, n_sim)
@printf("MC Est. : %f, Runtime: %s\n", optval, toc())


text:
#==============================================================================#
# module_jl_lambda_eval.jl
#
# Default AWS Lambda function for Julia.
#
# See http://docs.aws.amazon.com/lambda/latest/dg/API_Reference.html
#
# Copyright OC Technology Pty Ltd 2014 - All rights reserved
#==============================================================================#

__precompile__()

module module_jl_lambda_eval

lambda_function(f::Function) = Base.invokelatest(f)

lambda_function(e::Expr) = eval(Main, e)

function lambda_function(v::Vector)
    for e in v[1:end-1]
        eval(Main, e)
    end
    eval(Main, v[end])
end

lambda_function_with_event(event::String) = include_string(event)

precompile(lambda_function, (Function,))
precompile(lambda_function, (Expr,))
precompile(lambda_function, (Vector{Expr},))
precompile(lambda_function_with_event, (String,))

end

#==============================================================================#
# End of file.
#==============================================================================#


text:
struct MapLeaves{F}
    f::F
end
MapLeaves(f::MapLeaves) = f

struct MapLeavesIndex{F,Args}
    f::F
    args::Args
end

MapLeavesIndex(f::MapLeaves{F}, args::Args) where {F,Args} =
    MapLeavesIndex{F,Args}(f.f, args)
MapLeaves(f::MapLeavesIndex{F}) where {F} = MapLeaves{F}(f.f)

(f::MapLeavesIndex)(i::Integer) = MapLeaves(f)(getindex.(f.args, i)...)
(f::MapLeaves)(x...) = f.f(x...)
@generated function (f::MapLeaves)(x::Vararg{NTuple{N,T},K}) where {N,T,K}
    getters = map(Base.OneTo(N)) do i
        exprs = [:(x[$j][$i]) for j in Base.OneTo(K)]
        return :(f($(exprs...)))
    end
    return :(($(getters...),))
end
(f::MapLeaves)(x::AbstractArray...) = map(f, x...)

mapleaves(f::F, x...) where {F} = MapLeaves(f)(x...)

#####
##### Record Type for keeping track of what we're measuring
#####

"""
    Record{name}(data::T) where {T <: Union{Vector, NTuple}}

Create a named `Record` around `data`.
The elements of `data` may themselves be `Record`s, resulting in a hierarchical data structure.
"""
struct Record{name,T<:Union{Vector,NTuple}}
    data::T
end
Record{name}(data::T) where {name,T} = Record{name,T}(data)
Record{name}(data::Tuple{Tuple}) where {name} = Record{name}(data[1])

name(::Record{N}) where {N} = N

# Check if there are records below this record.
hassubrecord(R::Record) = eltype(R.data) <: Record

Base.getindex(R::Record, i...) = R.data[i...]
Base.length(R::Record) = length(R.data)

Base.iterate(R::Record) = iterate(R.data)
Base.iterate(R::Record, s) = iterate(R.data, s)

denest(x) = x
denest(x::Tuple{Tuple}) = denest(first(x))

# Recurse down the record stack, but apply the name to whatever is returned.
"""
    mapleaves(f, record::Record) -> Record

Apply `f` to each leaf element of `record`.
This will recursively descend through hierarchies of `Records` and only apply `f` to scalars.

The returned result will have the same hierarchical structure as `record`
"""
@generated function (f::MapLeaves)(R::Vararg{Record{name},N}) where {name,N}
    getters = [:(R[$i].data) for i in Base.OneTo(N)]
    return :(Record{name}(f($(getters...))))
end

# Applying a difference to two subsequent records
Base.:-(a::R, b::R) where {R<:Record} = mapleaves(-, a, b)
Base.:+(a::R, b::R) where {R<:Record} = mapleaves(+, a, b)
mapleaves(f) = (x...) -> mapleaves(f, x...)

# Aggregate across records
struct Aggregate{F}
    f::F
end

(f::Aggregate)(x) = x
(f::Aggregate)(x, y) = f.f.(x, y)
(f::Aggregate)(x::Record) = f(x.data)
const AggregateRecurseTypes =
    Union{AbstractVector,NTuple{<:Any,<:Record},NTuple{<:Any,<:NTuple}}
(f::Aggregate)(x::AggregateRecurseTypes) = reduce(f, map(f, x))

"""
    aggregate(record::Record)
    aggregate(f, record::Record)

Reduce over all the leaf (terminal) elements of `record`, applying `f` as the reduction function.
If `f` is not supplied, it will defult to `(x,y) -> x .+ y`.
"""
aggregate(f::F, x) where {F} = Aggregate{F}(f)(x)
aggregate(x) = aggregate(+, x)
#broadcast_add(x, y) = x .+ y

#####
##### Pretty Printing
#####

titleize(x) = titlecase(replace(String(x), "_" => " "))
function Base.show(io::IO, R::Record, pre = "")
    println(io, pre, titleize(name(R)), " Record with $(length(R)) entries:")
    post = pre * "   "
    if hassubrecord(R)
        show(io, first(R.data), post)
    else
        print(io, post, eltype(R.data))
    end

    # Only show contents from the top level
    if isempty(pre)
        println(io)
        println(io)
        showcontents(io, R)
    end
end

function showcontents(io::IO, R::Record, pre = "")
    post = pre * "   "
    print(io, pre, titleize(name(R)), ": ")
    if hassubrecord(R)
        print(io, "\n")
        showcontents.(Ref(io), R.data, post)
    else
        println(io, R.data)
    end
end


text:
using DrillHoles
using DataFrames
using Test

@testset "DrillHoles.jl" begin
    # dummy data
    collar = Collar(file=DataFrame(HOLEID=1:2, X=1:2, Y=1:2, Z=1:2))
    assays = Interval(file=DataFrame(HOLEID=[1,1,2], FROM=[1,3.5,0], TO=[3.5,8,7], A=[1,2,3]))
    lithos = Interval(file=DataFrame(HOLEID=[1,2,2], FROM=[0,0,4.4], TO=[8,4.4,8], L=["A","B","C"]))
    survey = Survey(file=DataFrame(HOLEID=[1,1,2,2], AT=[0,5,0,5], AZM=[0,1,20,21], DIP=[89,88,77,76]),
                    method=:mincurv, convention=:negativedownwards)

    # drill hole desurvey tests
    dh  = drillhole(collar, survey, [assays, lithos])
    tab = dh.table

    @test size(tab, 1) == 6
    @test size(tab, 2) == 9
    @test all(isapprox.(tab[!,:HOLEID], [1.0, 1.0, 1.0, 2.0, 2.0, 2.0]))
    @test all(isapprox.(tab[!,:FROM],   [0.0, 1.0, 3.5, 0.0, 4.4, 7.0]))
    @test all(isapprox.(tab[!,:TO],     [1.0, 3.5, 8.0, 4.4, 7.0, 8.0]))
    @test all(isequal.(tab[!,:L],       ["A", "A", "A", "B", "C", "C"]))
    @test all(isequal.(tab[!,:A],       [missing, 1, 2, 3, 3, missing]))
    @test all(round.(tab[!,:X], digits=1) .≈ [1.0, 1.0, 1.0, 2.2, 2.5, 2.6])
    @test all(round.(tab[!,:Y], digits=1) .≈ [1.0, 1.0, 1.2, 2.5, 3.3, 3.7])
    @test all(round.(tab[!,:Z], digits=1) .≈ [1.5, 3.2, 6.7, 4.1, 7.5, 9.3])

    # drill hole :equalcomp compositing tests
    comp1 = composite(dh; zone=:L, mode=:equalcomp)
    tab   = comp1.table

    @test size(tab, 1) == 16
    @test size(tab, 2) == 9
    @test all(isapprox.(tab[!,:HOLEID], [1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2]))
    @test all(isapprox.(tab[!,:FROM],   [0,1,2,3,4,5,6,7,0,1,2,3,4.4,5.4,6.4,7.4]))
    @test all(isapprox.(tab[!,:TO],     [1,2,3,4,5,6,7,8,1,2,3,4,5.4,6.4,7.4,8]))
    @test all(isequal.(tab[!,:L],       ["A","A","A","A","A","A","A","A","B","B","B","B","C","C","C","C"]))
    @test all(round.(tab[!,:X], digits=1) .≈ [1,1,1,1,1,1,1,1,2,2.1,2.2,2.3,2.4,2.5,2.6,2.6])
    @test all(round.(tab[!,:Y], digits=1) .≈ [1,1,1.1,1.1,1.1,1.1,1.2,1.2,2.1,2.3,2.5,2.8,3.1,3.3,3.5,3.7])
    @test all(round.(tab[!,:Z], digits=1) .≈ [1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,2.5,3.5,4.4,5.4,6.8,7.7,8.7,9.5])
    @test all(isequal.(round.(tab[!,:A], digits=1), [missing,1,1,1.5,2,2,2,2,3,3,3,3,3,3,3,missing]))

    # drill hole :nodiscard compositing tests
    comp2 = composite(dh; zone=:L, mode=:nodiscard)
    tab   = comp2.table

    @test size(tab, 1) == 16
    @test size(tab, 2) == 9
    @test all(isapprox.(tab[!,:HOLEID], [1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2]))
    @test all(isapprox.(tab[!,:FROM],   [0,1,2,3,4,5,6,7,0,1.1,2.2,3.3,4.4,5.3,6.2,7.1]))
    @test all(isapprox.(tab[!,:TO],     [1,2,3,4,5,6,7,8,1.1,2.2,3.3,4.4,5.3,6.2,7.1,8]))
    @test all(isequal.(tab[!,:L],       ["A","A","A","A","A","A","A","A","B","B","B","B","C","C","C","C"]))
    @test all(round.(tab[!,:X], digits=1) .≈ [1,1,1,1,1,1,1,1,2,2.1,2.2,2.3,2.4,2.5,2.6,2.6])
    @test all(round.(tab[!,:Y], digits=1) .≈ [1,1,1.1,1.1,1.1,1.1,1.2,1.2,2.1,2.4,2.6,2.8,3.1,3.3,3.5,3.7])
    @test all(round.(tab[!,:Z], digits=1) .≈ [1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,2.5,3.6,4.7,5.7,6.7,7.6,8.5,9.3])
    @test all(isequal.(round.(tab[!,:A], digits=1), [missing,1,1,1.5,2,2,2,2,3,3,3,3,3,3,3,missing]))

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

    # tests using different survey method and dip convention
    survey = Survey(file=DataFrame(HOLEID=[1,1,2,2], AT=[0,5,0,5], AZM=[0,1,20,21], DIP=[89,88,77,76]),
                    method=:tangential, convention=:auto)

    # drill hole desurvey tests
    dh  = drillhole(collar, survey, [assays, lithos])
    tab = dh.table

    @test dh.pars.invdip == true
    @test size(tab, 1) == 6
    @test size(tab, 2) == 9
    @test all(isapprox.(tab[!,:HOLEID], [1.0, 1.0, 1.0, 2.0, 2.0, 2.0]))
    @test all(isapprox.(tab[!,:FROM],   [0.0, 1.0, 3.5, 0.0, 4.4, 7.0]))
    @test all(isapprox.(tab[!,:TO],     [1.0, 3.5, 8.0, 4.4, 7.0, 8.0]))
    @test all(isequal.(tab[!,:L],       ["A", "A", "A", "B", "C", "C"]))
    @test all(isequal.(tab[!,:A],       [missing, 1, 2, 3, 3, missing]))
    @test all(round.(tab[!,:X], digits=1) .≈ [1.0, 1.0, 1.0, 2.2, 2.4, 2.6])
    @test all(round.(tab[!,:Y], digits=1) .≈ [1.0, 1.0, 1.1, 2.5, 3.2, 3.6])
    @test all(round.(tab[!,:Z], digits=1) .≈ [0.5,-1.2,-4.7,-0.1,-3.6,-5.3])

    # drill hole :equalcomp compositing tests
    comp1 = composite(dh; zone=:L, mode=:equalcomp)
    tab   = comp1.table

    @test size(tab, 1) == 16
    @test size(tab, 2) == 9
    @test all(isapprox.(tab[!,:HOLEID], [1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2]))
    @test all(isapprox.(tab[!,:FROM],   [0,1,2,3,4,5,6,7,0,1,2,3,4.4,5.4,6.4,7.4]))
    @test all(isapprox.(tab[!,:TO],     [1,2,3,4,5,6,7,8,1,2,3,4,5.4,6.4,7.4,8]))
    @test all(isequal.(tab[!,:L],       ["A","A","A","A","A","A","A","A","B","B","B","B","C","C","C","C"]))
    @test all(round.(tab[!,:X], digits=1) .≈ [1,1,1,1,1,1,1,1,2,2.1,2.2,2.3,2.4,2.5,2.5,2.6])
    @test all(round.(tab[!,:Y], digits=1) .≈ [1,1,1,1.1,1.1,1.1,1.1,1.2,2.1,2.3,2.5,2.7,3,3.3,3.5,3.7])
    @test all(round.(tab[!,:Z], digits=1) .≈ [0.5,-0.5,-1.5,-2.5,-3.5,-4.5,-5.5,-6.5,1.5,0.5,-0.4,-1.4,-2.8,-3.7,-4.7,-5.5])
    @test all(isequal.(round.(tab[!,:A], digits=1), [missing,1,1,1.5,2,2,2,2,3,3,3,3,3,3,3,missing]))

    # drill hole :nodiscard compositing tests
    comp2 = composite(dh; zone=:L, mode=:nodiscard)
    tab   = comp2.table

    @test size(tab, 1) == 16
    @test size(tab, 2) == 9
    @test all(isapprox.(tab[!,:HOLEID], [1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2]))
    @test all(isapprox.(tab[!,:FROM],   [0,1,2,3,4,5,6,7,0,1.1,2.2,3.3,4.4,5.3,6.2,7.1]))
    @test all(isapprox.(tab[!,:TO],     [1,2,3,4,5,6,7,8,1.1,2.2,3.3,4.4,5.3,6.2,7.1,8]))
    @test all(isequal.(tab[!,:L],       ["A","A","A","A","A","A","A","A","B","B","B","B","C","C","C","C"]))
    @test all(round.(tab[!,:X], digits=1) .≈ [1,1,1,1,1,1,1,1,2,2.1,2.2,2.3,2.4,2.4,2.5,2.6])
    @test all(round.(tab[!,:Y], digits=1) .≈ [1,1,1,1.1,1.1,1.1,1.1,1.2,2.1,2.3,2.6,2.8,3,3.2,3.4,3.6])
    @test all(round.(tab[!,:Z], digits=1) .≈ [0.5,-0.5,-1.5,-2.5,-3.5,-4.5,-5.5,-6.5,1.5,0.4,-0.7,-1.8,-2.7,-3.6,-4.5,-5.3])
    @test all(isequal.(round.(tab[!,:A], digits=1), [missing,1,1,1.5,2,2,2,2,3,3,3,3,3,3,3,missing]))
end


text:
@recipe function f(A::AbstractInterpolation; plotdensity = 10_000, denseplot = true)
  t = sort(A.t)
  start = t[1]; stop = t[end]
  if denseplot
    plott = collect(range(start,stop=stop,length=plotdensity))
  else
    plott = t
  end
  output = A.(plott)
  if !(eltype(output) <: AbstractVector)
    output = [[x] for x in output]
  end
  DiffEqArray(output, plott)
end

@recipe function f(A::GPInterpolation; plotdensity = 10_000, denseplot = true)
  t = sort(A.t)
  start = t[1]; stop = t[end]
  if denseplot
    plott = collect(range(start,stop=stop,length=plotdensity))
  else
    plott = t
  end
  plott, A(plott)
end


text:
n, p = 500, 5
((X, y, θ), (X1, y1, θ1)) = generate_continuous(n, p; seed=525)

# adding some outliers (both positive and negative)
y1a = outlify(y1, 0.1)

@testset "QuantileReg" begin
    δ = 0.5 # effectively LAD regression
    λ = 1.0
    rr = QuantileRegression(δ, lambda=λ, penalize_intercept=true,
                               scale_penalty_with_samples = false)
    J = objective(rr, X, y1a)
    o = RobustLoss(Quantile(δ)) + λ * L2Penalty()
    @test J(θ1) ≈ o(y1a, X1*θ1, θ1)
    ls = LinearRegression()
    θ_ls    = fit(ls, X, y1a)
    θ_lbfgs = fit(rr, X, y1a, solver=LBFGS())
    θ_iwls  = fit(rr, X, y1a, solver=IWLSCG())
    @test isapprox(J(θ1),      412.20773, rtol=1e-5)
    @test isapprox(J(θ_ls),    508.02443, rtol=1e-5)  # LS is crap bc outliers
    @test isapprox(J(θ_lbfgs), 411.98228, rtol=1e-5)
    @test isapprox(J(θ_iwls),  411.98,    rtol=1e-4)

    # NOTE: newton and newton-cg not available because ϕ = 0 identically
    # will throw an error if called.
    @test_throws ErrorException fit(rr, X, y1, solver=Newton())
    @test_throws ErrorException fit(rr, X, y1, solver=NewtonCG())

    # don't penalize intercept
    rr = QuantileRegression(δ, lambda=λ, scale_penalty_with_samples = false)
    J = objective(rr, X, y1a)
    ls = LinearRegression()
    θ_ls    = fit(ls, X, y1a)
    θ_lbfgs = fit(rr, X, y1a, solver=LBFGS())
    θ_iwls  = fit(rr, X, y1a, solver=IWLSCG())
    @test isapprox(J(θ1),      412.18594, rtol=1e-5)
    @test isapprox(J(θ_ls),    508.00993, rtol=1e-5)  # note that LS is crap due to outliers
    @test isapprox(J(θ_lbfgs), 411.95990, rtol=1e-5)
    @test isapprox(J(θ_iwls),  411.98,    rtol=1e-4)

    if DO_COMPARISONS
        # Compare with R's QuantReg package
        # NOTE: QuantReg doesn't allow for penalties so re-fitting with λ=0
        rr = QuantileRegression(δ, lambda=0)
        J  = objective(rr, X, y1a)
        θ_ls     = fit(LinearRegression(), X, y1a)
        θ_lbfgs  = fit(rr, X, y1a, solver=LBFGS())
        θ_iwls   = fit(rr, X, y1a, solver=IWLSCG())
        θ_qr_br  = rcopy(QUANTREG.rq_fit_br(X1, y1a))[:coefficients]
        θ_qr_fnb = rcopy(QUANTREG.rq_fit_fnb(X1, y1a))[:coefficients]
        # NOTE: we take θ_qr_br as reference point
        @test isapprox(J(θ_ls), 505.45286,  rtol=1e-5)
        @test J(θ_qr_br) ≈      409.570777 # <- ref value
        # Their IP algorithm essentially gives the same answer
        @test (J(θ_qr_fnb) - J(θ_qr_br)) ≤ 1e-10
        # Our algorithms are close enough
        @test isapprox(J(θ_lbfgs), 409.57154, rtol=1e-5)
        @test isapprox(J(θ_iwls),  409.59,    rtol=1e-4)
    end
end

###########################
## With Sparsity penalty ##
###########################

n, p = 500, 100
((X, y, θ), (X1, y1, θ1)) = generate_continuous(n, p;  seed=51112, sparse=0.1)
# pepper with outliers
y1a  = outlify(y1, 0.1)

@testset "LAD+L1" begin
    λ = 5.0
    γ = 10.0
    rr = LADRegression(λ, γ; penalize_intercept=true, scale_penalty_with_samples = false)
    J  = objective(rr, X, y1a)
    θ_ls    = X1 \ y1a
    θ_fista = fit(rr, X, y1a, solver=FISTA())
    θ_ista  = fit(rr, X, y1a, solver=ISTA())
    @test isapprox(J(θ_ls),    1058.6737, rtol=1e-5)
    @test isapprox(J(θ_fista), 454.23322, rtol=1e-5)
    @test isapprox(J(θ_ista),  454.27774, rtol=1e-5)
    @test nnz(θ_fista) == 43
    @test nnz(θ_ista)  == 43

    if DO_COMPARISONS
        # Compare with R's QuantReg package
        # NOTE: QuantReg doesn't apply the penalty on the intercept
        rr = LADRegression(5.0; penalty=:l1, scale_penalty_with_samples = false)
        J  = objective(rr, X, y1a)
        θ_ls       = X1 \ y1a
        θ_fista    = fit(rr, X, y1a, solver=FISTA())
        θ_ista     = fit(rr, X, y1a, solver=ISTA())
        θ_qr_lasso = rcopy(QUANTREG.rq_fit_lasso(X1, y1a))[:coefficients]
        @test isapprox(J(θ_ls),       888.3748, rtol=1e-5)
        @test isapprox(J(θ_qr_lasso), 425.5264, rtol=1e-5)
        # Our algorithms are close enough
        @test isapprox(J(θ_fista),    425.0526, rtol=1e-5)
        @test isapprox(J(θ_ista),     425.4113, rtol=1e-5)
        # in this case we do a fair bit better
        @test nnz(θ_qr_lasso) == 101
        @test nnz(θ_fista)    == 88
        @test nnz(θ_ista)     == 82
        # in this case fista is best
        @test J(θ_fista) < J(θ_ista) < J(θ_qr_lasso)
    end
end


text:
# Logging helping functions based on ConsoleLogger
# It's goal is to hide the lines where the log is coming from.

function ns_mfmt(level, _module, group, id, file, line)
    color = default_logcolor(level)
    prefix = (level == Logging.Warn ? "Warning" : string(level))*':'
    suffix = ""
    return color, prefix, suffix
end

nosource_debuglogger() = ConsoleLogger(stderr, Logging.Debug, meta_formatter=ns_mfmt)
nosource_infologger() = ConsoleLogger(stderr, Logging.Info, meta_formatter=ns_mfmt)
sourced_debuglogger() = ConsoleLogger(stderr, Logging.Debug)
sourced_infologger() = ConsoleLogger(stderr, Logging.Info)

# custom LogLevels

const IterInf = Logging.Info-5
const IterLow1 = Logging.Info-10
const IterLow2 = Logging.Info-15
const Compute1 = Logging.Info-20
const Compute2 = Logging.Info-25

# custom loggers for custom Loglevels

nosource_IterInflog() = ConsoleLogger(stderr, IterInf, meta_formatter=ns_mfmt)
nosource_IterLow1log() = ConsoleLogger(stderr, IterLow1, meta_formatter=ns_mfmt)
nosource_IterLow2log() = ConsoleLogger(stderr, IterLow2, meta_formatter=ns_mfmt)
nosource_Compute1log() = ConsoleLogger(stderr, Compute1, meta_formatter=ns_mfmt)
nosource_Compute2log() = ConsoleLogger(stderr, Compute2, meta_formatter=ns_mfmt)
nosource_Errorlog() = ConsoleLogger(stderr, Logging.Error, meta_formatter=ns_mfmt)


text:
## Top-k pooling

function topk_index(y::AbstractVector, k::Integer)
    v = nlargest(k, y)
    return collect(1:length(y))[y .>= v[end]]
end

topk_index(y::Adjoint, k::Integer) = topk_index(y', k)

## Get feature with defaults

get_feature(::Nothing, i) = nothing
get_feature(A::Fill{T,2,Axes}, i::Integer) where {T,Axes} = view(A, :, 1)
get_feature(A::AbstractMatrix, i::Integer) = view(A, :, i)

"""
    bypass_graph(nf_func, ef_func, gf_func)

Bypassing graph in FeaturedGraph and let other layer process (node, edge and global)features only.
"""
function bypass_graph(nf_func=identity, ef_func=identity, gf_func=identity)
    return function (fg::FeaturedGraph)
        FeaturedGraph(graph(fg), nf_func(node_feature(fg)), ef_func(edge_feature(fg)),
                      gf_func(global_feature(fg)))
    end
end

function add_self_loop!(adj::AbstractVector{T}, n::Int=length(adj)) where {T<:AbstractVector}
    for i = 1:n
        i in adj[i] || push!(adj[i], i)
    end
    adj
end


text:
function constraint_voltage_product_converter(pm::_PM.AbstractWRModel, wr, wi, w_fr, w_to)
    InfrastructureModels.relaxation_complex_product(pm.model, w_fr, w_to, wr, wi)
end

function constraint_voltage_product_converter(pm::_PM.AbstractWRConicModel, wr, wi, w_fr, w_to)
    InfrastructureModels.relaxation_complex_product_conic(pm.model, w_fr, w_to, wr, wi)
end
"""
Links converter power & current

```
pconv_ac[i]^2 + pconv_dc[i]^2 <= wc[i] * iconv_ac_sq[i]
pconv_ac[i]^2 + pconv_dc[i]^2 <= (Umax)^2 * (iconv_ac[i])^2
```
"""
function constraint_converter_current(pm::_PM.AbstractWRModel, n::Int, i::Int, Umax, Imax)
    wc = _PM.var(pm, n,  :wc_ac, i)
    pconv_ac = _PM.var(pm, n,  :pconv_ac, i)
    qconv_ac = _PM.var(pm, n,  :qconv_ac, i)
    iconv = _PM.var(pm, n,  :iconv_ac, i)
    iconv_sq = _PM.var(pm, n,  :iconv_ac_sq, i)

    JuMP.@constraint(pm.model, pconv_ac^2 + qconv_ac^2 <= wc * iconv_sq)
    JuMP.@constraint(pm.model, pconv_ac^2 + qconv_ac^2 <= (Umax)^2 * iconv^2)
    JuMP.@constraint(pm.model, iconv^2 <= iconv_sq)
    JuMP.@constraint(pm.model, iconv_sq <= iconv*Imax)
end

function constraint_converter_current(pm::_PM.AbstractWRConicModel, n::Int,  i::Int, Umax, Imax)
    wc = _PM.var(pm, n,  :wc_ac, i)
    pconv_ac = _PM.var(pm, n,  :pconv_ac, i)
    qconv_ac = _PM.var(pm, n,  :qconv_ac, i)
    iconv = _PM.var(pm, n,  :iconv_ac, i)
    iconv_sq = _PM.var(pm, n,  :iconv_ac_sq, i)

    JuMP.@constraint(pm.model, [wc/sqrt(2), iconv_sq/sqrt(2), pconv_ac, qconv_ac] in JuMP.RotatedSecondOrderCone())
    JuMP.@constraint(pm.model, [Umax * iconv/sqrt(2), Umax * iconv/sqrt(2), pconv_ac, qconv_ac] in JuMP.RotatedSecondOrderCone())
    JuMP.@constraint(pm.model, iconv_sq <= iconv*Imax)
end

################ TNEP constraints ##################
function constraint_voltage_product_converter_ne(pm::_PM.AbstractWRConicModel, wr, wi, w_fr, w_to, z)
    JuMP.@constraint(pm.model, wr >= z*JuMP.lower_bound(wr))
    JuMP.@constraint(pm.model, wr <= z*JuMP.upper_bound(wr))
    JuMP.@constraint(pm.model, wi >= z*JuMP.lower_bound(wi))
    JuMP.@constraint(pm.model, wi <= z*JuMP.upper_bound(wi))
    relaxation_complex_product_conic_on_off(pm.model, w_fr, w_to, wr, wi, z)
end

function constraint_voltage_product_converter_ne(pm::_PM.AbstractWRModel, wr, wi, w_fr, w_to, z)
    JuMP.@constraint(pm.model, wr >= z*JuMP.lower_bound(wr))
    JuMP.@constraint(pm.model, wr <= z*JuMP.upper_bound(wr))
    JuMP.@constraint(pm.model, wi >= z*JuMP.lower_bound(wi))
    JuMP.@constraint(pm.model, wi <= z*JuMP.upper_bound(wi))
    _IM.relaxation_complex_product_on_off(pm.model, w_fr, w_to, wr, wi, z)
end

"""
Links converter power & current

```
pconv_ac[i]^2 + pconv_dc[i]^2 <= wc[i] * iconv_ac_sq[i]
pconv_ac[i]^2 + pconv_dc[i]^2 <= (Umax)^2 * (iconv_ac[i])^2
```
"""
function constraint_converter_current_ne(pm::_PM.AbstractWRConicModel, n::Int, i::Int, Umax, Imax)
    wc = _PM.var(pm, n, :wc_ac_ne, i)
    pconv_ac = _PM.var(pm, n, :pconv_ac_ne, i)
    qconv_ac = _PM.var(pm, n, :qconv_ac_ne, i)
    iconv = _PM.var(pm, n, :iconv_ac_ne, i)
    iconv_sq = _PM.var(pm, n, :iconv_ac_sq_ne, i)
    JuMP.@constraint(pm.model, [wc/sqrt(2), iconv_sq/sqrt(2), pconv_ac, qconv_ac] in JuMP.RotatedSecondOrderCone())
    JuMP.@constraint(pm.model, [Umax * iconv/sqrt(2), Umax * iconv/sqrt(2), pconv_ac, qconv_ac] in JuMP.RotatedSecondOrderCone())
    JuMP.@constraint(pm.model, iconv_sq <= iconv*Imax)
end

function constraint_converter_current_ne(pm::_PM.AbstractWRModel, n::Int, i::Int, Umax, Imax)
    wc = _PM.var(pm, n, :wc_ac_ne, i)
    pconv_ac = _PM.var(pm, n, :pconv_ac_ne, i)
    qconv_ac = _PM.var(pm, n, :qconv_ac_ne, i)
    iconv = _PM.var(pm, n, :iconv_ac_ne, i)
    iconv_sq = _PM.var(pm, n, :iconv_ac_sq_ne, i)

    JuMP.@constraint(pm.model, pconv_ac^2 + qconv_ac^2 <= wc * iconv_sq)
    JuMP.@constraint(pm.model, pconv_ac^2 + qconv_ac^2 <= (Umax)^2 * iconv^2)
    JuMP.@constraint(pm.model,  iconv^2 <= iconv_sq)
    JuMP.@constraint(pm.model, iconv_sq <= iconv*Imax)
end


text:
using BlackBoxOptimizationBenchmarking, Optim
import BlackBoxOptimizationBenchmarking: minimizer, minimum, optimize
import Base.string

const BBOB = BlackBoxOptimizationBenchmarking

box(D) = fill((-5.5, 5.5),D)
pinit(D) = 10*rand(D).-5

optimize(opt::Optim.SAMIN,f,D,run_length) =
    Optim.optimize(f, fill(-5.5,D), fill(5.5,D), pinit(D), opt, Optim.Options(f_calls_limit=run_length,g_tol=1e-120,iterations=run_length))

run_length = 30_000
dimensions = 3
Ntrials = 20
Δf = 1e-6

# test on Sphere
f = BlackBoxOptimizationBenchmarking.F1
success_rate, x_dist, f_dist, runtime = BBOB.benchmark(Optim.SAMIN(), f, run_length, Ntrials, dimensions, Δf)

#test on all functions
perf = [(f=>BBOB.benchmark(Optim.SAMIN(), f, run_length, Ntrials, dimensions, Δf)[1][1]) for f in enumerate(BBOBFunction)]


text:
"""
Applies a function to each timeseries of a datacube present in a mask.
# Example:
```
julia> datacube = rand(1:10.0, 10,10,10)
julia> mask = rand(0:1, 10, 10)
julia> long_apply(sum, datacube)
```
"""
function long_apply(f, datacube::Array, mask = ones(size(datacube)[1], size(datacube)[2]))
    if ndims(f(datacube[1, 1, :])) ==0
        new_matrix = Array{Union{Missing, Float16}, 3}(undef, size(datacube, 1), size(datacube, 2))
        for i in 1:size(datacube)[1]
            for j in 1:size(datacube)[2]
                new_matrix[i, j] = f(datacube[i, j, :])
            end
        end
        return new_matrix
    else
        new_matrix = Array{Union{Missing, Float16}}(undef, size(datacube, 1), size(datacube, 2),size(datacube)[3])
        @showprogress for i in 1:size(datacube)[1]
            for j in 1:size(datacube)[2]
                if mask[i, j] == 1
                    new_matrix[i, j, :] = f(datacube[i, j, :])
                else
                    new_matrix[i, j, :] = datacube[i, j, :]
                end
            end
        end
        return new_matrix
    end
end
"""
Applies a function to each cross section of a datacube present.
# Example:
```
julia> datacube = rand(1:10.0, 10,10,10)
julia> mask = rand(0:1, 10, 10)
julia> long_apply(sum, datacube)
```
"""
function cross_apply(f, datacube, mask = ones(size(datacube)[1], size(datacube)[2]))
    if ndims(f(datacube[:, :, 1])) == 2
        new_matrix = Array{Union{Missing, Float16}}(undef, size(datacube, 1), size(datacube, 2),size(datacube)[3])
        for i in 1:size(datacube)[3]
            new_matrix[:, :, i] = f(datacube[:, :, i])
        end
        return new_matrix
    else
        array = []
        for i in 1:size(datacube)[3]
            push!(array, f(datacube[:, :, i]))
        end
        return array
    end
end


text:
#== HyperParameter Type ==#

mutable struct HyperParameter{T<:Real}
    value::T
    interval::Interval{T}
    fixed::Bool
    opt::Optimizer
    # function HyperParameter{T}(x::T,I::Interval{T};fixed::Bool=false,opt::Optimizer=Momentum(η=0.01)) where {T<:Real}
    function HyperParameter{T}(x::T,I::Interval{T};fixed::Bool=false,opt::Optimizer=Adam(α=0.01)) where {T<:Real}
    # function HyperParameter{T}(x::T,I::Interval{T};fixed::Bool=false,opt::Optimizer=VanillaGradDescent(η=0.001)) where {T<:Real}
        checkvalue(I, x) || error("Value $(x) must be in range " * string(I))
        new{T}(x, I, fixed, opt)
    end
end
HyperParameter(x::T, I::Interval{T} = interval(T); fixed::Bool = false, opt::Optimizer = Adam(α=0.01)) where {T<:Real} = HyperParameter{T}(x, I, fixed, opt)

eltype(::HyperParameter{T}) where {T} = T

@inline getvalue(θ::HyperParameter{T}) where {T} = θ.value

function setvalue!(θ::HyperParameter{T}, x::T) where {T}
    checkvalue(θ.interval, x) || error("Value $(x) must be in range " * string(θ.interval))
    θ.value = x
    return θ
end

function setparamoptimizer!(θ::HyperParameter{T},opt::Optimizer) where {T}
    θ.opt = copy(opt)
end

checkvalue(θ::HyperParameter{T}, x::T) where {T} = checkvalue(θ.interval, x)

convert(::Type{HyperParameter{T}}, θ::HyperParameter{T}) where {T<:Real} = θ
function convert(::Type{HyperParameter{T}}, θ::HyperParameter) where {T<:Real}
    HyperParameter{T}(convert(T, getvalue(θ)), convert(Interval{T}, θ.bounds))
end
convert(::Type{T}, θ::HyperParameter{T}) where {T<:Number} = T(θ.value)

function show(io::IO, θ::HyperParameter{T}) where {T}
    print(io, string("HyperParameter(", getvalue(θ), ",", string(θ.interval), ")"))
end

gettheta(θ::HyperParameter) = theta(θ.interval, getvalue(θ))

settheta!(θ::HyperParameter, x::T) where {T}= setvalue!(θ, eta(θ.interval,x))

checktheta(θ::HyperParameter, x::T) where {T} = checktheta(θ.interval, x)

getderiv_eta(θ::HyperParameter) = deriv_eta(θ.interval, getvalue(θ))

for op in (:isless, :(==), :+, :-, :*, :/)
    @eval begin
        $op(θ1::HyperParameter, θ2::HyperParameter) = $op(getvalue(θ1), getvalue(θ2))
        $op(a::Number, θ::HyperParameter) = $op(a, getvalue(θ))
        $op(θ::HyperParameter, a::Number) = $op(getvalue(θ), a)
    end
end
###Old version not using reparametrization
update!(param::HyperParameter{T},grad::T) where {T} = begin
    # println("Correc : $(getderiv_eta(param)), Grad : $(GradDescent.update(param.opt,grad)), theta : $(gettheta(param))")
    isfree(param) ? settheta!(param, gettheta(param) + update(param.opt,getderiv_eta(param)*grad)) : nothing
    # isfree(param) ? setvalue!(param, getvalue(param) + GradDescent.update(param.opt,grad)) : nothing
end

isfree(θ::HyperParameter) = !θ.fixed

setfixed!(θ::HyperParameter) = θ.fixed = true

setfree!(θ::HyperParameter) = θ.fixed = false

mutable struct HyperParameters{T<:Real}
    hyperparameters::Array{HyperParameter{T},1}
    function HyperParameters{T}(θ::Vector{T},intervals::Array{Interval{T,A,B}}) where {A<:Bound{T},B<:Bound{T}} where {T<:Real}
        this = new(Vector{HyperParameter{T}}(undef,length(θ)))
        for i in 1:length(θ)
            this.hyperparameters[i] = HyperParameter{T}(θ[i],intervals[i])
        end
        return this
    end
end

function HyperParameters(θ::Vector{T},intervals::Vector{Interval{T,A,B}}) where {A<:Bound{T},B<:Bound{T}} where {T<:Real}
    HyperParameters{T}(θ,intervals)
end

function setparamoptimizer!(θs::HyperParameters{T},opt::Optimizer) where {T}
    for θ in θs.hyperparameters
        θ.opt = copy(opt)
    end
end

@inline getvalue(θ::HyperParameters{T}) where T = broadcast(getvalue,θ.hyperparameters)

function Base.getindex(p::HyperParameters{T},it::Integer) where T
    return p.hyperparameters[it]
end

function update!(param::HyperParameters{T},grad::Vector{T}) where T
    update!.(param.hyperparameters,grad)
end

setfixed!(θ::HyperParameters{T}) where T = setfixed!.(θ.hyperparameters)

setfree!(θ::HyperParameters{T}) where T = setfree!.(θ.hyperparameters)