text:
using YaoBase, SparseArrays, BitBasis, YaoArrayRegister
export @ket_str, @bra_str

function parse_str(s::String)
    v = 0
    k = 1
    for each in reverse(filter(x -> x != '_', s))
        if each == '1'
            v += 1 << (k - 1)
            k += 1
        elseif each == '0'
            k += 1
        elseif each == '_'
            continue
        else
            error("expect 0 or 1, got $each at $k-th bit")
        end
    end
    return v, k - 1
end

function ket_m(s)
    error(
        "Please install SymEngine and type using SymEngine explicitly to use symbolic functionality of Yao",
    )
end

function bra_m(s)
    adjoint(ket_m(s))
end

"""
    @ket_str

Create a ket register. See also [`@bra_str`](@ref).

# Example

a symbolic quantum state can be created simply by

```jldoctest
julia> ket"110" + 2ket"111"
|110⟩ + 2|111⟩
```

qubits can be partially actived by [`focus!`](@ref)

```jldoctest
julia> ket"100" + ket"111" |> focus!(1:2)
|100⟩ + |111⟩
```

"""
macro ket_str(s)
    ket_m(s)
end

"""
    @bra_str

Create a bra register. See also [`@ket_str`](@ref).

# Example

Similar to `@ket_str` literal, a symbolic quantum state can be created
by

```jldoctest
julia> bra"111" + 2bra"101"
2⟨101| + ⟨111|

julia> bra"111" * (ket"101" + ket"111")
1
```
"""
macro bra_str(s)
    bra_m(s)
end

function print_braket(f, io::IO, r)
    print(io, "|")
    f()
    print(io, "⟩")
end

function print_braket(f, io::IO, r::AdjointArrayReg)
    print(io, "⟨")
    f()
    print(io, "|")
end

function print_basis(io, active::Int, remain::Int, r)
    print_braket(io, r) do
        printstyled(io, string(remain, base = 2, pad = nremain(r)), color = :light_black)
        print(io, string(active, base = 2, pad = nactive(r)))
    end
end

function print_amp(io::IO, x)
    if !isone(x)
        if isempty(free_symbols(x))
            print(io, round(N(x); digits = 2))
        else
            print(io, x)
        end
    end
end

function print_sym_state(io::IO, r::ArrayReg{1})
    st = state(r)
    m, n = size(st)
    isfirst_nonzero = true
    amp = st[1, 1]
    if !iszero(amp)
        print_amp(io, amp)
        print_basis(io, 0, 0, r)
        isfirst_nonzero = false
    end

    for j in 1:n, i in 1:m
        i == 1 && j == 1 && continue
        amp = st[i, j]
        if iszero(amp)
            continue
        end

        isfirst_nonzero || print(io, " + ")
        print_amp(io, amp)
        print_basis(io, i - 1, j - 1, r)
        isfirst_nonzero = false
    end
end

function print_sym_state(io::IO, r::AdjointArrayReg{1})
    st = state(r)
    m, n = size(st)
    isfirst_nonzero = true
    amp = st[1, 1]
    if !iszero(amp)
        print_amp(io, amp)
        print_basis(io, 0, 0, r)
        isfirst_nonzero = false
    end

    for j in 1:n, i in 1:m
        i == 1 && j == 1 && continue
        amp = st[i, j]
        if iszero(amp)
            continue
        end

        isfirst_nonzero || print(io, " + ")
        print_amp(io, amp)
        print_basis(io, j - 1, i - 1, r)
        isfirst_nonzero = false
    end
end

# TODO: need a more extensible implementation
# Base.:(*)(x::SymReg{B, MT}, y::SymReg{B, MT}) where {B, MT} = SymReg{B, MT}(kron(state(x), state(y)))
# Base.:(^)(x::SymReg{B, MT}, n::Int) where {B, MT} = SymReg{B, MT}(kron(state(x) for _ in 1:n))

# Base.:(*)(x::AdjointSymReg{B, MT}, y::AdjointSymReg{B, MT}) where {B, MT} = adjoint(parent(x) * parent(y))
# Base.:(^)(x::AdjointSymReg{B, MT}, n::Int) where {B, MT} = adjoint(parent(x)^n)


text:
using BenchmarkTools, PhyloNetworks, DataFrames, Logging

#suppresses @warn and @info for benchmarks
logger = SimpleLogger(stderr, Logging.Error);
old_logger = global_logger(logger);

# Define a parent BenchmarkGroup to contain our SUITE
const SUITE = BenchmarkGroup()

SUITE["nasm"] = BenchmarkGroup(["JC69", "HKY85"])
SUITE["fitDiscreteFixed"] = BenchmarkGroup(["ERSM", "BTSM", "JC69", "HKY85"])
SUITE["fitdiscrete"] = BenchmarkGroup(["ERSM", "BTSM", "JC69", "HKY85"])

# Add benchmarks to nasm group
SUITE["nasm"]["JC69"] = @benchmarkable JC69([0.5])
m1 = HKY85([.5], [0.25, 0.25, 0.25, 0.25]);
SUITE["nasm"]["HKY85"] = @benchmarkable P!(P(m1, 1.0), m1, 3.0)

# fitDiscreteFixed benchmarks
net_dat = readTopology("(((A:2.0,(B:1.0)#H1:0.1::0.9):1.5,(C:0.6,#H1:1.0::0.1):1.0):0.5,D:2.0);")
species_alone = ["C","A","B","D"]
dat_alone = DataFrame(trait=["hi","lo","lo","hi"])
SUITE["fitDiscreteFixed"]["ERSM"] = @benchmarkable fitdiscrete(net_dat, :ERSM, species_alone, dat_alone; optimizeQ=false, optimizeRVAS=false)
SUITE["fitDiscreteFixed"]["BTSM"] = @benchmarkable fitdiscrete(net_dat, :BTSM, species_alone, dat_alone; optimizeQ=false, optimizeRVAS=false)

fastafile = joinpath(@__DIR__, "..", "examples", "Ae_bicornis_Tr406_Contig10132.aln")
dna_dat, dna_weights = readfastatodna(fastafile, true);
net_dna = readTopology("((((((((((((((Ae_caudata_Tr275,Ae_caudata_Tr276),Ae_caudata_Tr139))#H1,#H2),(((Ae_umbellulata_Tr266,Ae_umbellulata_Tr257),Ae_umbellulata_Tr268),#H1)),((Ae_comosa_Tr271,Ae_comosa_Tr272),(((Ae_uniaristata_Tr403,Ae_uniaristata_Tr357),Ae_uniaristata_Tr402),Ae_uniaristata_Tr404))),(((Ae_tauschii_Tr352,Ae_tauschii_Tr351),(Ae_tauschii_Tr180,Ae_tauschii_Tr125)),(((((((Ae_longissima_Tr241,Ae_longissima_Tr242),Ae_longissima_Tr355),(Ae_sharonensis_Tr265,Ae_sharonensis_Tr264)),((Ae_bicornis_Tr408,Ae_bicornis_Tr407),Ae_bicornis_Tr406)),((Ae_searsii_Tr164,Ae_searsii_Tr165),Ae_searsii_Tr161)))#H2,#H4))),(((T_boeoticum_TS8,(T_boeoticum_TS10,T_boeoticum_TS3)),T_boeoticum_TS4),((T_urartu_Tr315,T_urartu_Tr232),(T_urartu_Tr317,T_urartu_Tr309)))),(((((Ae_speltoides_Tr320,Ae_speltoides_Tr323),Ae_speltoides_Tr223),Ae_speltoides_Tr251))H3,((((Ae_mutica_Tr237,Ae_mutica_Tr329),Ae_mutica_Tr244),Ae_mutica_Tr332))#H4))),Ta_caputMedusae_TB2),S_vavilovii_Tr279),Er_bonaepartis_TB1),H_vulgare_HVens23);");
for edge in net_dna.edge #adds branch lengths
    setLength!(edge,1.0)
    if edge.gamma < 0
        setGamma!(edge, 0.5)
    end
end
SUITE["fitDiscreteFixed"]["JC69"] = @benchmarkable fitdiscrete(net_dna, :JC69, dna_dat, dna_weights; optimizeQ=false, optimizeRVAS=false)
SUITE["fitDiscreteFixed"]["HKY85"] = @benchmarkable fitdiscrete(net_dna, :HKY85, dna_dat, dna_weights; optimizeQ=false, optimizeRVAS=false)

## fitdiscrete benchmarks
SUITE["fitdiscrete"]["ERSM"] = @benchmarkable fitdiscrete(net_dat, :ERSM, species_alone, dat_alone; optimizeQ=true, optimizeRVAS=true)
SUITE["fitdiscrete"]["BTSM"] = @benchmarkable fitdiscrete(net_dat, :BTSM, species_alone, dat_alone; optimizeQ=true, optimizeRVAS=true)
SUITE["fitdiscrete"]["JC69"] = @benchmarkable fitdiscrete(net_dna, :JC69, dna_dat, dna_weights; optimizeQ=true, optimizeRVAS=true)
SUITE["fitdiscrete"]["HKY85"] = @benchmarkable fitdiscrete(net_dna, :HKY85, dna_dat, dna_weights; optimizeQ=true, optimizeRVAS=true)

# If a cache of tuned parameters already exists, use it, otherwise, tune and cache
# the benchmark parameters. Reusing cached parameters is faster and more reliable
# than re-tuning `SUITE` every time the file is included.
paramspath = joinpath(dirname(@__FILE__), "params.json")

if isfile(paramspath)
    loadparams!(SUITE, BenchmarkTools.load(paramspath)[1], :evals);
else
    tune!(SUITE)
    BenchmarkTools.save(paramspath, params(SUITE));
end

global_logger(old_logger) #restores typical logging at end of benchmarks

ext: jl


text:
@require Temporal begin
using TableTraits
using DataValues

immutable TSIterator{T, S}
    source::S
end

TableTraits.isiterable(x::Temporal.TS) = true
TableTraits.isiterabletable(x::Temporal.TS) = true

function TableTraits.getiterator{S<:Temporal.TS}(ta::S)
    col_expressions = Array{Expr,1}()
    df_columns_tuple_type = Expr(:curly, :Tuple)

    # Add column for index
    push!(col_expressions, Expr(:(::), :Index, S.parameters[2]))
    push!(df_columns_tuple_type.args, S.parameters[2])

    etype = eltype(ta.values)
    if ndims(ta.values)==1
        push!(col_expressions, Expr(:(::), ta.fields[1], etype))
        push!(df_columns_tuple_type.args, etype)
    else
        for i in 1:size(ta.values,2)
            push!(col_expressions, Expr(:(::), Symbol(ta.fields[i]), etype))
            push!(df_columns_tuple_type.args, etype)
        end
    end
    t_expr = NamedTuples.make_tuple(col_expressions)

    t2 = :(TSIterator{Float64,Float64})
    t2.args[2] = t_expr
    t2.args[3] = S

    t = eval(t2)

    e_ta = t(ta)

    return e_ta
end

function Base.length{T,TS}(iter::TSIterator{T,TS})
    return length(iter.source)
end

function Base.eltype{T,TS}(iter::TSIterator{T,TS})
    return T
end

Base.eltype(::Type{TSIterator{T,TS}}) where {T,TS} = T

function Base.start{T,TS}(iter::TSIterator{T,TS})
    return 1
end

@generated function Base.next{T,S}(iter::TSIterator{T,S}, state)
    constructor_call = Expr(:call, :($T))

    # Add index column
    push!(constructor_call.args, :(iter.source.index[i]))

    for i in 1:length(T.parameters)-1
        push!(constructor_call.args, :(iter.source.values[i,$i]))
    end

    quote
        i = state
        a = $constructor_call
        return a, state+1
    end
end

function Base.done{T,TS}(iter::TSIterator{T,TS}, state)
    return state>length(iter.source.index)
end

# Sink

# TODO This is a terribly inefficient implementation. Minimally it
# should be changed to be more type stable.
function Temporal.TS(x; index_column::Symbol=:Index)
    isiterabletable(x) || error()

    iter = getiterator(x)

    if TableTraits.column_count(iter)<2
        error("Need at least two columns")
    end

    names = TableTraits.column_names(iter)

    timestep_col_index = findfirst(names, index_column)

    if timestep_col_index==0
        error("No index column found.")
    end

    col_types = TableTraits.column_types(iter)

    data_columns = collect(Iterators.filter(i->i[2][1]!=index_column, enumerate(zip(names, col_types))))

    orig_data_type = data_columns[1][2][2]

    data_type = orig_data_type <: DataValue ? orig_data_type.parameters[1] : orig_data_type

    orig_timestep_type = col_types[timestep_col_index]

    timestep_type = orig_timestep_type <: DataValue ? orig_timestep_type.parameters[1] : orig_timestep_type

    if any(i->i[2][2]!=orig_data_type, data_columns)
        error("All data columns need to be of the same type.")
    end

    t_column = Array{timestep_type,1}()
    d_array = Array{Array{data_type,1},1}()
    for i in data_columns
        push!(d_array, Array{data_type,1}())
    end

    for v in iter
        if orig_timestep_type <: DataValue
            push!(t_column, get(v[timestep_col_index]))
        else
            push!(t_column, v[timestep_col_index])
        end

        if orig_data_type <: DataValue
            for (i,c) in enumerate(data_columns)
                push!(d_array[i],get(v[c[1]]))
            end
        else
            for (i,c) in enumerate(data_columns)
                push!(d_array[i],v[c[1]])
            end
        end
    end

    d_array = hcat(d_array...)

    ta = Temporal.TS(d_array, t_column,[i[2][1] for i in data_columns])
    return ta
end

end


text:
# Autogenerated wrapper script for SARSOP_jll for x86_64-linux-gnu-cxx03
export polgraph, pomdpconvert, pomdpeval, pomdpsim, pomdpsol

JLLWrappers.@generate_wrapper_header("SARSOP")
JLLWrappers.@declare_executable_product(polgraph)
JLLWrappers.@declare_executable_product(pomdpconvert)
JLLWrappers.@declare_executable_product(pomdpeval)
JLLWrappers.@declare_executable_product(pomdpsim)
JLLWrappers.@declare_executable_product(pomdpsol)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        polgraph,
        "bin/polgraph",
    )

    JLLWrappers.@init_executable_product(
        pomdpconvert,
        "bin/pomdpconvert",
    )

    JLLWrappers.@init_executable_product(
        pomdpeval,
        "bin/pomdpeval",
    )

    JLLWrappers.@init_executable_product(
        pomdpsim,
        "bin/pomdpsim",
    )

    JLLWrappers.@init_executable_product(
        pomdpsol,
        "bin/pomdpsol",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# test greeks

let

    # Setup
    d1 = today()
    d2 = d1 + Day(150)
    c1 = EuropeanCall(d2, SingleStock(), 56.5USD)
    core1 = CoreModel(d1, 47.32USD, 0.0, 0.0)
    c2 = EuropeanCall(d2, SingleStock(), 56.5)
    core2 = CoreModel(d1, 47.32, 0.0, 0.0)
    gbm1 = GeomBMModel(core1, 0.1)
    gbm2 = GeomBMModel(core2, 0.1)
    Random.seed!(1)
    mc1 = montecarlo(gbm1, d1:Day(1):d2, 10000)
    Random.seed!(1)
    mc2 = montecarlo(gbm2, d1:Day(1):d2, 10000)

    # Compare with math formulae
    @test delta(gbm1, c1; autodiff = true) ≈ delta(gbm1, c1)
    @test vega(gbm1, c1; autodiff = true) ≈ vega(gbm1, c1)
    @test delta(gbm2, c2; autodiff = true) ≈ delta(gbm2, c2)
    @test vega(gbm2, c2; autodiff = true) ≈ vega(gbm2, c2)
    @test greek(gbm1, c1; metric = :vol) ≈ vega(gbm1, c1)
    @test greek(gbm2, c2; metric = :vol) ≈ vega(gbm2, c2)
    @test greek(gbm1, c1; metric = :strike) ≈ delta(gbm1, c1)
    @test greek(gbm2, c2; metric = :strike) ≈ delta(gbm2, c2)

    # Monte Carlo models greeks - vega
    let
        vol = gbm1.volatility + 1e-8
        gbm3 = GeomBMModel(gbm1.core, vol)
        Random.seed!(1)
        mc3 = montecarlo(gbm3, mc1.dates, size(mc1.paths, 1))
        x = (value(mc3, c1) - value(mc1, c1)) / 1e-8
        Random.seed!(1)
        @test isapprox(x, vega(gbm1, c1, MonteCarloModel, d1:Day(1):d2, 10000), rtol = 1e-4)
        #srand(1)
        #@test isapprox(x, greek(mc1, c1; metric = :vol), rtol = 1e-4)

        vol = gbm2.volatility + 1e-8
        gbm4 = GeomBMModel(gbm2.core, vol)
        Random.seed!(1)
        mc4 = montecarlo(gbm4, mc2.dates, size(mc2.paths, 1))
        x = (value(mc4, c2) - value(mc2, c2)) / 1e-8
        Random.seed!(1)
        @test isapprox(x, vega(gbm2, c2, MonteCarloModel, d1:Day(1):d2, 10000), rtol = 1e-4)
        #srand(1)
        #@test isapprox(x, greek(mc2, c2; metric = :vol), rtol = 1e-4)
    end

    # Monte Carlo models greeks - delta
    let
        s = gbm1.core.startprice + (1e-8 * USD)
        vol = gbm1.volatility
        gbm3 = GeomBMModel(gbm1.core.carrycurve.reference_date, s,
                            gbm1.core.yieldcurve.rate, gbm1.core.carrycurve.rate, vol)
        Random.seed!(1)
        mc3 = montecarlo(gbm3, mc1.dates, size(mc1.paths, 1))
        x = (value(mc3, c1) - value(mc1, c1)) / 1e-8
        Random.seed!(1)
        @test isapprox(x.val, delta(gbm1, c1, MonteCarloModel, d1:Day(1):d2, 10000), rtol = 1e-4)
        #srand(1)
        #@test isapprox(x.val, greek(mc1, c1; metric = :strike), rtol = 1e-4)

        s = gbm2.core.startprice + 1e-8
        vol = gbm2.volatility
        gbm4 = GeomBMModel(gbm2.core.carrycurve.reference_date, s,
                            gbm2.core.yieldcurve.rate, gbm2.core.carrycurve.rate, vol)
        Random.seed!(1)
        mc4 = montecarlo(gbm4, mc2.dates, size(mc2.paths, 1))
        x = (value(mc4, c2) - value(mc2, c2)) / 1e-8
        Random.seed!(1)
        @test isapprox(x, delta(gbm2, c2, MonteCarloModel, d1:Day(1):d2, 10000), rtol = 1e-4)
        #srand(1)
        #@test isapprox(x, greek(mc2, c2; metric = :strike), rtol = 1e-4)
    end

    # Compare rho with forward difference
    let
        r = gbm1.core.yieldcurve.rate + 1e-8
        gbm3 = GeomBMModel(gbm1.core.carrycurve.reference_date, gbm1.core.startprice,
                                 r, gbm1.core.carrycurve.rate, gbm1.volatility)
        x = (value(gbm3, c1) - value(gbm1,c1)) / 1e-8
        @test isapprox(x, rho(gbm1, c1), rtol = 1e-4)
        @test isapprox(x, greek(gbm1, c1; metric = :interest), rtol = 1e-4)

        Random.seed!(1)
        mc3 = montecarlo(gbm3, mc1.dates, size(mc1.paths,1))
        x = (value(mc3, c1) - value(mc1,c1)) / 1e-8
        Random.seed!(1)
        @test isapprox(x, rho(gbm1, c1, MonteCarloModel, d1:Day(1):d2, 10000), rtol = 1e-4)
        #srand(1)
        #@test isapprox(x, greek(mc1, c1; metric = :interest), rtol = 1e-4)

        gbm4 = GeomBMModel(gbm2.core.carrycurve.reference_date, gbm2.core.startprice,
                                 r, gbm2.core.carrycurve.rate, gbm2.volatility)
        x = (value(gbm4, c2) - value(gbm2,c2)) / 1e-8
        @test isapprox(x, rho(gbm2, c2), rtol = 1e-4)
        @test isapprox(x, greek(gbm2, c2; metric = :interest), rtol = 1e-4)

        Random.seed!(1)
        mc4 = montecarlo(gbm4, mc2.dates, size(mc2.paths,1))
        x = (value(mc4, c2) - value(mc2,c2)) / 1e-8
        Random.seed!(1)
        @test isapprox(x, rho(gbm2, c2, MonteCarloModel, d1:Day(1):d2, 10000), rtol = 1e-4)
        #srand(1)
        #@test isapprox(x, greek(mc2, c2; metric = :interest), rtol = 1e-4)
    end

    # Compare gamma with central difference
    let
        s1 = gbm1.core.startprice + (1e-4 * USD)
        s2 = gbm1.core.startprice - (1e-4 * USD)
        gbm31 = GeomBMModel(gbm1.core.carrycurve.reference_date, s1,
                            gbm1.core.yieldcurve.rate, gbm1.core.carrycurve.rate, gbm1.volatility)
        gbm32 = GeomBMModel(gbm1.core.carrycurve.reference_date, s2,
                            gbm1.core.yieldcurve.rate, gbm1.core.carrycurve.rate, gbm1.volatility)
        x = (value(gbm31,c1) + value(gbm32,c1) - (2 * value(gbm1, c1))) / (1e-4)^2
        @test isapprox(x.val, Miletus.gamma(gbm1, c1), rtol = 1e-4)
        @test isapprox(x.val, greek(gbm1, c1; metric = :strike, n = 2), rtol = 1e-4)

        Random.seed!(1)
        mc31 = montecarlo(gbm31, mc1.dates, size(mc1.paths, 1))
        Random.seed!(1)
        mc32 = montecarlo(gbm32, mc1.dates, size(mc1.paths, 1))
        x = (value(mc31,c1) + value(mc32,c1) - (2 * value(mc1, c1))) / (1e-4)^2
        Random.seed!(1)
        @test x.val - Miletus.gamma(gbm1, c1, MonteCarloModel, d1:Day(1):d2, 10000) < 1e-8
        Random.seed!(1)
        @test x.val - greek(gbm1, c1; metric = :strike, n = 2) < 1e-8

        s1 = gbm2.core.startprice + 1e-4
        s2 = gbm2.core.startprice - 1e-4
        gbm41 = GeomBMModel(gbm2.core.carrycurve.reference_date, s1,
                            gbm2.core.yieldcurve.rate, gbm2.core.carrycurve.rate, gbm2.volatility)
        gbm42 = GeomBMModel(gbm2.core.carrycurve.reference_date, s2,
                            gbm2.core.yieldcurve.rate, gbm2.core.carrycurve.rate, gbm2.volatility)
        x = (value(gbm41,c2) + value(gbm42,c2) - (2 * value(gbm2, c2))) / (1e-4)^2
        @test isapprox(x, Miletus.gamma(gbm2, c2), rtol = 1e-4)
        @test isapprox(x, greek(gbm2, c2; metric = :strike, n = 2), rtol = 1e-4)

        Random.seed!(1)
        mc41 = montecarlo(gbm41, mc2.dates, size(mc2.paths, 1))
        Random.seed!(1)
        mc42 = montecarlo(gbm42, mc2.dates, size(mc2.paths, 1))
        x = (value(mc41,c2) + value(mc42,c2) - (2 * value(mc2, c2))) / (1e-4)^2
        Random.seed!(1)
        @test x - Miletus.gamma(gbm2, c2, MonteCarloModel, d1:Day(1):d2, 10000) < 1e-8
        Random.seed!(1)
        @test x - greek(gbm2, c2; metric = :strike, n = 2) < 1e-8
    end

end


text:
module ML

using JuliaDB, Dagger, OnlineStats, PooledArrays, Statistics

import Dagger: ArrayOp, DArray, treereduce
import JuliaDB: Dataset, DDataset, nrows
import Base: merge

# Core schema types

schema(xs::AbstractArray) = nothing # catch-all
schema(xs::AbstractArray, ::Nothing) = nothing # catch-all
merge(::Nothing, ::Nothing) = nothing
width(::Nothing) = 0
featuremat!(A, ::Nothing, xs) = A

schema(xs::ArrayOp) = schema(compute(xs))
schema(xs::ArrayOp, T) = schema(compute(xs), T)

#-----------------------------------------------------------------------# Continuous
struct Continuous
    stat::Variance
end

function schema(xs::AbstractArray, ::Type{Continuous})
    Continuous(fit!(Variance(), xs))
end
function schema(xs::AbstractArray{<:Real})
    schema(xs, Continuous)
end
width(::Continuous) = 1
function merge(c1::Continuous, c2::Continuous)
    Continuous(merge(c1.stat, c2.stat))
end
Statistics.mean(c::Continuous)::Float64 = mean(c.stat)
Statistics.std(c::Continuous)::Float64 = std(c.stat)
function Base.show(io::IO, c::Continuous)
    write(io, "Continous(μ=$(mean(c)), σ=$(std(c)))")
end
Base.@propagate_inbounds function featuremat!(A, c::Continuous, xs, dropmissing=Val(false))
    m = mean(c)
    s = std(c)
    for i in 1:length(xs)
        x =  xs[i]
        if dropmissing isa Val{true}
            if ismissing(x)
                continue
            else
                # A[i, 1] = (get(x) - m) / s
                A[i, 1] = (x - m) / s
            end
        else
            A[i, 1] = (x - m) / s
        end
    end
    A
end

#-----------------------------------------------------------------------# Categorical
struct Categorical
    stat::CountMap
end
function Categorical(xs::AbstractArray{T}) where {T}
    Categorical(CountMap(Dict(zip(xs, 1:length(xs)))))
end

function schema(xs::AbstractArray{T}, ::Type{Categorical}) where {T}
    Categorical(fit!(CountMap(T), xs))
end

function schema(xs::PooledArray)
    schema(xs, Categorical)
end

Base.keys(c::Categorical) = keys(c.stat)
width(c::Categorical) = length(keys(c))
merge(c1::Categorical, c2::Categorical) = Categorical(merge(c1.stat, c2.stat))
function Base.show(io::IO, c::Categorical)
    write(io, "Categorical($(collect(keys(c))))")
end
Base.@propagate_inbounds function featuremat!(A, c::Categorical, xs, dropmissing=Val(false))
    ks = keys(c)
    labeldict = Dict{eltype(ks), Int}(zip(ks, 1:length(ks)))
    for i = 1:length(xs)
        if dropmissing isa Val{true} && ismissing(xs[i])
            continue
        end
        A[i, labeldict[xs[i]]] = one(eltype(A))
    end
    A
end

# distributed schema calculation
function schema(xs::DArray)
    collect(treereduce(delayed(merge), delayedmap(x -> schema(x), xs.chunks)))
end
function schema(xs::DArray, ::Type{Continuous})
    collect(treereduce(delayed(merge), delayedmap(x -> schema(x, Continuous), xs.chunks)))
end
function schema(xs::DArray, ::Type{Categorical})
    collect(treereduce(delayed(merge), delayedmap(x -> schema(x, Categorical), xs.chunks)))
end

struct Maybe{T}
    feature::T
end
function schema(xs::Vector{Union{Missing,S}}, T::Type) where {S}
    Maybe(schema(collect(skipmissing(xs)), T))
end
schema(xs::Vector{Union{T,Missing}}) where {T} = Maybe(schema(collect(skipmissing(xs))))

width(c::Maybe) = width(c.feature) + 1
merge(m1::Maybe, m2::Maybe) = Maybe(merge(m1.feature, m2.feature))
nulls(xs) = map(ismissing, xs)
Base.@propagate_inbounds function featuremat!(A, c::Maybe, xs, dropmissing=Val(true))
    copyto!(A, CartesianIndices((1:length(xs), 1:1)), reshape(nulls(xs), (length(xs), 1)), CartesianIndices((1:length(xs), 1:1)))
    featuremat!(view(A, 1:length(xs), 2:size(A, 2)), c.feature, xs, Val(true))
    A
end

# Schema inference

const Schema = Dict{Symbol,Any}

# vecTs: type of column vectors in each chunk
function schema(cols, names; hints=Dict())
    d = Schema()
    for (col, name) in zip(cols, names)
        if haskey(hints, name)
            d[name] = schema(col, hints[name])
        else
            d[name] = schema(col)
        end
    end
    d
end

function schema(t::Union{Dataset, DDataset}; hints=Dict())
    schema(collect(columns(t)), colnames(t), hints=hints)
end

width(sch::Schema) = sum(width(s) for s in values(sch))
function featuremat!(A, schemas::Schema, t::Dataset)
    j = 0
    for col in keys(schemas)
        schema = schemas[col]
        featuremat!(view(A, 1:length(t), j+1:j+width(schema)),
                   schema, column(t, col))

        j += width(schema)
    end
    A
end

splitschema(xs::Schema, ks...) =
    filter((k,v) -> k ∉ ks, xs),
    filter((k,v) -> k ∈ ks, xs)

function featuremat(sch, xs)
    featuremat!(zeros(Float32, length(xs), width(sch)), sch, xs)'
end
featuremat(t) = featuremat(schema(t), t)

function featuremat(s, t::DDataset)
    t = compute(t)
    w = width(s)
    h = length(t)
    lengths = get.(nrows.(t.domains))
    domains = Dagger.DomainBlocks((1,1), ([w], cumsum(lengths)))

    DArray(
        Float32,
        Dagger.ArrayDomain(w, sum(lengths)),
        domains,
        reshape(delayedmap(x->featuremat(s,x), t.chunks), (1, length(t.chunks)))
    )
end

end


text:
using Test
using Random

function michaelis_menten(u, p, t)
    [0.6 - 1.5u[1]/(0.3+u[1])] # Should be 0.6*0.3+0.6u[1] - 1.5u[1] = u[2]*u[1]-0.3*u[2]
end

u0 = [0.5]
tspan = (0.0, 4.0)

problem_1 = ODEProblem(michaelis_menten, u0, tspan)
solution_1 = solve(problem_1, Tsit5(), saveat = 0.1)
problem_2 = ODEProblem(michaelis_menten, 2*u0, tspan)
solution_2 = solve(problem_2, Tsit5(), saveat = 0.1)
X = [solution_1[:,:] solution_2[:,:]]
ts = [solution_1.t; solution_2.t]

DX = similar(X)
for (i, xi) in enumerate(eachcol(X))
    DX[:, i] = michaelis_menten(xi, [], ts[i])
end

@parameters t
@variables u[1:2]
u = collect(u)
h = [monomial_basis(u[1:1], 4)...]
basis = Basis([h; h .* u[2]], u)

@testset "Ideal data" begin

    prob = ContinuousDataDrivenProblem(X, ts, DX)

    opts = [ImplicitOptimizer(5e-1);ImplicitOptimizer(0.4:0.1:0.7)]
    for opt in opts
        res = solve(prob, basis, opt,u[2:2] ,normalize = false, denoise = false, maxiter = 10000)
        m = metrics(res)
        @test all(m[:L₂] .< 1e-1)
        @test all(m[:AIC] .> 1000.0)
        @test all(m[:R²] .> 0.6)
    end

end

Random.seed!(2345)
X = X .+ 1e-3*randn(size(X))

@testset "Noisy data" begin

    prob = ContinuousDataDrivenProblem(X, ts, GaussianKernel())

    for opt in [ImplicitOptimizer(3e-1);ImplicitOptimizer(3e-1:0.1:7e-1)]
        res = solve(prob, basis, opt, u[2:2], normalize = false, denoise = true)
        m = metrics(res)
        @test all(m[:L₂] .< 1e-1)
        @test all(m[:AIC] .> 1000.0)
        @test all(m[:R²] .> 0.9)
    end

end


text:
export GraphBuilder, nlines, maxorder
export assign_label!, add!, tn_attach!, finish!
export generate_eingraph

mutable struct GraphBuilder{T}
    lines::Vector{Int}
    tensors::Vector{AbstractArray{T,N} where N}
    labels::Vector{Tuple}
    label_counter::Int
end

function swaplines!(gb::GraphBuilder, i::Int, j::Int)
    li = gb.lines[i]
    gb.lines[i] = gb.lines[j]
    gb.lines[j] = li
    return gb
end

function GraphBuilder(::Type{T}, config::Union{Vector{Int}, BitStr}) where T
    n = length(config)
    tensors = (AbstractArray{T,N} where N)[c == 0 ? T[1, 0.0] : T[0.0, 1] for c in config]
    GraphBuilder(collect(1:n), tensors, Vector{Tuple}([(i,) for i=1:n]), n)
end
GraphBuilder(config::Union{Vector{Int}, BitStr}) = GraphBuilder(ComplexF64, config)
GraphBuilder(n::Int) = GraphBuilder(zeros(Int,n))

nlines(gb::GraphBuilder) = length(gb.lines)
ntensors(gb::GraphBuilder) = length(gb.tensors)
maxorder(gb::GraphBuilder) = maximum(ndims, gb.tensors)
function assign_label!(gb::GraphBuilder, lineno::Int)
    olabel = gb.lines[lineno]
    nlabel = gb.label_counter+=1
    gb.lines[lineno] = nlabel
    olabel, nlabel
end

assign_label!(gb::GraphBuilder) = gb.label_counter+=1

function add!(gb::GraphBuilder, tensor::AbstractArray{T,N}, label::NTuple{N,Int}) where {T,N}
    push!(gb.tensors, tensor)
    push!(gb.labels, label)
end

function tn_attach!(gb::GraphBuilder, gate::Pair{Int,<:AbstractArray{T2,2}}) where {T2}
    a, b = assign_label!(gb, gate.first)
    add!(gb, gate.second, (a,b))
    gb
end

function tn_attach!(gb::GraphBuilder, gate::Pair{NTuple{M,Int},<:AbstractArray}; pre::Int=-1, hasnext::Bool=false) where {M}
    @assert M == log2dim1(gate.second)
    as, bs = zip([assign_label!(gb, loc) for loc in gate.first]...)
    prenex = Int[]
    pre > 0 && push!(prenex, pre)
    hasnext && push!(prenex, assign_label!(gb))
    tensor = reshape(gate.second, fill(2,2*M+length(prenex))...)
    add!(gb, tensor, (as..., bs..., prenex...))
    gb
end

function tn_attach!(gb::GraphBuilder, gates::Pair...)
    ng = length(gates)
    ng == 0 && return gb

    local b
    for i=1:ng
        gate = gates[i]
        tn_attach!(gb, gate; pre=i==1 ? -1 : gb.label_counter, hasnext=i!=ng)
    end
    gb
end

function finish!(gb::GraphBuilder, config::Union{Vector{Int}, BitStr})
    @assert length(config) == nlines(gb)
    for i=1:nlines(gb)
        add!(gb, config[i]==1 ? [0.0im, 1] : [1, 0.0im], (gb.lines[i],))
    end
    gb.lines .= -1
    gb
end

function tn_pop!(gb::GraphBuilder)
    t = pop!(gb.tensors)
    l = pop!(gb.labels)
    # tensors labels must be arrange as `left` first, then `right` and finally `up` and `down`.
    nconnect = findfirst(x->x in gb.lines, l) - 1
    if !(nconnect isa Nothing)
        for i=1:nconnect
            # find and restore the line
            for j=1:nlines(gb)
                if gb.lines[j] == l[nconnect+i]
                    gb.lines[j] = l[i]
                end
            end
        end
    end
    return (t, l)
end

function reset_final!(gb::GraphBuilder, config::Union{Vector{Int},BitStr})
    n = nlines(gb)
    nt = ntensors(gb)
    gb.lines = [(@assert length(gb.labels[nt-n+i]) == 1; gb.labels[nt-n+i][1]) for i=1:n]
    gb.tensors = gb.tensors[1:end-n]
    gb.labels = gb.labels[1:end-n]
    finish!(gb, config)
end

function generate_eingraph(gb::GraphBuilder)
    EinGraph([gb.tensors...], gb.labels)
end
maxorder(gb::EinGraph) = maximum(ndims, gb.tensors)


text:
"""
RhoPorosity(args...)

Calculate the porosity from density log

`` { \\phi  }_{ \\rho  }=\\frac { { \\rho  }_{ ma }-{ \\rho  }_{ log } }{ { \\rho  }_{ ma }-{ \\rho  }_{ Fluid } } ``

## Input Example

```julia
RhoPorosity(RhoLog,RhoMa,RhoF)
```
"""
function RhoPorosity(RhoLog,RhoMa,RhoF)
    phi=(RhoMa.-RhoLog)./(RhoMa.-RhoF)
end


text:
module Kabsch
export calc_centroid, kabsch, rotate, rmsd, translate_points, kabsch_rmsd
	# Calculate root mean square deviation of two matrices A, B
	# http://en.wikipedia.org/wiki/Root-mean-square_deviation_of_atomic_positions
	function rmsd(A::Array{Float64,2}, B::Array{Float64,2})

		RMSD::Float64 = 0.0

		# D pairs of equivalent atoms
		D::Int = size(A)[1]::Int # <- oddly _only_ changing this to Int makes it work on 32-bit systems.
		# N coordinates
		N::Int = length(A)::Int

		for i::Int64 = 1:N
			RMSD += (A[i]::Float64 - B[i]::Float64)^2
		end
		return sqrt(RMSD / D)
	end

	# calculate a centroid of a matrix
	function calc_centroid(m::Array{Float64,2})

		sum_m::Array{Float64,2} = sum(m,1)
		size_m::Int64 = size(m)[1]

		return map(x -> x/size_m, sum_m)
	end

	# Translate P, Q so centroids are equal to the origin of the coordinate system
	# Translation der Massenzentren, so dass beide Zentren im Ursprung des Koordinatensystems liegen
	function translate_points(P::Array{Float64,2}, Q::Array{Float64,2})
		# Calculate centroids P, Q
		# Die Massenzentren der Proteine
		centroid_p::Array{Float64,2} = calc_centroid(P)
		centroid_q::Array{Float64,2} = calc_centroid(Q)

		P = broadcast(-,P, centroid_p)

		Q = broadcast(-,Q, centroid_q)

		return P, Q, centroid_p, centroid_q
	end

	# Input: Two sets of points: reference, coords as Nx3 Matrices (so)
	# returns optimally rotated matrix
	function kabsch(reference::Array{Float64,2}, coords::Array{Float64,2})

		centered_reference::Array{Float64,2}, centered_coords::Array{Float64,2}, centroid_p::Array{Float64,2}, centroid_q::Array{Float64,2}  = translate_points(reference, coords)
		# Compute covariance matrix A
		A::Array{Float64,2} = *(centered_coords', centered_reference)

		# Calculate Singular Value Decomposition (SVD) of A
		u::Array{Float64,2}, d::Array{Float64,1}, vt::Array{Float64,2} = svd(A)

		# check for reflection
		f::Int64 = sign(det(vt) * det(u))
		m::Array{Int64,2} = [1 0 0; 0 1 0; 0 0 f]

		# Calculate the optimal rotation matrix _and_ superimpose it
		return broadcast(+, *(centered_coords, u, m, vt'), centroid_p)

	end

	# directly return RMSD for matrices P, Q for convenience
	function kabsch_rmsd(P::Array{Float64,2}, Q::Array{Float64,2})
		return rmsd(P,kabsch(P,Q))
	end

end


text:
"""
    pass1(project)

 - Remove `pass1/`
 - Copy files from `original/` into `pass1`
 - Convert file extensions (.py, .pyx, .pxd) to .jl, move into `src/`
 - Merge .pxd and .pyx to single file
 - Create `src/project.jl` and include all other files
"""
function pass1(project)
    println("Performing pass1...")
    println("   Copying project...")
    # Clean
    rm(joinpath(project.paths[:example_root], "pass1", "src");recursive=true,force=true)
    # Convert file extension
    # If .pxd and .pyx have the same name,
    # merge the files into one .jl
    files = [joinpath(root, f) for
            (root, dirs, files) in
            Base.Filesystem.walkdir(
            project.paths[:local_original])
            for f in files]

    mkpath(joinpath(project.paths[:example_root], "pass1", "src"))

    println("   Converting file extensions...")
    cd(joinpath(project.paths[:example_root], "pass1", "src")) do
        for file in files
            file_new = basename(file)
            for ext in extensions_from(project.language_from)
                if endswith(file_new, ext)
                    file_new = replace(file_new, ext=>extension_to(project.language_to))
                    cp(file, file_new;force=true)
                end
            end
        end
    end
    println("   Merging files...")
    cd(joinpath(project.paths[:example_root], "pass1", "src")) do
        for file in files
            file_new = basename(file)
            ext = extension_from_append(project.language_from)
            if endswith(file_new, ext)
                file_new = replace(file_new, ext=>extension_to(project.language_to))
                contents_new = ""*open(f->read(f, String), file)*"\n #pyx files\n"
                if isfile(file_new)
                    open(file_new, "a") do io
                        print(io, contents_new)
                    end
                end
            end
        end
    end
    files = [joinpath(root, f) for
            (root, dirs, files) in
            Base.Filesystem.walkdir(joinpath(project.paths[:example_root], "pass1", "src")) for f in files]
    project_src =
        joinpath(project.paths[:example_root],
            "pass1", "src",project.name)

    # TODO: include order requires topological sort of all files
    if !isfile(project_src)
        println("   Creating new project files...")
        open(project_src, "w") do io
            println(io, "module $(project.base_name)\n\n")
            for f in files
                println(io, "include(\"$(basename(f))\")")
            end
            println(io, "\n\nend")
        end
    end
end


text:
# Author: Mathias Louboutin, mlouboutin3@gatech.edu
# Date: December 2020

mutable struct BregmanParams
    verbose
    progTol
    maxIter
    store_trace
    antichatter
    quantile
    alpha
    spg
end

"""
    bregman_options(;verbose=1, optTol=1e-6, progTol=1e-8, maxIter=20
                    store_trace=false, linesearch=false, alpha=.25, spg=false)

Options structure for the bregman iteration algorithm

# Arguments

- `verbose`: level of verbosity (0: no output, 1: final, 2: iter (default), 3: debug)
- `progTol`: tolerance used to check for lack of progress (default: 1e-9)
- `maxIter`: maximum number of iterations (default: 20)
- `store_trace`: Whether to store the trace/history of x (default: false)
- `antichatter`: Whether to use anti-chatter step length correction
- `quantile`: Thresholding level as quantile value, (default=.95 i.e thresholds 95% of the vector)
- `alpha`: Strong convexity modulus. (step length is ``α \\frac{||r||_2^2}{||g||_2^2}``)

"""
bregman_options(;verbose=1, progTol=1e-8, maxIter=20, store_trace=false, antichatter=true, quantile=.95, alpha=.5, spg=false) =
                BregmanParams(verbose, progTol, maxIter, store_trace, antichatter, quantile, alpha, spg)

"""
    bregman(A, TD, x, b, options)

Linearized bregman iteration for the system

``\\frac{1}{2} ||TD \\ x||_2^2 + λ ||TD \\ x||_1  \\ \\ \\ s.t Ax = b``

For example, for sparsity promoting denoising (i.e LSRTM)

# Arguments

- `TD`: curvelet transform
- `A`: Forward operator (J or preconditioned J for LSRTM)
- `b`: observed data
- `x`: Initial guess
"""
function bregman(A, TD, x::Array{T}, b, options) where {T}
    # residual function wrapper
    function obj(x)
        d = A*x
        fun = .5*norm(d - b)^2
        grad = A'*(d - b)
        return fun, grad
    end

    return bregman(obj, x, options, TD)
end

"""
    bregman(fun, TD, x, b, options)

Linearized bregman iteration for the system

``\\frac{1}{2} ||TD \\ x||_2^2 + λ ||TD \\ x||_1  \\ \\ \\ s.t Ax = b``

For example, for sparsity promoting denoising (i.e LSRTM)

# Arguments

- `TD`: curvelet transform
- `fun`: residual function, return the tuple (``f = \\frac{1}{2}||Ax - b||_2``, ``g = A^T(Ax - b)``)
- `b`: observed data
- `x`: Initial guess

"""
function bregman(funobj::Function, x::AbstractArray{T}, options::BregmanParams, TD=nothing) where {T}
    # Output Parameter Settings
    if options.verbose > 0
        @printf("Running linearized bregman...\n");
        @printf("Progress tolerance: %.2e\n",options.progTol)
        @printf("Maximum number of iterations: %d\n",options.maxIter)
        @printf("Anti-chatter correction: %d\n",options.antichatter)
    end
    isnothing(TD) && (TD = Matrix{T}(I, length(x), length(x)))
    # Intitalize variables
    z = TD*x
    d = similar(z)
    options.spg && (gold = similar(x); xold=similar(x))
    if options.antichatter
        tk = 0 * z
    end

    # Result structure
    sol = breglog(x, z)
    # Initialize λ
    λ = T(0)

    # Output Log
    if options.verbose > λ
        @printf("%10s %15s %15s %15s %5s\n","Iteration","Step Length", "L1-2    ", "||A*x - b||_2^2", "λ")
    end

    # Iterations
    for i=1:options.maxIter
        f, g = funobj(x)
        # Preconditionned ipdate direction
        d .= -TD*g
        # Step length
        t = (options.spg && i> 1) ? T(dot(x-xold, x-xold)/dot(x-xold, g-gold)) : T(options.alpha*f/norm(d)^2)

        # Anti-chatter
        if options.antichatter
            @. tk = tk - sign.(d)
            # Chatter correction
            inds_z = findall(abs.(z) .> λ)
            mul!(d, d, t)
            d[:] .*= abs.(tk[:])/i
        end
        # Update z variable
        @. z = z + d
        # Get λ at first iteration
        i == 1 && (λ = T(quantile(abs.(z), options.quantile)))
        # Save curent state
        options.spg && (gold .= g; xold .= x)
        # Update x
        x = TD'*soft_thresholding(z, λ)

        obj_fun = λ * norm(z, 1) + .5 * norm(z, 2)^2
        if options.verbose > 0
            @printf("%10d %15.5e %15.5e %15.5e %15.5e \n",i, t, obj_fun, f, λ)
        end
        norm(x - sol.x) < options.progTol && (@printf("Step size below progTol\n"); break;)
        update!(sol; iter=i, ϕ=obj_fun, residual=f, x=x, z=z, g=g, store_trace=options.store_trace)
    end
    return sol
end

# Utility functions
"""
Simplified Quantile from Statistics.jl since we only need simplified version of it.
"""
function quantile(u::AbstractVector, p::Real)
    0 <= p <= 1 || throw(ArgumentError("input probability out of [0,1] range"))
    n = length(u)
    v = sort(u; alg=Base.QuickSort)

    m = 1 - p
    aleph = n*p + oftype(p, m)
    j = clamp(trunc(Int, aleph), 1, n-1)
    γ = clamp(aleph - j, 0, 1)

    n == 1 ? a = v[1] : a = v[j]
    n == 1 ? b = v[1] : b = v[j+1]

    (isfinite(a) && isfinite(b)) ? q = a + γ*(b-a) : q = (1-γ)*a + γ*b
    return q
end

"""
Bregman result structure
"""
mutable struct BregmanIterations
    x
    z
    g
    ϕ
    residual
    ϕ_trace
    r_trace
    x_trace
    z_trace
end

function update!(r::BregmanIterations; x=nothing, z=nothing, ϕ=nothing, residual=nothing, g=nothing, iter=1, store_trace=false)
    ~isnothing(x) && copyto!(r.x, x)
    ~isnothing(z) && copyto!(r.z, z)
    ~isnothing(ϕ) && (r.ϕ = ϕ)
    ~isnothing(residual) && (r.residual = residual)
    ~isnothing(g) && copyto!(r.g, g)
    (~isnothing(x) && length(r.x_trace) == iter-1 && store_trace) && (push!(r.x_trace, x))
    (~isnothing(z) && length(r.z_trace) == iter-1 && store_trace) && (push!(r.z_trace, z))
    (~isnothing(ϕ) && length(r.ϕ_trace) == iter-1) && (push!(r.ϕ_trace, ϕ))
    (~isnothing(residual) && length(r.r_trace) == iter-1) && (push!(r.r_trace, residual))
end

function breglog(init_x, init_z; f0=0, obj0=0)
    return BregmanIterations(1*init_x, 1*init_z, 0*init_z, f0, obj0, Vector{}(), Vector{}(), Vector{}(), Vector{}())
end


text:
#-----------------------------------#
# DIFFERENTIAL CORRECTION FUNCTIONS #
#-----------------------------------#

function corrector_solve end
function corrector_callback end

function DiffEqBase.__solve(state::State, corrector::Abstract_DifferentialCorrector;
                            dc_maxiters=10, dc_tolerance=1e-6, dc_xtol=dc_tolerance^2,
                            dc_method=:newton, dc_linesearch=LineSearches.BackTracking(),
                            verbose=false, kwargs...)

    # Run once with strict=true to retrieve the initial guess.
    correction_state = convert_to_frame(state, corrector.frame)
    guess_sol, guess_tspan = corrector_solve(corrector, correction_state; strict=true, verbose, kwargs...)

    # The state we will correct begins with a tspan up to the x-crossing
    correction_state = remake(correction_state; tspan=guess_tspan)
    u0 = deepcopy(correction_state.u0)

    # Solve iteratively to find a corrected.
    sol = nlsolve(
        # Function for computing residuals (F) and Jacobian (J)
        only_fj!((F, J, x) -> corrector(correction_state, F, J, x; verbose, kwargs...)),
        # Initial guess (x) -- includes the initial guess tspan's final time
        [u0[corrector.u0_free]..., 0.];
        # Options
        method=dc_method, linesearch=dc_linesearch, ftol=dc_tolerance, xtol=dc_xtol,
        show_trace=verbose, extended_trace=verbose, iterations=dc_maxiters
    )
    verbose && @info sol

    # Propagate out the last solution.
    u0[corrector.u0_free] .= sol.zero[1:length(corrector.u0_free)]
    tspan = (guess_tspan[begin], guess_tspan[end] + sol.zero[end])
    final_state = remake(correction_state; u0, tspan)
    traj = solve(final_state, corrector.alg; callback=corrector_callback(corrector, final_state.model), kwargs...)

    return traj
end

#----------#
# INCLUDES #
#----------#

include("correctors/axisymmetric.jl")
include("correctors/continuation.jl")


text:
"""
    vocab_reduction(wv::WordVectors, seed, nn)

Produces a reduced vocabulary version of `wv` by removing all but
the `nn` nearest neighbors of each word present in the
vocabulary `seed`.
"""
function vocab_reduction(wv::WordVectors, seed, nn)
    idxs = Int[]
    for word in unique(seed)
        if word in wv.vocab
            n = _get_neighbors_number(word, nn)
            pos, _ = cosine_vec(wv, Word2Vec.get_vector(wv, word), n)
            push!(idxs, pos...)
        end
    end
    unique!(idxs)
    return WordVectors(wv.vocab[idxs], wv.vectors[:, idxs])
end

_get_neighbors_number(word, nn::Int) = nn  # fixed number of neighbors
_get_neighbors_number(word::T, nn::Dict{T, Int}) where {T} = nn[word]  # different neighbor number for each vocab entry

max_stars_repo_name: zgornel/EmbeddingsAnalysis

max_issues_repo_name: zgornel/EmbeddingsAnalysis

max_forks_repo_name: zgornel/EmbeddingsAnalysis


text:
"""
A "simple" model that merely wraps around an underlying `MathOptInterface`
model.

    struct MOIModel{T <: MOI.AbstractOptimizer} <: AbstractModel
        opt::T
    end

Some pieces of the [functional interface](@ref functional_interface) are
implemented:

    optimize!(model::MOIModel)
    tear_down(model::MOIModel, instance::AbstractInstance, ::Type{MILPResult})
"""
struct MOIModel{T<:MOI.AbstractOptimizer} <: AbstractModel
    opt::T
end

function optimize!(model::MOIModel)
    MOI.optimize!(model.opt)
    return nothing
end

function tear_down(model::MOIModel, instance::AbstractInstance, ::Type{MILPResult})
    opt = model.opt
    result = MILPResult(
        termination_status = MOI.get(opt, MOI.TerminationStatus()),
        primal_status = MOI.get(opt, MOI.PrimalStatus()),
        solve_time_sec = MOI.get(opt, MOI.SolveTime()),
    )
    if result.primal_status == MOI.FEASIBLE_POINT
        num_vars = MOI.get(opt, MOI.NumberOfVariables())
        result.x = [
            MOI.get(opt, MOI.VariablePrimal(), MOI.VariableIndex(i)) for
            i in 1:num_vars
        ]
    end
    try
        primal_bound = MOI.get(opt, MOI.ObjectiveValue())
        result.primal_bound = primal_bound
    catch ArgumentError
        @warn "Primal bound not available; proceeding anyway."
    end
    try
        dual_bound = MOI.get(opt, MOI.ObjectiveBound())
        result.dual_bound = dual_bound
    catch ArgumentError
        @warn "Dual bound not available; proceeding anyway."
    end
    try
        node_count = MOI.get(opt, MOI.NodeCount())
        result.node_count = node_count
    catch ArgumentError
        @warn "Node count not available; proceeding anyway."
    end
    try
        simplex_iters = MOI.get(opt, MOI.SimplexIterations())
        result.simplex_iters = simplex_iters
    catch ArgumentError
        @warn "Simplex iteration count not available; proceeding anyway."
    end
    return result
end


text:
### estimators for π0 (pi0) ###

## Storey estimator ##

"""
Storey π0 estimator

**Parameters**

- λ : tuning parameter, FloatingPoint, default: 0.1

**Examples**

```julia
Storey()
Storey(0.1)
```

**References**

Storey, JD (2002). "A Direct Approach to False Discovery Rates." Journal of the
Royal Statistical Society, doi:10.1111/1467-9868.00346

"""
type Storey <: Pi0Estimator
    λ::AbstractFloat

    Storey(λ) = isin(λ, 0., 1.) ? new(λ) : throw(DomainError())
end

Storey() = Storey(0.1)

function estimate_pi0{T<:AbstractFloat}(pValues::Vector{T}, pi0estimator::Storey)
    storey_pi0(pValues, pi0estimator.λ)
end

function storey_pi0{T<:AbstractFloat}(pValues::Vector{T}, lambda::T)
    validPValues(pValues)
    pi0 = (sum(pValues .>= lambda) / length(pValues)) / (1.-lambda)
    pi0 = min(pi0, 1.)
    return pi0
end

## Storey bootstrap estimator ##
"""
Storey closed-form bootstrap π0 estimator

StoreyBootstrap(λseq, q)

Reference: David Robinson, 2012
"""
type StoreyBootstrap <: Pi0Estimator
    ## check range of arguments
    λseq::Vector{AbstractFloat}
    q   ::AbstractFloat

    StoreyBootstrap(λseq, q) =
        isin(λseq, 0., 1.) && isin(q, 0., 1.) ? new(λseq, q) : throw(DomainError())
end

StoreyBootstrap() = StoreyBootstrap(collect(0.05:0.05:0.95), 0.1)

function estimate_pi0{T<:AbstractFloat}(pValues::Vector{T}, pi0estimator::StoreyBootstrap)
    bootstrap_pi0(pValues, pi0estimator.λseq, pi0estimator.q)
end

function bootstrap_pi0{T<:AbstractFloat,S<:AbstractFloat}(pValues::Vector{T}, lambda::Vector{S} = [0.05:0.05:0.95;], q::S = 0.1)
    #validPValues(pValues)
    n = length(pValues)
    w = Int[sum(pValues .>= l) for l in lambda] ## TODO: check >= or >
    pi0 = w ./ n ./ (1. - lambda)
    min_pi0 = quantile(pi0, q)
    mse = (w ./ (n.^2 .* (1. - lambda).^2 )) .* (1. - w/n) + (pi0 - min_pi0).^2
    pi0 = min(pi0[indmin(mse)], 1.)
    pi0
end

## Least SLope (LSL) estimator ##

"""
Least SLope (LSL) π0 estimator

LeastSlope()
"""
type LeastSlope <: Pi0Estimator
end

function estimate_pi0{T<:AbstractFloat}(pValues::Vector{T}, pi0estimator::LeastSlope)
    lsl_pi0(pValues)
end

function lsl_pi0{T<:AbstractFloat}(pValues::Vector{T})
    n = length(pValues)
    if !issorted(pValues)
        pValues = sort(pValues)
    end
    s0 = lsl_slope(1, n, pValues)
    sx = 0.
    for i in 2:n
        s1 = lsl_slope(i, n, pValues)
        if (s1 - s0) < 0.
            sx = s1
            break
        end
        s0 = s1
    end
    pi0 = min( 1/sx + 1, n ) / n
    return(pi0)
end

function lsl_slope{T<:AbstractFloat}(i::Int, n::Int, pval::Vector{T})
    s = (1 - pval[i]) / (n - i + 1)
    return s
end

## alternative, vectorized version
## used for comparison and compactness
function lsl_pi0_vec{T<:AbstractFloat}(pValues::Vector{T})
    n = length(pValues)
    if !issorted(pValues)
        pValues = sort(pValues)
    end
    s = (1 - pValues) ./ (n - collect(1:n) + 1)
    d = diff(s) .< 0
    idx = findfirst(d) + 1
    pi0 = min( 1/s[idx] + 1, n ) / n
    return(pi0)
end

## Oracle
"""
Oracle π0

Oracle(π0)
"""
type Oracle <: Pi0Estimator
    π0::AbstractFloat

    Oracle(π0) = isin(π0, 0., 1.) ? new(π0) : throw(DomainError())
end

Oracle() = Oracle(1.0)

function estimate_pi0{T<:AbstractFloat}(pValues::Vector{T}, pi0estimator::Oracle)
    pi0estimator.π0
end

## Two-Step estimator: Benjamini, Krieger and Yekutieli (2006) ##

"""
Two step π0 estimator

TwoStep(α)

Reference: Benjamini, Krieger and Yekutieli, 2006
"""
type TwoStep <: Pi0Estimator
    α::AbstractFloat
    ## method::PValueAdjustmentMethod

    TwoStep(α) = isin(α, 0., 1.) ? new(α) : throw(DomainError())
end

TwoStep() = TwoStep(0.05)

function estimate_pi0{T<:AbstractFloat}(pValues::Vector{T}, pi0estimator::TwoStep)
    twostep_pi0(pValues, pi0estimator.α)
end

function twostep_pi0{T<:AbstractFloat}(pValues::Vector{T}, alpha::T, method::PValueAdjustmentMethod = BenjaminiHochberg())
    padj = adjust(pValues, method)
    pi0 = sum(padj .>= (alpha/(1+alpha))) / length(padj)
    return(pi0)
end

# RightBoundary procedure as defined in Definition 2 of Liang and Nettleton 2012
# "Adaptive and dynamic adaptive procedures for false discovery rate control and estimation"
"""
Right boundary π0 estimator

RightBoundary(λseq)
"""
type RightBoundary <: Pi0Estimator
    ## check range of arguments
    λseq::Vector{Float64}

    RightBoundary(λseq) =
        isin(λseq, 0., 1.) ? new(λseq) : throw(DomainError())
end

# λseq used in Liang, Nettleton 2012
RightBoundary() = RightBoundary([collect(0.02:0.02:0.1); collect(0.15:0.05:0.95)])

function estimate_pi0{T<:AbstractFloat}(pValues::Vector{T}, pi0estimator::RightBoundary)
    rightboundary_pi0(pValues, pi0estimator.λseq)
end

function rightboundary_pi0{T<:AbstractFloat}(pValues::Vector{T}, λseq)
    validPValues(pValues)
    n = length(pValues)
    if !issorted(λseq)
        λseq = sort(λseq)
    end
    # make sure we catch p-values equal to 1 despite left closure
    # use closed=:left because we have been using >= convention in this package
    # note that original paper uses > convention.
    h = fit(Histogram, pValues, [λseq; Inf], closed=:left)
    pi0_estimates = reverse(cumsum(reverse(h.weights)))./(1.-λseq)./n
    pi0_decrease = diff(pi0_estimates) .>= 0
    pi0_decrease[end] = true
    pi0 = pi0_estimates[findfirst(pi0_decrease, true) + 1]
    return(min(pi0,1))
end

## Censored BUM
"""
Censored BUM π0 estimator

CensoredBUM(γ0, λ, xtol, maxiter)
"""
type CensoredBUM <: Pi0Estimator
    γ0::Float64
    λ::Float64
    xtol::Float64
    maxiter::Int64

    function CensoredBUM(γ0, λ, xtol, maxiter)
        if isin(γ0, 0., 1.) && isin(λ, 0., 1.) && isin(xtol, 0., 1.) && maxiter > 0
            new(γ0, λ, xtol, maxiter)
        else
            throw(DomainError())
        end
    end
end

CensoredBUM() = CensoredBUM(0.5, 0.05, 1e-6, 10000)
CensoredBUM(γ0, λ) = CensoredBUM(γ0, λ, 1e-6, 10000)

type CensoredBUMFit <: Pi0Fit
    π0::Float64
    param::Vector{Float64}
    is_converged::Bool
end

function fit(pi0estimator::CensoredBUM, pValues;
             kw...)
    π0, param, is_converged = cbum_pi0(pValues, pi0estimator.γ0, pi0estimator.λ,
                                       pi0estimator.xtol, pi0estimator.maxiter;
                                       kw...)
    return CensoredBUMFit(π0, param, is_converged)
end

function estimate_pi0{T<:AbstractFloat}(pValues::Vector{T}, pi0estimator::CensoredBUM)
    estimate_pi0(fit(pi0estimator, pValues))
end

function estimate_pi0(pi0fit::CensoredBUMFit)
    π0 = pi0fit.is_converged ? π0 = pi0fit.π0 : NaN
    return π0
end

function cbum_pi0{T<:AbstractFloat}(pValues::Vector{T}, γ0::T = 0.5, λ::T = 0.05,
                                    xtol::Float64 = 1e-6, maxiter::Int = 10000)
    n = length(pValues)
    idx_right = pValues .>= λ
    n2 = sum(idx_right)
    n1 = n - n2
    sz = (1-γ0)*n
    szr = (1-γ0)*n2
    szl = sz - szr
    ## compute constant values only once
    pr = pValues[idx_right]
    lpr = log(pr)
    ll = log(λ)
    zr = fill(1-γ0, n2)
    pi0_old = γ0 = α = γ = Inf
    for i in 1:maxiter
        γ = 1 - sz/n
        α = -szr / ( ll * szl + sum(zr .* lpr) )
        γ = max(min(γ, 1.), 0.)
        α = max(min(α, 1.), 0.)
        xl = (1-γ) * (λ^α)
        szl = (xl ./ (γ*λ + xl)) * n1
        xr = (1-γ) * α * pr.^(α-1)
        zr = xr ./ (γ + xr)
        szr = sum(zr)
        sz = szl + szr
        pi0_new = γ + (1-γ)*α
        if abs(pi0_new - pi0_old) <= xtol
            return pi0_new, [γ, α], true
        end
        γ0 = γ
        pi0_old = pi0_new
    end
    return NaN, [γ, α], false
end

function cbum_pi0_naive(pValues, γ0 = 0.5, λ = 0.05,
                        xtol = 1e-6, maxiter = 10000)
    n = length(pValues)
    z = fill(1-γ0, n)
    idx_left = pValues .< λ
    idx_right = !idx_left
    pi0_old = γ0 = α = γ = Inf
    ## compute constant values only once
    lpr = log(pValues[idx_right])
    ll = log(λ)
    for i in 1:maxiter
        γ = sum(1-z) / n
        α = -sum(z[idx_right])
        α = α / ( ll * sum(z[idx_left]) + sum(z[idx_right] .* lpr) )
        xl = (1-γ) * (λ^α)
        z[idx_left] = xl ./ (γ*λ + xl)
        xr = (1-γ) * α * pValues[idx_right].^(α-1)
        z[idx_right] = xr ./ (γ + xr)
        pi0_new = γ + (1-γ)*α
        if abs(pi0_new - pi0_old) <= xtol
            return pi0_new, [γ, α], true
        end
        γ0 = γ
        pi0_old = pi0_new
    end
    return NaN, [γ, α], false
end

## BUM
"""
BUM π0 estimator

BUM(γ0, xtol, maxiter)
"""
type BUM <: Pi0Estimator
    γ0::Float64
    xtol::Float64
    maxiter::Int64

    function BUM(γ0, xtol, maxiter)
        if isin(γ0, 0., 1.) && isin(xtol, 0., 1.)
            new(γ0, xtol, maxiter)
        else
            throw(DomainError())
        end
    end
end

BUM() = BUM(0.5, 1e-6, 10000)

BUM(y0::Float64) = BUM(y0, 1e-6, 10000)

type BUMFit <: Pi0Fit
    π0::Float64
    param::Vector{Float64}
    is_converged::Bool
end

function fit(pi0estimator::BUM, pValues;
             kw...)
    π0, param, is_converged = cbum_pi0(pValues, pi0estimator.γ0, eps(),
                                       pi0estimator.xtol, pi0estimator.maxiter;
                                       kw...)
    return BUMFit(π0, param, is_converged)
end

function estimate_pi0{T<:AbstractFloat}(pValues::Vector{T}, pi0estimator::BUM)
    estimate_pi0(fit(pi0estimator, pValues))
end

function estimate_pi0(pi0fit::BUMFit)
    π0 = pi0fit.is_converged ? π0 = pi0fit.π0 : NaN
    return π0
end

## Longest constant interval in the Grenander estimator: Langaas et al., 2005
# Leave unexported until thoroughly tested

"""
Flat Grenander π0 estimator

FlatGrenander()

Estimates π0 by the longest constant interval in the Grenander estimator

Reference: Langaas et al., 2005: section 4.3
"""
type FlatGrenander <: Pi0Estimator
end

function estimate_pi0{T<:AbstractFloat}(pValues::Vector{T}, pi0estimator::FlatGrenander)
    flat_grenander_pi0(pValues)
end

function flat_grenander_pi0(pv::Vector{Float64})
    p, f, F = grenander(pv)
    pi0 = longest_constant_interval(p, f)
    return pi0
end

function longest_constant_interval(p::Vector{Float64}, f::Vector{Float64})
    p = [0.0; p]
    f = [Inf; f]
    i2 = length(f)
    Δp_max = Δp = 0.0
    pi0 = 1.0
    for i1 in (length(f)-1):-1:1
        if f[i2] ≈ f[i1] ## within constant interval
            Δp = p[i2] - p[i1]
        else
            if Δp >= Δp_max
                Δp_max = Δp
                pi0 = f[i2]
            end
            i2 = i1
            Δp = 0.0
        end
        if f[i1] > 1.0
            break
        end
    end
    return pi0
end


text:
module Rebugger

using UUIDs, InteractiveUtils
using REPL
import REPL.LineEdit, REPL.Terminals
using REPL.LineEdit: buffer, bufend, content, edit_splice!
using REPL.LineEdit: transition, terminal, mode, state

using CodeTracking, Revise, JuliaInterpreter, HeaderREPLs
using Revise: RelocatableExpr, striplines!, printf_maxsize, whichtt, hasfile, unwrap
using JuliaInterpreter: FrameCode, scopeof
using Base.Meta: isexpr
using Core: CodeInfo

# Reexports
export @breakpoint, breakpoint, enable, disable, remove, break_on, break_off

const msgs = []  # for debugging. The REPL-magic can sometimes overprint error messages

include("debug.jl")
include("ui.jl")
include("printing.jl")
include("deepcopy.jl")

# Set up keys that enter rebug mode from the regular Julia REPL
# This should be called from your ~/.julia/config/startup.jl file
function repl_init(repl)
    repl.interface = REPL.setup_interface(repl; extra_repl_keymap = get_rebugger_modeswitch_dict())
end

function rebugrepl_init()
    # Set up the Rebugger REPL mode with all of its key bindings
    repl_inited = isdefined(Base, :active_repl)
    while !isdefined(Base, :active_repl)
        sleep(0.05)
    end
    sleep(0.1) # for extra safety
    # Set up the custom "rebug" REPL
    iprompt, eprompt = rebugrepl_init(Base.active_repl, repl_inited)
    interpret_prompt_ref[] = iprompt
    rebug_prompt_ref[] = eprompt
    return nothing
end

function rebugrepl_init(main_repl, repl_inited)
    irepl = HeaderREPL(main_repl, InterpretHeader())
    interface = REPL.setup_interface(irepl; extra_repl_keymap=Dict[])
    iprompt = interface.modes[end]
    erepl = HeaderREPL(main_repl, RebugHeader())
    interface = REPL.setup_interface(erepl; extra_repl_keymap=[get_rebugger_modeswitch_dict(), rebugger_keys])
    eprompt = interface.modes[end]
    add_keybindings(main_repl; override=repl_inited, keybindings...)
    return iprompt, eprompt
end

function __init__()
    # schedule(Task(rebugrepl_init))
    task = Task() do
        try
            rebugrepl_init()
        catch exception
            @error "Rebugger initialization failed" exception=(exception, catch_backtrace())
        end
    end
    schedule(task)
end

include("precompile.jl")
_precompile_()

end # module


text:
# helper functions for mcc172

struct Status
    hardwareoverrun::Bool 		# The board address.
    bufferoverrun::Bool 		# The product ID, one of [HatIDs](@ref HatIDs)
    triggered::Bool 			# The hardware version
    running::Bool				# The product name
end

"""
	mcc172_status_decode(returncode::UInt16)

This function returns a structure of the meaning of the status of the calls: mcc172_a_in_scan_status and mcc172_a_in_scan_read.

	struct Status
		hardwareoverrun::Bool 		# The board address.
		bufferoverrun::Bool 		# The product ID, one of [HatIDs](@ref HatIDs)
		triggered::Bool 			# The hardware version
		running::Bool				# The product name
	end
"""
function mcc172_status_decode(returncode::UInt16)
	# made return code to an Array of descriptive strings
	status = Status(returncode & 0b1 == 0b1 ? true : false,
					returncode & 0b10 == 0b10 ? true : false,
					returncode & 0b100 == 0b100 ? true : false,
					returncode & 0b1000 == 0b1000 ? true : false)
	return status
end

"""
	function chandict()

Dictionary of channel definitions (addresses)
"""
function chandict()
    chdict = Dict{Symbol, UInt8}(
        :CHAN0 => 0x01 << 0,
        :CHAN1 => 0x01 << 1,
        :CHAN2 => 0x01 << 2,
        :CHAN3 => 0x01 << 3,
        :CHAN4 => 0x01 << 4,
        :CHAN5 => 0x01 << 5,
        :CHAN6 => 0x01 << 6,
        :CHAN7 => 0x01 << 7)
end

max_stars_repo_name: Spectrum-Tec/MccDaqHats

max_issues_repo_name: Spectrum-Tec/MccDaqHats

max_forks_repo_name: Spectrum-Tec/MccDaqHats


text:
# ---------------- #
# Fitting routines #
# ---------------- #

# Tensor representation, single function method; calls function that computes coefficients below below
function get_coefs{T}(basis::Basis, bs::BasisMatrix{Tensor}, y::Vector{T})
    _get_coefs_deep(basis,bs,y)[:,1]
end

# Tensor representation, multiple function method; calls function that computes coefficients below
function get_coefs{T}(basis::Basis, bs::BasisMatrix{Tensor}, y::Matrix{T})
    _get_coefs_deep(basis,bs,y)
end

function _get_coefs_deep(basis::Basis,bs::BasisMatrix{Tensor},y)
    if any(bs.order[1, :] .!= 0)
        error("invalid basis structure - first elements must be order 0")
    end
    to_kron = bs.vals[1, :]  # 68
    ckronxi(to_kron, y, ndims(basis):-1:1)  # 66
end

# convert to expanded and call the method below
get_coefs(basis::Basis, bs::BasisMatrix{Direct}, y) =
    get_coefs(basis, convert(Expanded, bs), y)

get_coefs(basis::Basis, bs::BasisMatrix{Expanded}, y) = bs.vals[1] \ y

# common checks to be run at the top of each funfit
function check_funfit(basis::Basis, x, y)
    m = size(y, 1)
    length(basis) > m && error("Can't be more basis funcs than points in y")
    return m
end

# get_coefs(::Basis, ::BasisMatrix, ::Array) does almost all the work for
# these methods
function funfitxy(basis::Basis, bs::BasisMatrix, y)
    m = check_funfit(basis, bs, y)
    c = get_coefs(basis, bs, y)
    c, bs
end

function funfitxy{T}(basis::Basis, x::Vector{Vector{T}}, y)
    m = check_funfit(basis, x, y)

    # additional checks for Array{Any} array
    mm = prod([size(xi, 1) for xi in x])
    mm != m && error("x and y are incompatible")

    # get Tensor form -- most efficient
    bs = BasisMatrix(basis, Tensor(), x, 0)

    # get coefs and return
    c = get_coefs(basis, bs, y)
    c, bs
end

function funfitxy(basis::Basis, x, y)
    # check input sizes
    m = check_funfit(basis, x, y)

    # additional check
    size(x, 1) != m && error("x and y are incompatible")

    bs = BasisMatrix(basis, Direct(), x, 0)
    c = get_coefs(basis, bs, y)
    c, bs
end

function funfitf(basis::Basis, f::Function, args...)
    x = nodes(basis)[1]
    y = f(x, args...)
    funfitxy(basis, x, y)[1]
end

# ---------- #
# Evaluation #
# ---------- #

# funeval wants to evaluate at a matrix. As a stop-gap until I find some
# time, this method makes a scalar x into a 1x1 matrix
funeval(c, basis::Basis{1}, x::Real, order=0) =
    funeval(c, basis, fill(x, 1, 1), order)

# when x is a vector and we have a univariate interpoland, we don't want to
# return just the first element (see next method)
funeval{T<:Number}(c, basis::Basis{1}, x::Vector{T}, order=0) =
    funeval(c, basis, x[:, :], order)

# Here we want only the first element because we have a N>1 dimensional basis
# and we passed only a 1 dimensional array of points at which to evaluate,
# meaning we must have passed a single point.
funeval{N,T<:Number}(c, basis::Basis{N}, x::Vector{T}, order=0) =
    funeval(c, basis, reshape(x, 1, N), order)[1]

function funeval{N,T}(c, basis::Basis{N}, x::Vector{Vector{T}}, order=0)
    # check inputs
    size(x, 1) == N ||  error("x must have d=$N elements")
    order =_check_order(N, order)

    # construct tensor form BasisMatrix
    bs = BasisMatrix(basis, Tensor(), x, order)  # 67

    # pass of to specialized method below
    funeval(c, bs, order)
end

# helper method to construct BasisMatrix{Direct}, then pass to the method
# below below
function funeval{N}(c, basis::Basis{N}, x::Matrix, order=0)
    # check that inputs are conformable
    size(x, 2) == N || error("x must have d=$(N) columns")  # 62
    order =_check_order(N, order)

    # construct BasisMatrix in Direct form
    bs = BasisMatrix(basis, Direct(), x, order)  # 67

    # pass of to specialized method below
    funeval(c, bs, order)
end

# funeval1
function funeval(c, bs::BasisMatrix{Tensor},
                 order::Matrix{Int}=fill(0, 1, size(bs.order, 2)))
    kk, d = size(order)  # 95

    # 98 reverse the order of evaluation: B(d) × B(d-1) × ⋯ × B(1)
    order = flipdim(order .+ (size(bs.vals, 1)*(0:d-1)' - bs.order+1), 2)

    # 99
    nx = prod([size(bs.vals[1, j], 1) for j=1:d])

    f = Array(eltype(c), nx, size(c, 2), kk)  # 100

    for i=1:kk
        f[:, :, i] = ckronx(bs.vals, c, order[i, :])  # 102
    end
    return squeeze_trail(f)
end

function funeval(c, bs::BasisMatrix{Direct},
                 order::Matrix{Int}=fill(0, 1, size(bs.order, 2)))  # funeval2
    kk, d = size(order)  # 95
    # 114 reverse the order of evaluation: B(d)xB(d-1)x...xB(1)
    order = flipdim(order .+ (size(bs.vals, 1)*(0:d-1)' - bs.order+1), 2)

    f = Array(eltype(c), size(bs.vals[1], 1), size(c, 2), kk)  # 116

    for i=1:kk
        f[:, :, i] = cdprodx(bs.vals, c, order[i, :])  # 118
    end
    return squeeze_trail(f)
end

function funeval(c, bs::BasisMatrix{Expanded},
                 order::Matrix{Int}=fill(0, 1, size(bs.order, 2)))  # funeval3
    nx = size(bs.vals[1], 1)
    kk = size(order, 1)
    f = Array(promote_type(eltype(c), eltype(bs.vals[1])), nx, size(c, 2), kk)
    for i=1:kk
        this_order = order[i, :]
        ind = findfirst(x->bs.order[x, :] == this_order, 1:kk)
        if ind == 0
            msg = string("Requested order $(this_order) not in BasisMatrix ",
                         "with order $(bs.order)")
            error(msg)
        end
        f[:, :, i] = bs.vals[ind]*c  # 154
    end

    return squeeze_trail(f)
end

# ------------------------------ #
# Convenience `Interpoland` type #
# ------------------------------ #

type Interpoland{T,N,BST<:ABSR}
    basis::Basis{N}               # the basis -- can't change
    coefs::Vector{T}              # coefficients -- might change
    bmat::BasisMatrix{BST}  # BasisMatrix at nodes of `b`
end

function Interpoland(basis::Basis, bs::BasisMatrix, y::AbstractVector)
    c = get_coefs(basis, bs, y)[:, 1]  # get first (only) column as `Vector`
    Interpoland(basis, c, bs)
end

function Interpoland(basis::Basis, x::Array, y::AbstractVector)
    c, bs = funfitxy(basis, x, y)
    c = c[:, 1]
    Interpoland(basis, c, bs)
end

Interpoland(p::BasisParams, x, y) = Interpoland(Basis(p), x, y)

function Interpoland(basis::Basis, f::Function)
    # TODO: Decide if I want to do this or if I would rather do
    #       x, xd = nodes(basis); y = f(x); Interpoland(basis, xd, y)
    #       to get the BasisMatrix in Tensor format (potentially more
    #       efficient)
    x = nodes(basis)[1]
    y = f(x)
    Interpoland(basis, x, y)
end

Interpoland(p::BasisParams, f::Function) = Interpoland(Basis(p), f)

# let funeval take care of order and such. This just exists to make it so the
# user doesn't have to keep track of the coefficient vector
evaluate(interp::Interpoland, x; order=0) =
    funeval(interp.coefs, interp.basis, x, order)

#=
TODO: I can add an

@stagedfunction evaluate{N,T}(interp::Interpoland{Basis{N}}, x::NTuple{N,Vector{T}})

that calls gridmake for me. Then they can pass points along individual dimensions
=#

# now, given a new vector of `y` data we construct a new coefficient vector
function update_coefs!(interp::Interpoland, y::Vector)
    # leverage the BasisMatrix we kept around
    c = funfitxy(interp.basis, interp.bmat, y)[1]
    copy!(interp.coefs, c)  # update c inplace b/c Interpoland is immutable
end

# similar for a function -- just hand off to above
update_coefs!(interp::Interpoland, f::Function) =
    update_coefs!(interp, f(nodes(interp.basis)[1]))

# alias update_coefs! to fit!
fit!(interp::Interpoland, y::Vector) = update_coefs!(interp, y)
fit!(interp::Interpoland, f::Function) = update_coefs!(interp, f)

Base.show{T,N,BST<:ABSR}(io::IO, ::Interpoland{T,N,BST}) =
    print(io, "$N dimensional interpoland with $BST BasisMatrix")


text:
using Polyhedra

"""
Test `center_of_mass` in a case where the center of mass is different from the
centroid (namely, a square pyramid).
"""
function comsquarepyramidtest(lib::Polyhedra.Library)
    poly = polyhedron(convexhull([0, 0, 0],
                                 [1, 0, 0],
                                 [0, 1, 0],
                                 [1, 1, 0],
                                 [1/2, 1/2, 1]),
                      lib)
    # Centroid is [1/2, 1/2, 1/5].
    # Center of mass is [1/2, 1/2, 1/4].
    @test center_of_mass(poly) ≈ [1/2, 1/2, 1/4]
end

max_stars_repo_name: blegat/Polyhedra

max_issues_repo_name: blegat/Polyhedra

max_forks_repo_name: blegat/Polyhedra


text:
"""
    TypeZ()

Lorem ipsum dolor sit amet.
"""
struct TypeZ

end


text:
using SymbolicWedderburn
using PermutationGroups
using Cyclotomics

using DynamicPolynomials
using SumOfSquares
using SCS

@polyvar x y z

motzkin = x^4*y^2 + y^4*x^2 - 3*x^2*y^2 + 1

m = let m = SOSModel(optimizer_with_attributes(SCS.Optimizer, "eps"=>1e-5, "acceleration_lookback"=>10))

    @variable m t
    @objective m Max t
    @variable m sos SOSPoly(monomials([x, y], 0:5))
    @constraint m motzkin - t == sos
    optimize!(m)
    m
end

using DynamicPolynomials
using MultivariatePolynomials
const MP = MultivariatePolynomials
using MultivariateBases

function SymbolicWedderburn.OnPoints(basis::Union{<:MonomialVector, AbstractVector{<:Monomial}})
    basis_exps = Vector{Vector{Int}}(undef, length(basis))
    basis_dict = Dict{Vector{Int}, Int}()
    sizehint!(basis_dict, length(basis))

    for (i, b) in enumerate(basis)
        e = MP.exponents(b) # so that we allocate exponents only once
        basis_exps[i] = e
        basis_dict[e] = i
    end

    return SymbolicWedderburn.OnPoints(basis_exps, basis_dict)
end

include(joinpath(@__DIR__, "..", "test", "smallgroups.jl"));
G = PermGroup([perm"(1,2)"])
chars_vars = SymbolicWedderburn.characters_dixon(G)

mvec = monomials([x,y], 0:3)

chars_mvec = let chars = chars_vars, basis = mvec

    @assert all(χ.inv_of == first(chars).inv_of for χ in chars)

    induced_action = SymbolicWedderburn.OnPoints(basis)
    ccG_large = induced_action.(conjugacy_classes(first(chars)))

    # double check:
    let ccls = ccG_large, large_gens = induced_action.(gens(G))
        G_large = PermGroup(large_gens)
        ccG_large = conjugacy_classes(G_large)
        @assert all(Set.(collect.(ccG_large)) .== Set.(collect.(ccls)))
    end

    [SymbolicWedderburn.Character(values(χ), χ.inv_of, ccG_large) for χ in chars]
end

vr_chars = SymbolicWedderburn.real_vchars(chars_mvec)

U = filter!(x->all(!iszero, x), [SymbolicWedderburn.matrix_projection(χ) for χ in vr_chars])

R = map(U) do c_u
    u = last(c_u)
    if all(isreal, u)
        image_coeffs, pivots = SymbolicWedderburn.row_echelon_form(float.(u))
        dim = length(pivots)
        image_coeffs[1:dim, :]
    else
        throw("Not Implemented")
    end
end

R = filter!(!iszero ∘ first ∘ size, R)

msym = let msym = SOSModel(optimizer_with_attributes(SCS.Optimizer,
    "eps"=>1e-12,
    "max_iters"=>100_000,
    "acceleration_lookback"=>0,
    "alpha"=>1.95,
    ))

    @variable   msym t
    @objective  msym Max t

    @variable   msym sos1 SOSPoly(FixedPolynomialBasis(R[1]*mvec))
    @variable   msym sos2 SOSPoly(FixedPolynomialBasis(R[2]*mvec))

    @variable   msym sos3 SOSPoly(monomials([x,y], 0:1))

    @constraint msym motzkin - t + 0.125 == sos1 + sos2 + (2 - x^2 - y^2)*sos3
    @constraint msym t <= -0.05
    optimize!(msym)
    msym
end

termination_status(msym)


text:
export make_graph

"""
    make_graph(A::Matrix, verbose::Bool=false)

Use the matrix `A` to build a graph. Ideally, `A` is an adjacency matrix.
But if not we (1) zero out the diagonal, (2) make it symmetric, and (3)
make it zero-one.

Optional second argument `verbose=true` reports what's happening.
"""
function make_graph(A::Matrix, verbose::Bool = false)::SimpleGraph
    n, c = size(A)
    if n != c
        error("Matrix must be square; can't fix that.")
    end

    A = copy(A)

    verbose && println("Making sure the matrix is symmetric")
    B = transpose(A)
    if A != B
        A += B
    end

    verbose && println("Making sure the diagonal is all zeros")
    for j = 1:n
        @inbounds A[j, j] = 0
    end

    verbose && println("Making sure the entries are zeros and ones")
    A = Int.(A .!= 0)

    verbose && println("Creating the graph")
    return SimpleGraph(A)

end


text:
"""
    PWM(i::Int32)
This device allows for low level PWM control of selected pins. The valid pins
dictionary pwm_pins relates to memory adresses in of the AM3359 chip, see p.182
in www.ti.com/product/AM3359/technicaldocuments.
"""

type PWM <: IO_Object
    i::Int32
    pin::String
    chip::String
    basedir::String
    filestreams::Array{IOStream,1}
    function PWM(i::Int32)

        (i < 1 || i > length(pwm_pins)) && error("Invalid PWM index: $i")

        # Export the PWM pin
        basedir = export_pwm(i)

        # Setup filestreams
        pins = collect(keys(pwm_pins))
        pin = pins[i]
        chip = pwm_pins[pin][2]

        enable_filestream = open("$(basedir)/$(pwm_pins[pin][2])/pwm$(pwm_pins[pin][3])/enable","r+")
        period_filestream = open("$(basedir)/$(pwm_pins[pin][2])/pwm$(pwm_pins[pin][3])/period","r+")
        duty_cycle_filestream = open("$(basedir)/$(pwm_pins[pin][2])/pwm$(pwm_pins[pin][3])/duty_cycle","r+")
        polarity_filestream = open("$(basedir)/$(pwm_pins[pin][2])/pwm$(pwm_pins[pin][3])/polarity","r+")
        return new(i, pin, chip, basedir, [enable_filestream, period_filestream, duty_cycle_filestream, polarity_filestream])
    end
end

"""
    write!(pwm::PWM, args::Tuple{Int32,String}, debug::Bool=false)
Writes an entry to an operation on the PWM, of the form args = (operation, entry).
"""
function write!(pwm::PWM, args::Tuple{Int32,String}, debug::Bool=false)
    debug && return

    operation, entry = args[1], args[2]
    (operation < 1 || operation > length(pwm.filestreams)) && error("Invalid PWM operation: $operation")

    # Input data check
    assert_pwm_write(operation, entry)

    # Write to file
    seekstart(pwm.filestreams[operation])
    write(pwm.filestreams[operation], "$entry\n")
    flush(pwm.filestreams[operation])
end

"""
    assert_pwm_write(operation::Int32, entry::String)
Assertsion for the PWM input data.
"""
function assert_pwm_write(operation::Int32, entry::String)
  if operation == "1"
    entry ∉ ["0", "1"] && error("Invalid SysLED entry $(entry), valid options are 0 and 1 ::String")
  else
    number = try
      parse(Int32, entry)
    catch
      error("Invalid SysLED entry $(entry), cannot parse as Int32")
    end
    (number < 0 || number > 100000000) && error("Invalid SysLED entry $(entry), not in the range [0,100000000]")
  end
end

"""
    l = read(pwm::PWM, operation::Int32, debug::Bool=false)
Reads the current value from an operation on a GPIO.
"""
function read(pwm::PWM, operation::Int32, debug::Bool=false)
  debug && return
  # Filestreams 1, 2 and 3 are readable
  operation ∉ [1,2,3,4] && error("Invalid GPIO operation: $operation for reading")
  seekstart(pwm.filestreams[operation])
  l = readline(pwm.filestreams[operation])
  return l
end

"""
    teardown!(pwd::PWM)
Closes all open streams on the PWM, and unexports it from the file system
"""
function teardown(pwm::PWM, debug::Bool=false)
  debug && return

  #Close all IOStreams
  for stream in pwm.filestreams
    close(stream)
  end

  if isdefined(:RUNNING_TESTS)
    # Remove the dummy file system for testing
    try
      rm("$(pwm.basedir)/$(pwm_pins[pwm.pin][2])/pwm$(pwm_pins[pwm.pin][3])"; recursive=true)
    catch
      error("Could not remove the requested GPIO testfiles for channel $(pwm_pins[pwm.pin][2])/pwm$(pwm_pins[pwm.pin][3]).")
    end
  else
    #Unexport filestructure
    filename = "/sys/class/pwm/$(pwm.chip)/unexport"
    export_number = pwm_pins[pwm.pin][3]
    write(filename, export_number)
  end
end

"""
    export_gpio(i::Int32, debug::Bool=false)
Export the GPIO file system, either for real-time or testing usecases.
"""
function export_pwm(i::Int32)
  # Find chip and export number
  pins = collect(keys(pwm_pins))
  pin = pins[i]
  chip = pwm_pins[pin][2]

  if isdefined(:RUNNING_TESTS)
    # Export a dummy file system for testing
    basedir = "$(pwd())/testfilesystem/pwm"
    complete_path = "$(basedir)/$(pwm_pins[pin][2])/pwm$(pwm_pins[pin][3])"
    try
      mkpath(complete_path)
    catch
      error("Could not export the PWM device for $(pwm_pins[pin][2]) for testing as the directory $(basedir)/$(pwm_pins[pin][2])/pwm$(pwm_pins[pin][3]) already exists.")
    end
    try
      f = open("$(complete_path)/enable", "w"); write(f,"0"); close(f);
      f = open("$(complete_path)/period", "w"); write(f,"0"); close(f);
      f = open("$(complete_path)/duty_cycle", "w"); write(f,"0"); close(f);
      f = open("$(complete_path)/polarity", "w"); write(f,"0"); close(f);
    catch
      error("Could not open the requested GPIO testfiles for $(complete_path).")
    end
  else
    basedir = "/sys/class/pwm"

    # Configure the pin to run PWM if possible
    Base.run(`config-pin $(pin) pwm`)

    # Export the filestructure of the corresponding chip
    filename = "/sys/class/pwm/$(chip)/export"
    exportNumber = pwm_pins[pin][3]
    write(filename, exportNumber)
  end
  return basedir
end

"""
    to_string(pwm::PWM,, debug::Bool=false)
Generates a string representation of the GPIO device.
"""
function to_string(pwm::PWM, debug::Bool=false)
  debug && return
  message = "\nID: $(pwm.i)\n\nAvailable filestreams:\n"
  for ii = 1:length(pwm.filestreams)
    message = string(message, "  index=$(ii) - name=$(pwm.filestreams[ii].name) - write/read=$(iswritable(pwm.filestreams[ii]))/$(isreadable(pwm.filestreams[ii]))\n")
  end
  return message
end


text:
function strength(g::MetaDiGraph, v::Integer; normalize::Bool=true)
    instrength(g, v; normalize=normalize) + outstrength(g, v; normalize=normalize)
end
function strength(g::MetaGraph, v::Integer; normalize::Bool=true)
    instrength(g, v; normalize=normalize)
end

function strength(g, vs::AbstractVector=vertices(g); normalize::Bool=true)
    [strength(g, v_i; normalize=normalize) for v_i in vs]
end

function instrength(g, vs::AbstractVector=vertices(g); normalize::Bool=true)
    [instrength(g, v_i; normalize=normalize) for v_i in vs]
end
function instrength(g::AbstractMetaGraph{T,M}, v::Integer; normalize::Bool=true) where {T,M}
    v_v = inneighbors(g, v)
    out = zero(weighttype(M))
    for e in getfield(v_v, :v)
        out += metadata(e, :weight)
    end
    if normalize
        return out / (length(v_v) - 1)
    else
        return out
    end
end

function outstrength(g, vs::AbstractVector=vertices(g); normalize::Bool=true)
    [outstrength(g, v_i; normalize=normalize) for v_i in vs]
end
function outstrength(g::AbstractMetaGraph{T,M}, v::Integer; normalize::Bool=true) where {T,M}
    v_v = outneighbors(g, v)
    out = zero(weighttype(M))
    for e in getfield(v_v, :v)
        out += metadata(e, :weight)
    end
    if normalize
        return out / (length(v_v) - 1)
    else
        return out
    end
end


text:
### A Pluto.jl notebook ###
# v0.12.20

using Markdown
using InteractiveUtils

# ╔═╡ 92fa7440-6b38-11eb-0f13-dd8cc31baaee
using GMT

# ╔═╡ ca2e6c40-6b39-11eb-2416-5f9501d31726
md"# Bar plots

A simple bar plot showing color and bar width (in data units) assignement.
"

# ╔═╡ e6c6d3b2-6b39-11eb-0392-593da0f77eb6
bar(1:5, (20, 35, 30, 35, 27), width=0.5, color=:lightblue, limits=(0.5,5.5,0,40), show=true)

# ╔═╡ d7ddbbb0-6bf3-11eb-233e-3f2155a88415
md"### -

A colored bar plot with colors proportional to bar heights. In this case we let the plot limits be determined from data. We also plot a colorbar by using the **colorbar**=*true* option."

# ╔═╡ 099787d0-6bf4-11eb-0a54-e3133dd45400
bar(rand(15), color=:turbo, figsize=(14,8), title="Colored bars", colorbar=true, show=true)

# ╔═╡ 524f7d00-6bf6-11eb-0d30-531720224196
md"### -
Example showing how to plot bar groups and at same time assign variable transparency to each of the group's *band* usinf the **fillalpha** option. Here, each row on the input data represents a bar group that has as many *bands* as n_columns - 1. -1 because first column must hold the *xx* coordinates of each group. The colors come from the automatic cyclic scheme.
"

# ╔═╡ 3e636a70-6b3a-11eb-183c-7fcb8a761de1
bar([0. 1 2 3; 1 2 3 4], fillalpha=[0.3 0.5 0.7], show=true)

# ╔═╡ 6cb27160-6bf7-11eb-029a-0b54adf510fe
md"### -
Next example shows how to plot error bars in a grouped bar. Similar to this [mapplotlib's example](https://matplotlib.org/3.1.1/gallery/lines_bars_and_markers/barchart.html#sphx-glr-gallery-lines-bars-and-markers-barchart-py) (labels will come later).
"

# ╔═╡ 2c820e2e-6b3d-11eb-39aa-630c8d82953b
bar(1:5, [20 25; 35 32; 30 34; 35 20; 27 25], width=0.7, fill=["lightblue", "brown"],
    error_bars=(y=[2 3; 3 5; 4 2; 1 3; 2 3],), xticks=(:G1, :G2, :G3, :G4, :G5), yaxis=(annot=5,label=:Scores), frame=(title="Scores by group and gender", axes=:WSrt), show=true)

# ╔═╡ 43b175e2-6b3e-11eb-141e-71bb8c1253c2
md"### -
Example of a verticaly stacked bar plot. In this exampled we pass the *xx* coordinates as first argument and the individual bar heights in a matrix with smae number of rows as the number of elements in the *x* vector. To make it plot a stracked bar we used the option **stacked**=*true*.
"

# ╔═╡ 52167540-6b3e-11eb-343a-596dd9483996
bar(1:3,[-5 -15 20; 17 10 21; 10 5 15], stacked=1, show=1)

# ╔═╡ 92f28e50-6be3-11eb-2e4a-5798e1617080
md"To create an horizontal bar plot we use the **hbar**=*true* option"

# ╔═╡ bb6ca5a0-6be3-11eb-268d-d510fa533101
bar([0. 1 2 3; 1 2 3 4], hbar=true, show=true, name="bars5.png")

# ╔═╡ d47d45de-6be3-11eb-12c0-edfff29dced9
md"### -
And one horizontal and stacked but this time we pick the colors.
"

# ╔═╡ 9a94db40-6bed-11eb-35b6-b1ac3148df98
bar([0. 1 2 3; 1 2 3 4], stack=true, hbar=true, fill=["red", "green", "blue"], show=true)

# ╔═╡ Cell order:
# ╟─ca2e6c40-6b39-11eb-2416-5f9501d31726
# ╠═92fa7440-6b38-11eb-0f13-dd8cc31baaee
# ╠═e6c6d3b2-6b39-11eb-0392-593da0f77eb6
# ╟─d7ddbbb0-6bf3-11eb-233e-3f2155a88415
# ╠═099787d0-6bf4-11eb-0a54-e3133dd45400
# ╟─524f7d00-6bf6-11eb-0d30-531720224196
# ╠═3e636a70-6b3a-11eb-183c-7fcb8a761de1
# ╟─6cb27160-6bf7-11eb-029a-0b54adf510fe
# ╠═2c820e2e-6b3d-11eb-39aa-630c8d82953b
# ╟─43b175e2-6b3e-11eb-141e-71bb8c1253c2
# ╠═52167540-6b3e-11eb-343a-596dd9483996
# ╟─92f28e50-6be3-11eb-2e4a-5798e1617080
# ╠═bb6ca5a0-6be3-11eb-268d-d510fa533101
# ╟─d47d45de-6be3-11eb-12c0-edfff29dced9
# ╠═9a94db40-6bed-11eb-35b6-b1ac3148df98


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

include("mh_sampler.jl")
include("mh_tuner.jl")


text:
module test_octo_backends

using Test # @test_throws
using Octo.Backends

module UnsupportedAdapter
end

@test_throws Backends.UnsupportedError Backends.backend(UnsupportedAdapter)

end # module test_octo_backends


text:
import TurbulenceConvection
import CLIMAParameters
const CP = CLIMAParameters

struct EarthParameterSet{NT} <: CP.AbstractEarthParameterSet
    nt::NT
end

CLIMAParameters.Planet.MSLP(ps::EarthParameterSet) = ps.nt.MSLP
CLIMAParameters.Planet.cp_d(ps::EarthParameterSet) = ps.nt.cp_d
CLIMAParameters.Planet.cp_v(ps::EarthParameterSet) = ps.nt.cp_v
CLIMAParameters.Planet.R_d(ps::EarthParameterSet) = ps.nt.R_d
CLIMAParameters.Planet.R_v(ps::EarthParameterSet) = ps.nt.R_v
CLIMAParameters.Planet.molmass_ratio(ps::EarthParameterSet) = ps.nt.molmass_ratio
# microphysics 0-moment parameters
CLIMAParameters.Atmos.Microphysics_0M.τ_precip(ps::EarthParameterSet) = ps.nt.τ_precip
# microphysics 1-moment parameters
CLIMAParameters.Atmos.Microphysics.τ_acnv_rai(ps::EarthParameterSet) = ps.nt.τ_acnv_rai
# entrainment/detrainment parameters
CLIMAParameters.Atmos.EDMF.c_ε(ps::EarthParameterSet) = ps.nt.c_ε # factor multiplyer for dry term in entrainment/detrainment
CLIMAParameters.Atmos.EDMF.α_b(ps::EarthParameterSet) = ps.nt.α_b # factor multiplyer for pressure buoyancy terms (effective buoyancy is (1-α_b))
CLIMAParameters.Atmos.EDMF.α_a(ps::EarthParameterSet) = ps.nt.α_a # factor multiplyer for pressure advection
CLIMAParameters.Atmos.EDMF.α_d(ps::EarthParameterSet) = ps.nt.α_d # factor multiplyer for pressure drag
CLIMAParameters.Atmos.EDMF.H_up_min(ps::EarthParameterSet) = ps.nt.H_up_min # minimum updarft top to avoid zero devision in pressure drag and turb-entr
CLIMAParameters.Atmos.EDMF.ω_pr(ps::EarthParameterSet) = ps.nt.ω_pr # experiment scale factor for turbulent Prandtl number
CLIMAParameters.Atmos.EDMF.c_δ(ps::EarthParameterSet) = ps.nt.c_δ # factor multiplyer for moist term in entrainment/detrainment
CLIMAParameters.Atmos.EDMF.β(ps::EarthParameterSet) = ps.nt.β # sorting power for ad-hoc moisture detrainment function
CLIMAParameters.Atmos.EDMF.χ(ps::EarthParameterSet) = ps.nt.χ # fraction of updraft air for buoyancy mixing in entrainment/detrainment (0≤χ≤1)
CLIMAParameters.Atmos.EDMF.c_t(ps::EarthParameterSet) = ps.nt.c_t # factor multiplyer for turbulent term in entrainment/detrainment
CLIMAParameters.Atmos.EDMF.c_λ(ps::EarthParameterSet) = ps.nt.c_λ # scaling factor for TKE in entrainment scale calculations
CLIMAParameters.Atmos.EDMF.w_min(ps::EarthParameterSet) = ps.nt.w_min # minimum updraft velocity to aviod zero division in b/w²
CLIMAParameters.Atmos.EDMF.μ_0(ps::EarthParameterSet) = ps.nt.μ_0 # dimentional scale logistic function in the dry term in entrainment/detrainment
# mixing length parameters
CLIMAParameters.Atmos.EDMF.c_m(ps::EarthParameterSet) = ps.nt.c_m # tke diffusivity coefficient
CLIMAParameters.Atmos.EDMF.c_d(ps::EarthParameterSet) = ps.nt.c_d # tke dissipation coefficient
CLIMAParameters.Atmos.EDMF.Pr_n(ps::EarthParameterSet) = ps.nt.Pr_n # tke dissipation coefficient
CLIMAParameters.Atmos.EDMF.smin_ub(ps::EarthParameterSet) = ps.nt.smin_ub #  lower limit for smin function
CLIMAParameters.Atmos.EDMF.smin_rm(ps::EarthParameterSet) = ps.nt.smin_rm #  upper ratio limit for smin function

#! format: off
function create_parameter_set(namelist)
    TC = TurbulenceConvection
    nt = (;
        MSLP = 100000.0, # or grab from, e.g., namelist[""][...]
        cp_d = 1004.0,
        cp_v = 1859.0,
        R_d = 287.1,
        R_v = 461.5,
        molmass_ratio = 461.5/287.1,
        τ_precip = TC.parse_namelist(namelist, "microphysics", "τ_precip"; default = 1000.0),
        τ_acnv_rai = TC.parse_namelist(namelist, "microphysics", "τ_acnv_rai"; default = 1000.0),
        c_ε = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "entrainment_factor"),
        c_div = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "entrainment_massflux_div_factor"; default = 0.0),
        α_b = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "pressure_normalmode_buoy_coeff1"),
        α_a = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "pressure_normalmode_adv_coeff"),
        α_d = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "pressure_normalmode_drag_coeff"),
        H_up_min = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "min_updraft_top"),
        ω_pr = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "Prandtl_number_scale"),
        c_δ = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "detrainment_factor"),
        c_gen = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "general_ent_params"),
        β = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "sorting_power"),
        χ = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "updraft_mixing_frac"),
        c_t = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "turbulent_entrainment_factor"),
        c_λ = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "entrainment_smin_tke_coeff"),
        w_min = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "min_upd_velocity"),
        μ_0 = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "entrainment_scale"),
        μ = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "entrainment_sigma"),
        γ_lim = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "area_limiter_scale"),
        β_lim = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "area_limiter_power"),
        c_m = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "tke_ed_coeff"),
        c_d = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "tke_diss_coeff"),
        c_b = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "static_stab_coeff"; default = 0.4), # this is here due to an value error in CliMAParmaeters.jl
        Pr_n = TC.parse_namelist(namelist, "turbulence", "prandtl_number_0"), # this is here due to an value error in CliMAParmaeters.jl
        smin_ub = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "smin_ub"),
        smin_rm = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "smin_rm"),
        l_max = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "l_max"; default = 1.0e6),

        ## Stochastic parameters
        # lognormal model
        stoch_ε_lognormal_var = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "stochastic", "entr_lognormal_var"; default = 0.0),
        stoch_δ_lognormal_var = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "stochastic", "detr_lognormal_var"; default = 0.0),
        # sde model
        sde_ϵ_θ = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "stochastic", "sde_entr_theta"; default = 1.0),
        sde_ϵ_σ = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "stochastic", "sde_entr_std"; default = 0.0),
        sde_δ_θ = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "stochastic", "sde_detr_theta"; default = 1.0),
        sde_δ_σ = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "stochastic", "sde_detr_std"; default = 0.0),
    )
    param_set = EarthParameterSet(nt)
    if !isbits(param_set)
        @show param_set
        error("The parameter set MUST be isbits in order to be stack-allocated.")
    end
    return param_set
end
#! format: on


text:
if isdefined(FileIO, :action)
    # FileIO >= 1.6
    _format_stream(format, io) = FileIO.Stream{format}(io)
else
    _format_stream(format, io) = FileIO.Stream(format, io)
end

if VERSION < v"1.2.0"
    isnothing(x) = x === nothing
end


text:
bitmap(grid::Grid2D) = Bitmap2D(grid)
bitmap(grid::Grid3D) = Bitmap3D(grid)

empty(M::Bitmap2D) = Bitmap2D(M.grid)
empty(M::Bitmap3D) = Bitmap3D(M.grid)
empty(grid::Grid2D) = Bitmap2D(grid)
empty(grid::Grid3D) = Bitmap3D(grid)

grid(B::AbstractBitmap) = B.grid


text:
using Documenter, Latte

makedocs(
    modules = [Latte],
    doctest  = false
)

custom_deps() = run(`pip install --user pygments mkdocs mkdocs-material`)

deploydocs(
    deps = custom_deps,
    repo = "github.com/IntelLabs/Latte.jl.git",
    julia = "0.4"
)


text:
##
##  f z e r o . j l  Root Finding
##

function ridders(f::Function, a::Real, b::Real;
                 maxiter::Integer = 1000, tol::Real = eps())

    x1 = a;     x2 = b
    f1 = f(x1); f2 = f(x2)
    if f1 * f2 > 0
        error("f(a) and f(b) must have different signs.")
    elseif f1 == 0
        return x1
    elseif f2 == 0
        return f2
    end

    niter = 2
    while niter < maxiter && abs(x1 - x2) > tol
        xm = (x1 + x2)/2.0
        fm = f(xm)
        if fm == 0; return xm; end

        x3 = xm + (xm - x1) * sign(f1 - f2) * fm / sqrt(fm^2 - f1 * f2)
        f3 = f(x3)
        niter += 2
        if f3 == 0; return x3; end

        if fm * f3 < 0
            x1 = xm;  f1 = fm
            x2 = x3;  f2 = f3
        elseif f1 * f3 < 0
            x2 = x3;  f2 = f3
        elseif f2 * f3 < 0
            x1 = x3;  f1 = f3
        else
            error("Inform the maintainer: you should never get here.")
        end
    end

    return (x1 + x2) / 2.0
end

function brent_dekker(f::Function, a::Real, b::Real;
                       maxiter::Integer = 1000, tol::Real = eps())

	x1 = a; f1 = f(x1)
	if f1 == 0; return a; end
	x2 = b; f2 = f(x2)
	if f2 == 0; return b; end
	if (f1*f2 > 0.0)
	    error("Brent-Dekker: Root is not bracketed in [a, b]")
	end

	x3 = 0.5*(a+b)
	# Beginning of iterative loop
	niter = 1
	while (niter <= maxiter)
		f3 = f(x3)
		if abs(f3) < tol
		    x0 = x3
		    break
		end

		# Tighten brackets [a, b] on the root
		if f1*f3 < 0.0
		    b = x3
		else
		    a = x3
		end
		if (b-a) < tol*max(abs(b), 1.0)
		    x0 = 0.5*(a + b)
		    break
	    end

		# Try quadratic interpolation
		denom = (f2 - f1)*(f3 - f1)*(f2 - f3)
		numer = x3*(f1 - f2)*(f2 - f3 + f1) + f2*x1*(f2 - f3) + f1*x2*(f3 - f1)
		# if denom==0, push x out of bracket to force bisection
		if denom == 0
			dx = b - a
		else
			dx = f3*numer/denom
		end

		x = x3 + dx
		# If interpolation goes out of bracket, use bisection.
		if (b - x)*(x - a) < 0.0
			dx = 0.5*(b - a)
			x  = a + dx;
		end

		# Let x3 <-- x & choose new x1, x2 so that x1 < x3 < x2.
		if x1 < x3
			x2 = x3; f2 = f3
		else
			x1 = x3; f1 = f3
		end

		niter += 1
		if abs(x - x3) < tol
		    x0 = x
		    break
	    end
		x3 = x
	end

    if (niter > maxiter)
        println("Maximum number of iterations has been reached.")
    end

    prec = min(abs(x1-x3), abs(x2-x3))
    return x3
end


text:
function enqueue_transmissions!(state::SimState, ::Val{ConstantKernelContact}, source_id::Integer, params::SimParams)
  progression = progressionof(params, source_id)

  start_time = progression.incubation_time
  end_time =  if      !ismissing(progression.mild_symptoms_time);   progression.mild_symptoms_time
              elseif  !ismissing(progression.severe_symptoms_time); progression.severe_symptoms_time
              elseif  !ismissing(progression.recovery_time);        progression.recovery_time
              else    error("no recovery nor symptoms time defined")
              end

  strain = strainof(state, source_id)

  total_infection_rate = (end_time - start_time) * params.constant_kernel_param
  total_infection_rate *= straindata(params, strain).outside_multiplier
  total_infection_rate *= spreading(params, source_id)

  num_infections = rand(state.rng, Poisson(total_infection_rate))

  if num_infections == 0
    return
  end
  @assert start_time != end_time "pathologicaly short time for infections there shouldn't be any infections but are $num_infections, progression=$progression"

  time_dist = Uniform(state.time, end_time - start_time + state.time) # in global time reference frame

  num_individuals = numindividuals(params)

  for _ in 1:num_infections
    subject_id = sample(state.rng, 1:(num_individuals-1)) # exclude the source itself
    if subject_id >= source_id # restore the indexing
      subject_id +=1
    end

    if Healthy == health(state, subject_id)
      infection_time::TimePoint = rand(state.rng, time_dist) |> TimePoint
      @assert state.time <= infection_time <= (end_time-start_time + state.time)
      push!(state.queue, Event(Val(TransmissionEvent), infection_time, subject_id, source_id, ConstantKernelContact, strain))
    end
  end
  nothing
end

function enqueue_transmissions!(state::SimState, ::Val{HouseholdContact}, source_id::Integer, params::SimParams)
  household = householdof(params, source_id)

  if 1 == length(household)
    return
  end

  progression = progressionof(params, source_id)

  start_time = progression.incubation_time
  end_time =  if      !ismissing(progression.severe_symptoms_time); progression.severe_symptoms_time
              elseif  !ismissing(progression.recovery_time);        progression.recovery_time
              else    error("no recovery nor severe symptoms time defined")
              end

  max_time = time(state) - start_time + end_time

  strain = strainof(state, source_id)
  mean_infection_time = (length(household)-1) / params.household_kernel_param
  mean_infection_time /= straindata(params, strain).household_multiplier
  time_dist = Exponential(mean_infection_time)

  for subject_id in household
    if subject_id == source_id || Healthy != health(state, subject_id)
      continue
    end

    infection_time = time(state) + rand(state.rng, time_dist)
    if infection_time > max_time
      continue
    end

    @assert time(state) <= infection_time <= (end_time - start_time + time(state))
    push!(state.queue, Event(Val(TransmissionEvent),
      infection_time,
      subject_id,
      source_id,
      HouseholdContact,
      strain)
    )
  end
  nothing
end

function enqueue_transmissions!(state::SimState, ::Val{HospitalContact}, source_id::Integer, params::SimParams)
  if nothing === params.hospital_kernel_params
    return
  end
  progression = progressionof(params, source_id)

  @assert progression.severity in SA[Severe, Critical]
  @assert !ismissing(progression.severe_symptoms_time)
  @assert !ismissing(progression.recovery_time) || !ismissing(progression.death_time)

  start_time = progression.severe_symptoms_time
  end_time = ismissing(progression.recovery_time) ? progression.death_time : progression.recovery_time

  strain = strainof(state, source_id)

  total_infection_rate = (end_time - start_time) * params.hospital_kernel_params.kernel_constant
  total_infection_rate *= straindata(params, strain).outside_multiplier

  num_infections = rand(state.rng, Poisson(total_infection_rate))

  if num_infections == 0
    return
  end
  @assert start_time != end_time "pathologicaly short time for infections there shouldn't be any infections but are $num_infections, progression=$progression"

  time_dist = Uniform(state.time, end_time - start_time + state.time) # in global time reference frame

  for _ in 1:num_infections
    subject_id = sample(state.rng, params.hospital_kernel_params.hospital_staff_ids)
    if subject_id == source_id # self infection not possible
      continue
    end

    if Healthy == health(state, subject_id)
      infection_time::TimePoint = rand(state.rng, time_dist) |> TimePoint
      @assert state.time <= infection_time <= (end_time-start_time + state.time)
      push!(state.queue, Event(Val(TransmissionEvent), infection_time, subject_id, source_id, HospitalContact))
    end
  end
end

function enqueue_transmissions!(state::SimState, ::Val{AgeCouplingContact}, source_id::Integer, params::SimParams)
  if params.age_coupling_params === nothing
    return
  end

  progression = progressionof(params, source_id)

  start_time = progression.incubation_time
  end_time =  if      !ismissing(progression.mild_symptoms_time);   progression.mild_symptoms_time
              elseif  !ismissing(progression.severe_symptoms_time); progression.severe_symptoms_time
              elseif  !ismissing(progression.recovery_time);        progression.recovery_time
              else    error("no recovery nor symptoms time defined")
              end

  strain = strainof(state, source_id)

  total_infection_rate = (end_time - start_time) * straindata(params, strain).outside_multiplier
  total_infection_rate *= sourceweightof(params.age_coupling_params, source_id)
  total_infection_rate *= spreading(params, source_id)

  num_infections = rand(state.rng, Poisson(total_infection_rate))

  if num_infections == 0
    return
  end
  @assert start_time != end_time "pathologicaly short time for infections there shouldn't be any infections but there are $num_infections, progression=$progression"

  time_dist = Uniform(state.time, end_time - start_time + state.time) # in global time reference frame

  for _ in 1:num_infections
    subject_id = sample(state.rng, params.age_coupling_params.coupling, source_id)

    if subject_id == source_id # self-infection difficult to avoid at sampling
      break
    end

    if Healthy == health(state, subject_id)
      infection_time::TimePoint = rand(state.rng, time_dist) |> TimePoint
      @assert state.time <= infection_time <= (end_time-start_time + state.time)
      push!(state.queue, Event(Val(TransmissionEvent), infection_time, subject_id, source_id, AgeCouplingContact, strain))
    end
  end
  nothing
end