text:
function longest_common_prefix(prefixa, prefixb)
    maxplength = min(sizeof(prefixa), sizeof(prefixb))
    maxplength == 0 && return ""
    idx = findfirst(i -> (prefixa[i] != prefixb[i]), 1:maxplength)
    idx = idx == nothing ? maxplength : idx - 1
    prefixa[1:idx]
end

function skip_to_nl(str, idxend)
    while (idxend < sizeof(str)) && str[idxend] != '\n'
        idxend = nextind(str, idxend)
    end
    idxend > sizeof(str) ? prevind(str, idxend) : idxend
end

tostr(buf::IOBuffer) = _unescape_string(String(take!(buf)))

"""
parse_string_or_cmd(ps)

When trying to make an `INSTANCE` from a string token we must check for
interpolating opoerators.
"""
function parse_string_or_cmd(ps::ParseState, prefixed = false)
    sfullspan = ps.nt.startbyte - ps.t.startbyte
    sspan = 1 + ps.t.endbyte - ps.t.startbyte

    istrip = (ps.t.kind == Tokens.TRIPLE_STRING) || (ps.t.kind == Tokens.TRIPLE_CMD)
    iscmd = ps.t.kind == Tokens.CMD || ps.t.kind == Tokens.TRIPLE_CMD

    if ps.errored
        return ErrorToken()
    end

    lcp = nothing
    exprs_to_adjust = []
    function adjust_lcp(expr, last = false)
    end
    function adjust_lcp(expr::LITERAL, last = false)
        push!(exprs_to_adjust, expr)
        str = expr.val
        (isempty(str) || (lcp != nothing && isempty(lcp))) && return
        (last && str[end] == '\n') && return (lcp = "")
        idxstart, idxend = 2, 1
        while nextind(str, idxend) - 1 < sizeof(str) && (lcp == nothing || !isempty(lcp))
            idxend = skip_to_nl(str, idxend)
            idxstart = nextind(str, idxend)
            while nextind(str, idxend) - 1 < sizeof(str)
                c = str[nextind(str, idxend)]
                if c == ' ' || c == '\t'
                    idxend += 1
                elseif c == '\n'
                    # All whitespace lines in the middle are ignored
                    idxend += 1
                    idxstart = idxend + 1
                else
                    prefix = str[idxstart:idxend]
                    lcp = lcp === nothing ? prefix : longest_common_prefix(lcp, prefix)
                    break
                end
            end
        end
        if idxstart != nextind(str, idxend)
            prefix = str[idxstart:idxend]
            lcp = lcp === nothing ? prefix : longest_common_prefix(lcp, prefix)
        end
    end

    # there are interpolations in the string
    if prefixed != false || iscmd
        t_str = val(ps.t, ps)
        _val = istrip ? t_str[4:prevind(t_str, sizeof(t_str), 3)] : t_str[2:prevind(t_str, sizeof(t_str))]
        expr = LITERAL(sfullspan, sspan,
            iscmd ? replace(_val, "\\`" => "`") :
                    replace(_val, "\\\"" => "\""), ps.t.kind)
        if istrip
            adjust_lcp(expr)
            ret = EXPR{StringH}(Any[expr], sfullspan, sspan)
        else
            return expr
        end
    else
        ret = EXPR{StringH}(Any[], sfullspan, sspan)
        input = IOBuffer(val(ps.t, ps))
        startbytes = istrip ? 3 : 1
        seek(input, startbytes)
        b = IOBuffer()
        while true
            if eof(input)
                lspan = position(b)
                # str = tostr(b)
                if b.size == 0
                # if sizeof(str) == 0
                    ex = ErrorToken()
                elseif istrip
                    str = tostr(b)
                    str = str[1:prevind(str, prevind(str, sizeof(str), 2))]
                    ex = LITERAL(lspan + ps.nt.startbyte - ps.t.endbyte - 1 + startbytes, lspan + startbytes, str, Tokens.STRING)
                else
                    str = tostr(b)
                    str =  str[1:prevind(str, sizeof(str))]
                    ex = LITERAL(lspan + ps.nt.startbyte - ps.t.endbyte - 1 + startbytes, lspan + startbytes, str, Tokens.STRING)
                end
                push!(ret.args, ex)
                istrip && adjust_lcp(ex, true)
                break
            end
            c = read(input, Char)
            if c == '\\'
                write(b, c)
                write(b, read(input, Char))
            elseif c == '$'
                lspan = position(b)
                str = tostr(b)
                ex = LITERAL(lspan + startbytes, lspan + startbytes, str, Tokens.STRING)
                push!(ret.args, ex); istrip && adjust_lcp(ex)
                startbytes = 0
                op = OPERATOR(1, 1, Tokens.EX_OR, false)
                if peekchar(input) == '('
                    lparen = PUNCTUATION(Tokens.LPAREN, 1, 1)
                    rparen = PUNCTUATION(Tokens.RPAREN, 1, 1)
                    skip(input, 1)
                    ps1 = ParseState(input)

                    if ps1.nt.kind == Tokens.RPAREN
                        call = UnarySyntaxOpCall(op, EXPR{InvisBrackets}(Any[lparen, rparen]))
                        push!(ret.args, call)
                        skip(input, 1)
                    else
                        interp = @closer ps1 paren parse_expression(ps1)
                        call = UnarySyntaxOpCall(op, EXPR{InvisBrackets}(Any[lparen, interp, rparen]))
                        push!(ret.args, call)
                        seek(input, ps1.nt.startbyte + 1)
                    end
                    # Compared to flisp/JuliaParser, we have an extra lookahead token,
                    # so we need to back up one here
                else
                    pos = position(input)
                    ps1 = ParseState(input)
                    next(ps1)
                    if ps1.t.kind == Tokens.WHITESPACE
                        t = EXPR{ErrorToken}([], ps.t.startbyte, ps.t.endbyte - ps.t.startbyte + 1)
                    else
                        t = INSTANCE(ps1)
                    end
                    # Attribute trailing whitespace to the string
                    t = adjustspan(t)
                    call = UnarySyntaxOpCall(op, t)
                    push!(ret.args, call)
                    seek(input, pos + t.fullspan)
                end
            else
                write(b, c)
            end
        end
    end

    single_string_T = (Tokens.STRING,ps.t.kind)
    if istrip
        if lcp != nothing && !isempty(lcp)
            for expr in exprs_to_adjust
                for (i, a) in enumerate(ret.args)
                    if expr == a
                        ret.args[i] = typeof(a)(expr.fullspan, expr.span, replace(expr.val, "\n$lcp" => "\n"), expr.kind)
                    end
                end
            end
        end
        # Drop leading newline
        if ret.args[1] isa LITERAL && ret.args[1].kind in single_string_T &&
                !isempty(ret.args[1].val) && ret.args[1].val[1] == '\n'
            ret.args[1] = dropleadlingnewline(ret.args[1])
        end
    end

    if (length(ret.args) == 1 && ret.args[1] isa LITERAL && ret.args[1].kind in single_string_T)
        ret = ret.args[1]::LITERAL
    end
    update_span!(ret)

    return ret
end

adjustspan(x::IDENTIFIER) = IDENTIFIER(x.span, x.span, x.val)
adjustspan(x::KEYWORD)= KEYWORD(x.kind, x.span, x.span)
adjustspan(x::OPERATOR) = OPERATOR(x.span, x.span, x.kind, x.dot)
adjustspan(x::LITERAL) = LITERAL(x.span, x.span, x.val, x.kind)
adjustspan(x::PUNCTUATION) = PUNCTUATION(x.kind, x.span, x.span)
function adjustspan(x::EXPR)
    x.fullspan = x.span
    return x
end

dropleadlingnewline(x::LITERAL) = LITERAL(x.fullspan, x.span, x.val[2:end], x.kind)


text:
using DataFrames
using DataFramesMeta
import CSV
using StatsPlots
using Statistics
using EcologicalNetworks
using EcologicalNetworksPlots
using TSne
using GLM

# Get the main predictions
predictions = DataFrame(CSV.File("hpc/outputs/predictions.csv"))

logistic = (x) -> 1.0 / (1.0 + exp(-x))
logit = (p) -> log(p/(1.0-p))

predictions.evidence = (predictions.updated ./ predictions.initial) .- 1.0
predictions.P = logistic.(predictions.evidence)

# Get clover
clover_df = joinpath("data", "clover.csv") |> CSV.File |> DataFrame
hosts = sort(unique(clover_df.Host))
viruses = sort(unique(clover_df.Virus))
A = zeros(Bool, (length(viruses), length(hosts)))
clover = BipartiteNetwork(A, viruses, hosts)
for clover_row in eachrow(clover_df)
    clover[clover_row.Virus, clover_row.Host] = true
end

# Correct the dataframe
for r in eachrow(predictions)
    r.value = clover[r.virus_pos, r.host_pos]
end

imputed = @linq predictions |>
    where(:value .== false) |>
    where(:P .>= 0.846847) |>
    select(:value, :host, :virus, :evidence, :P) |>
    orderby(:evidence)

CSV.write("artifacts/imputed_associations.csv", select(imputed, Not(:value)))

zoonoses = @where(imputed, :host .== "Homo sapiens")
CSV.write("artifacts/zoonoses.csv", select(zoonoses, Not(:value)))

imputed_clover = copy(clover)
for r in eachrow(imputed)
    imputed_clover[r.virus, r.host] = true
end

# Embeddings for Nardus
nohuman = filter(!isequal("Homo sapiens"), species(clover, dims=2))
rnk = 12
Lc, Rc = rdpg(clover, rnk)
Lt, Rt = rdpg(imputed_clover, rnk)
Lnc, Rnc = rdpg(clover[:,nohuman], rnk)
Lnt, Rnt = rdpg(imputed_clover[:,nohuman], rnk)
df = DataFrame(virus = String[], rank = Int[], clover = Float64[], trefle = Float64[], clover_nohuman = Float64[], trefle_nohuman = Float64[])
for i in 1:size(Lc, 1)
    for j in 1:size(Lc, 2)
        push!(df, (species(clover; dims=1)[i], j, Lc[i,j], Lt[i,j], Lnc[i,j], Lnt[i,j]))
    end
end
CSV.write(joinpath("artifacts", "viral_subspace.csv"), df)

# Writing the predictions
d = DataFrame(virus = String[], host = String[])
for int in interactions(imputed_clover)
    push!(d, (int.from, int.to))
end
sort!(d, :virus)
CSV.write("artifacts/trefle.csv", d)

# Degree distribution
function Pk(N::T; dims::Union{Nothing,Integer}=nothing) where {T<:AbstractEcologicalNetwork}
    deg = collect(values(degree(N; dims=dims)))
    u = sort(unique(deg))
    p = zeros(Float64, length(u))
    for (i,x) in enumerate(u)
        p[i] = sum(deg.==x)
    end
    return (u, p./sum(p))
end

#=
plot(Pk(clover; dims=1)..., lab="Raw data", c=:grey, ls=:dash)
plot!(Pk(imputed_clover; dims=1)..., lab="Post imputation", c=:black)
xaxis!(:log, "Virus degree", (1, 1e3))
yaxis!(:log, (0.0005, 1.0), "Frequency")
savefig("mainfigs/degree-virus.png")

plot(Pk(clover; dims=2)..., lab="Raw data", c=:grey, ls=:dash)
plot!(Pk(imputed_clover; dims=2)..., lab="Post imputation", c=:black)
xaxis!(:log, "Host degree", (1, 1e3))
yaxis!(:log, (0.0005, 1.0), "Frequency")
savefig("mainfigs/degree-host.png")

plot(Pk(clover)..., lab="Raw data", c=:grey, ls=:dash)
plot!(Pk(imputed_clover)..., lab="Post imputation", c=:black)
xaxis!(:log, "Degree", (1, 1e3))
yaxis!(:log, (0.0005, 1.0), "Frequency")
savefig("mainfigs/degree-global.png")
=#

# Embedding
UCLOV = EcologicalNetworks.mirror(convert(UnipartiteNetwork, clover))
emb_clover = tsne(convert.(Float64, Array(UCLOV.edges)), 2, 0, 2000, 6)

IO = initial(RandomInitialLayout, UCLOV)
for (i,s) in enumerate(species(UCLOV))
    IO[s].x = emb_clover[i,1]
    IO[s].y = emb_clover[i,2]
end

UIMPT = EcologicalNetworks.mirror(convert(UnipartiteNetwork, imputed_clover))
emb_impt = tsne(convert.(Float64, Array(UIMPT.edges)), 2, 0, 2000, 6)

IM = initial(RandomInitialLayout, UIMPT)
for (i,s) in enumerate(species(UIMPT))
    IM[s].x = emb_impt[i,1]
    IM[s].y = emb_impt[i,2]
end

p_orig = scatter(IO, clover, bipartite=true, nodesize=degree(clover), msc=:skyblue, msw=0.5, aspectratio=1, dpi=300)
p_imput = scatter(IM, imputed_clover, bipartite=true, nodesize=degree(imputed_clover), msc=:teal, msw=0.5, aspectratio=1, dpi=300)

plot(p_orig, p_imput)
savefig("figures/before-after.png")

# Number of paths
top_zoo = zoonoses.virus#vec(zoonoses[end-19:end,:].virus)
TMP = UnipartiteQuantitativeNetwork(number_of_paths(UCLOV; n=3), EcologicalNetworks._species_objects(UCLOV)...)
i = interactions(TMP)
filter!(int -> int.from != int.to, i)
filter!(int -> "Homo sapiens" in [int.from], i)
filter!(int -> int.to in top_zoo, i)
df = DataFrame(virus = String[], paths = Int64[])
for int in i
    push!(df, (int.to, int.strength))
end

path_zoo = leftjoin(zoonoses, df, on=:virus)

# Residuals from number of paths
X = log.(path_zoo.evidence)
y = log.(path_zoo.paths)

regmodel = lm(permutedims(permutedims(X)), y)

path_zoo.residuals = residuals(regmodel)
sort!(path_zoo, :residuals, rev=true)
select!(path_zoo, Not(:value))
select!(path_zoo, Not(:host))

@df path_zoo scatter(:evidence, :paths, frame=:box, lab="", dpi=400, marker_z=path_zoo.residuals, aspectratio=1, clim=(-3,3), c=:PRGn)
xaxis!(:log, (0.9, 500), "Evidence of interaction")
yaxis!(:log, (0.9, 500), "Number of paths to human")
plot!(x -> coef(regmodel)[1]x, lab="", lw=2.0, c=:darkgrey, ls=:dash)
savefig("figures/number_of_paths.png")

tax = unique(select(clover_df, [:VirusClass, :VirusOrder, :VirusFamily, :VirusGenus, :Virus]))
path_zoo = leftjoin(path_zoo, tax, on = :virus => :Virus, matchmissing = :equal)

@df dropmissing(select(path_zoo, [:residuals, :VirusOrder])) dotplot(:VirusOrder, :residuals, group=:VirusOrder, legend=:outerleft, size=(1000, 400), frame=:zerolines)
xaxis!(rotation=90)

CSV.write(joinpath("artifacts", "number_of_paths.csv"), path_zoo)

# Overlap analysis
overlap_results = DataFrame(sp = String[], step = Symbol[], score = Float64[])

raw_ajs = AJS(clover; dims=2)
raw_human = filter(p -> "Homo sapiens" in p.first, raw_ajs)
for r in raw_human
    sp = filter(p -> p != "Homo sapiens", collect(r.first))[1]
    push!(overlap_results, (sp, :initial, r.second))
end

imp_ajs = AJS(imputed_clover; dims=2)
imp_human = filter(p -> "Homo sapiens" in p.first, imp_ajs)
for r in imp_human
    sp = filter(p -> p != "Homo sapiens", collect(r.first))[1]
    push!(overlap_results, (sp, :imputed, r.second))
end

top_10_initial = @linq overlap_results |>
    where(:step .== ^(:initial)) |>
    orderby(:score)
top_10_initial = last(top_10_initial, 10)

top_10_imputed = @linq overlap_results |>
    where(:step .== ^(:imputed)) |>
    orderby(:score)
top_10_imputed = last(top_10_imputed, 10)

plot(grid=false, xticks=((1,2),["Initial", "Imputed"]), legend=false, dpi=300)
raw_order = LinRange(minimum(top_10_initial.score)-0.02, maximum(top_10_initial.score)+0.02, 10)
for (i,r) in enumerate(raw_order)
    plot!([0.88, 1.0], [raw_order[i], top_10_initial.score[i]], lab="", c=:darkgrey)
end
annotate!([(0.85, raw_order[i], text(r.sp, :right, 6)) for (i,r) in enumerate(eachrow(top_10_initial))])
imp_order = LinRange(minimum(top_10_imputed.score)-0.02, maximum(top_10_imputed.score)+0.02, 10)
for (i,r) in enumerate(imp_order)
    plot!([2.0, 2.12], [top_10_imputed.score[i], r], lab="", c=:darkgrey)
end
annotate!([(2.15, imp_order[i], text(r.sp, :left, 6)) for (i,r) in enumerate(eachrow(top_10_imputed))])
for rinit in eachrow(top_10_initial)
    for rimpt in eachrow(top_10_imputed)
        if rinit.sp == rimpt.sp
            plot!([1, 2], [rinit.score, rimpt.score], c=:grey, ls=:dot)
        end
    end
end
@df top_10_initial scatter!(fill(1, 10), :score, c=:black, msw=0.0)
@df top_10_imputed scatter!(fill(2, 10), :score, c=:black, msw=0.0)
xaxis!((-0.1,3.1), false)
yaxis!((0.0,0.75), "Similarity")
savefig("figures/human-similarity.png")

# Co-occurrences in datasets
imputed.cooc = fill(false, size(imputed, 1))
for i in 1:size(imputed, 1)
    v = imputed.virus[i]
    h = imputed.host[i]
    dbhost = unique(@where(clover_df, :Host .== h).Database)
    dbvirus = unique(@where(clover_df, :Virus .== v).Database)
    imputed.cooc[i] = length(dbvirus ∩ dbhost) > 0
end

@df imputed density(:P, group=:cooc, frame=:box, legend=:topleft, dpi=300)
xaxis!((0.846847, 1.0), "Probability")
yaxis!((0, 45),)
savefig("figures/probability-by-cooccurrence.png")

# Phylogeny - rank correlation
phylodist = DataFrame(CSV.File("artifacts/phylo_distance_to_human.csv"))
rename!(phylodist, "Column1" => "sp")
phylodist.sp = map(n -> replace(n, "_" => " "), phylodist.sp)

sharing_results = DataFrame(sp = String[], step = Symbol[], count = Int64[])

raw_share = overlap(clover; dims=2)
for r in filter(p -> "Homo sapiens" in p.first, raw_share)
    sp = filter(p -> p != "Homo sapiens", collect(r.first))[1]
    push!(sharing_results, (sp, :initial, r.second))
end

imp_share = overlap(imputed_clover; dims=2)
for r in filter(p -> "Homo sapiens" in p.first, imp_share)
    sp = filter(p -> p != "Homo sapiens", collect(r.first))[1]
    push!(sharing_results, (sp, :imputed, r.second))
end

sharing_results = leftjoin(sharing_results, select(overlap_results, :sp, :score), on=:sp)
phyloverlap = leftjoin(sharing_results, phylodist, on=:sp)
CSV.write("artifacts/sharing-phylogeny.csv", phyloverlap)

ext: jl

max_stars_repo_name: viralemergence/trefle

max_issues_repo_name: viralemergence/trefle

max_forks_repo_name: viralemergence/trefle


text:
using EDA
launch()


text:
# Gradient algorithm
function est_gradient!(dθ, θ, g, t; regressor=[0.0], output=0.0)
    w, y = regressor, output
    e = w'*θ - output
    dθ .= -g*w*e
    return θ
end

gradient_estimator(g) = (dθ, θ, p, t; kwargs...) -> est_gradient!(dθ, θ, g, t; kwargs...)


text:
# Most algorithms are written keeping the concepts and pseudo code of
# CLRS 3ed as close to as in the book. There has been minor deviations
# taken at places for programming and architectural reasons.

# These data structures assumes strong trichotomy of keys exist
# k1 and k2 if are valid keys then one of the following relationship
# has to be true.
# 1. k1  < k2
# 2. k1 == k2
# 3. k1  > k2
#
# However, to the user of these data structures have to implement
# `Base.isless` only.
# No other operators are used in the code so that there are no implicit
# overheads due to fallback options.

abstract type AbstractNode{K, V} end
abstract type AbstractBST{K, V} end

Base.isless(n1::T, n2::T) where {T <:AbstractNode} = Base.isless(n1.k, n2.k)
Base.isless(n::AbstractNode{K, V}, k::K) where {K, V} = isless(n.k, k)
Base.isless(k::K, n::AbstractNode{K, V}) where {K, V} = isless(k, n.k)

_k(n::AbstractNode) = n.k
_v(n::AbstractNode) = n.v
_l(n::AbstractNode) = n.l
_r(n::AbstractNode) = n.r
_p(n::AbstractNode) = n.p
_l!(t::T, x::N, y::N) where {K, V,
                             N <: AbstractNode{K, V},
                             T <: AbstractBST{K, V}} =
                                 ((x.l, y.p) = (y, (isnil(t, y) ? y.p : x)))
_r!(t::T, x::N, y::N) where {K, V,
                             N <: AbstractNode{K, V},
                             T <: AbstractBST{K, V}} =
                                 ((x.r, y.p) = (y, (isnil(t, y) ? y.p : x)))
_p!(x::T, y::T) where {T <: AbstractNode} = (x.p = y)

@inline function _extremum(dir::Function, n::AbstractNode, t::AbstractBST)
    while true
        nn = dir(n)
        isnil(t, nn) && return n
        n = nn
    end
end

_maximum(n::AbstractNode, t::AbstractBST) = _extremum(_r, n, t)
_minimum(n::AbstractNode, t::AbstractBST) = _extremum(_l, n, t)

_successor(x::AbstractNode, t::AbstractBST)   = _pred_succ(_r, _minimum, x, t)
_predecessor(x::AbstractNode, t::AbstractBST) = _pred_succ(_l, _maximum, x, t)

@inline function _pred_succ(f::Function, g::Function,
                            x::AbstractNode, t::AbstractBST)
    !isnil(t, f(x)) && return g(f(x), t)
    y = _p(x)
    while !isnil(t, y) && x === f(y)
        x = y
        y = _p(y)
    end
    return y
end

function _inorder(f::Function, n::AbstractNode, t::AbstractBST)
    if !isnil(t, n)
        proceed = true
        proceed = (proceed && _inorder(f, n.l, t))
        res = f(n)
        if res isa Bool
            proceed = (proceed && res)
        end
        proceed = (proceed && _inorder(f, n.r, t))
        return proceed
    else
        return true
    end
end

function node_print(t::AbstractBST{K, V},
                    n::AbstractNode{K, V},
                    prefix::String,
                    left::Bool) where {K, V}
    isnil(t, n) && return
    prefix *= prefix
    node_print(t, n.l, prefix, true)
    RED="\033[0;31m"
    NC="\033[0m"
    if n.red
        println(prefix, RED, n.k, NC)
    else
        println(prefix, n.k)
    end
    node_print(t, n.r, prefix, false)
end

@inline function _search(t::AbstractBST, n::AbstractNode{K, V}, k::K) where {K, V}
    while true
        if k < _k(n)
            nn = _l(n)
            d  = -1
        elseif _k(n) < k
            nn = _r(n)
            d = 1
        else
            return (n, 0)
        end
        isnil(t, nn) && return n, d
        n = nn
    end
end

Base.length(t::AbstractBST) = t.n
Base.isempty(t::AbstractBST) = Base.length(t) == 0

function Base.maximum(t::AbstractBST)
    isempty(t) && error("Empty tree cannot have a maximum")
    n = _maximum(t.root, t)
    return n.k => n.v
end

function Base.minimum(t::AbstractBST)
    isempty(t) && error("Empty tree cannot have a minimum")
    n = _minimum(t.root, t)
    return n.k => n.v
end

@inline function left_rotate!(t::T,
                              x::N) where {T <: AbstractBST, N <: AbstractNode}
    y = x.r
    x.r = y.l
    !isnil(t, y.l) && (y.l.p = x)
    y.p = x.p
    if isnil(t, x.p)
        t.root = y
    elseif x === x.p.l
        x.p.l = y
    else
        x.p.r = y
    end
    y.l = x
    x.p = y
    return
end

@inline function right_rotate!(t::T,
                               y::N) where {T <: AbstractBST, N <: AbstractNode}
    x = y.l
    y.l = x.r
    !isnil(t, x.r) && (x.r.p = y)
    x.p = y.p
    if isnil(t, y.p)
        t.root = x
    elseif y === y.p.r
        y.p.r = x
    else
        y.p.l = x
    end
    x.r = y
    y.p = x
    return
end

@inline function Base.delete!(t::AbstractBST{K, V}, k::K) where {K, V}
    isempty(t) && error("Cannot delete from empty tree")
    n, d = _search(t, t.root, k)
    if d == 0
        _delete!(t, n)
        return n.k => n.v
    else
        return nothing
    end
end

mutable struct BSTNode{K, V} <: AbstractNode{K, V}
    k::K
    v::V
    l::BSTNode{K, V}
    r::BSTNode{K, V}
    p::BSTNode{K, V}
    function BSTNode{K, V}() where {K, V}
        self = new{K, V}()
        self.l = self.r = self.p = self
    end
end

mutable struct BinarySearchTree{K, V} <: AbstractBST{K, V}
    root::BSTNode{K, V}
    nil::BSTNode{K, V}
    n::Int
    unique::Bool
    function BinarySearchTree{K, V}() where {K, V}
        s = new{K, V}()
        nil = BSTNode{K, V}()
        @assert nil.l === nil && nil.r === nil && nil.p === nil
        s.n = 0
        s.unique = false
        s.root = s.nil = nil
        return s
    end
end

@inline function Base.get(t::AbstractBST{K, V}, k::K, v::V) where {K, V}
    isempty(t) && return v
    n, d = _search(t, t.root, k)
    d != 0 && return v
    return n.v
end

@inline function Base.get!(t::AbstractBST{K, V}, k::K, v::V) where {K, V}
    if isempty(t)
        insert!(t, k, v)
        return v
    end
    n, d = _search(t, t.root, k)
    if d != 0
        insert!(t, k, v)
        return v
    end
    return n.v
end

function BSTNode(t::BinarySearchTree{K, V}, k::K, v::V) where {K, V}
    s = BSTNode{K, V}()
    s.k, s.v, s.l, s.r, s.p = k, v, t.nil, t.nil, t.nil
    return s
end

isnil(t::BinarySearchTree, n::BSTNode) = n === t.nil
Base.empty!(t::BinarySearchTree) = ((t.root, t.n) = (t.nil, 0))

function Base.insert!(t::BinarySearchTree, k::K, v::V) where {K, V}
    t.root = t.n == 0 ? BSTNode(t, k, v) : _insert!(t, t.root, k, v, t.unique)
    t.n += 1
    return
end

@inline function _insert!(t::BinarySearchTree,
                          n::BSTNode{K, V},
                          k::K, v::V,
                          unique::Bool=true) where {K, V}
    tn = ni = n
    while !isnil(t, tn)
        ni = tn
        tn = k < ni.k ? ni.l : (!unique || (ni.k < k)) ? ni.r :
            unique && error("Key $k already exists.")
    end
    nn = BSTNode(t, k, v)
    if k < ni.k
        _l!(t, ni, nn)
    else
        _r!(t, ni, nn)
    end
    return n
end

@inline function  _transplant!(t::BinarySearchTree, u::BSTNode, v::BSTNode)
    if isnil(t, u.p)
        t.root = v
    elseif u === u.p.l
        u.p.l = v
    else
        u.p.r = v
    end
    v.p = u.p
end

function _delete!(t::BinarySearchTree, z::BSTNode)
    if isnil(t, z.l)
        _transplant!(t, z, z.r)
    elseif isnil(t, z.r)
        _transplant!(t, z, z.l)
    else
        y = _minimum(z.r, t)
        if y.p !== z
            _transplant!(t, y, y.r)
            y.r = z.r
            y.r.p = y
        end
        _transplant!(t, z, y)
        y.l = z.l
        y.l.p = y
    end
    t.n -= 1
    return z
end

#=
CLRS 3rd Ed. Chapter 13

1. Every node is either red or black.
2. The root is black.
3. Every leaf (NIL ) is black.
4. If a node is red, then both its children are black.
5. For each node, all simple paths from the node to descendant leaves contain the
same number of black nodes.

=#

mutable struct RBNode{K, V} <: AbstractNode{K, V}
    k::K
    v::V
    red::Bool
    l::RBNode{K, V}
    r::RBNode{K, V}
    p::RBNode{K, V}
    function RBNode{K, V}() where {K, V}
        self = new{K, V}()
        self.red = false
        self.l = self.r = self.p = self
    end
end

mutable struct RBTree{K, V} <: AbstractBST{K, V}
    root::RBNode{K, V}
    nil::RBNode{K, V}
    n::Int
    unique::Bool
    function RBTree{K, V}() where {K, V}
        s = new{K, V}()
        nil = RBNode{K, V}()
        @assert nil.l === nil && nil.r === nil && nil.p === nil
        s.n = 0
        s.unique = false
        s.root = s.nil = nil
        return s
    end
end

function RBNode(t::RBTree{K, V}, k::K, v::V) where {K, V}
    s = RBNode{K, V}()
    s.k, s.v, s.red, s.l, s.r, s.p = k, v, false, t.nil, t.nil, t.nil
    return s
end

function Base.show(io::IO, t::AbstractBST)
    println(io, "$(typeof(t)) Tree with $(length(t)) nodes.")
    !isnil(t, t.root) && println(io, "Root at: $(t.root.k).")
end

isnil(t::RBTree, n::RBNode) = n === t.nil
Base.empty!(t::RBTree) = (t.root = t.nil; t.n = 0; nothing)

function Base.insert!(t::RBTree{K, V}, k::K, v::V) where {K, V}
    z = RBNode(t, k, v)
    y::RBNode{K, V} = t.nil
    x::RBNode{K, V} = t.root

    while x !== t.nil
        y = x
        x = k < x.k ? x.l :
            !t.unique || (x.k < k) ? x.r :
            t.unique && error("Key $k already exists.")
    end
    z.p = y
    if y === t.nil
        t.root = z
    elseif k < y.k
        y.l = z
    else
        y.r = z
    end
    z.red = true
    _insert_fixup!(t, z)
    t.n += 1
    return t
end

@inline function _insert_fixup!(t::RBTree, z::RBNode)
    while z.p.red
        if z.p === z.p.p.l
            y = z.p.p.r
            if y.red
                z.p.red = false
                y.red = false
                z.p.p.red = true
                z = z.p.p
            else
                if z === z.p.r
                    z = z.p
                    left_rotate!(t, z)
                end
                z.p.red = false
                z.p.p.red = true
                right_rotate!(t, z.p.p)
            end
        else
            y = z.p.p.l
            if y.red
                z.p.red = false
                y.red = false
                z.p.p.red = true
                z = z.p.p
            else
                if z === z.p.l
                    z = z.p
                    right_rotate!(t, z)
                end
                z.p.red = false
                z.p.p.red = true
                left_rotate!(t, z.p.p)
            end
        end
    end
    t.root.red = false
end

# Intermediate step used in delete!

@inline function  _transplant!(t::RBTree, u::RBNode, v::RBNode)
    if u.p === t.nil
        t.root = v
    elseif u === u.p.l
        u.p.l = v
    else
        u.p.r = v
    end
    v.p = u.p
end

@inline function _delete!(t::RBTree, z::RBNode)
    y = z
    y_is_red = y.red
    if z.l === t.nil
        x = z.r
        _transplant!(t, z, z.r)
    elseif z.r === t.nil
        x = z.l
        _transplant!(t, z, z.l)
    else
        y = _minimum(z.r, t)
        y_is_red = y.red
        x = y.r
        if y.p === z
            x.p = y
        else
            _transplant!(t, y, y.r)
            y.r = z.r
            y.r.p = y
        end
        _transplant!(t, z, y)
        y.l = z.l
        y.l.p = y
        y.red = z.red
    end
    if !y_is_red
        _delete_fixup!(t, x)
    end
    t.n -= 1
    return z
end

function _delete_fixup!(t::RBTree, x::RBNode)
    while x !== t.root && !x.red
        if x === x.p.l
            w = x.p.r
            if w.red
                w.red = false
                x.p.red = true
                left_rotate!(t, x.p)
                w = x.p.r
            end
            if !w.l.red && !w.r.red
                w.red = true
                x = x.p
            else
                if !w.r.red
                    w.l.red = false
                    w.red = true
                    right_rotate!(t, w)
                    w = x.p.r
                end
                w.red = x.p.red
                x.p.red = false
                w.r.red = false
                left_rotate!(t, x.p)
                x = t.root
            end
        else
            w = x.p.l
            if w.red
                w.red = false
                x.p.red = true
                right_rotate!(t, x.p)
                w = x.p.l
            end
            if !w.l.red && !w.r.red
                w.red = true
                x = x.p
            else
                if !w.l.red
                    w.r.red = false
                    w.red = true
                    left_rotate!(t, w)
                    w = x.p.l
                end
                w.red = x.p.red
                x.p.red = false
                w.l.red = false
                right_rotate!(t, x.p)
                x = t.root
            end
        end
    end
    x.red = false
end

mutable struct Iterator{K, V, T, N}
    tree::T
    from::N
    to::N
    function Iterator{K, V, T, N}(t::T,
                                  from::N,
                                  to::N) where {K, V,
                                                T <: AbstractBST{K, V},
                                                N <: AbstractNode{K, V}}
        @assert isnil(t, from) || isnil(t, to) || !(to.k < from.k)
        "`from` value cannot be more that the `to` value"
        new{K, V, T, N}(t, from, to)
    end
end

function Iterator(t::T, from::N=_minimum(t.root, t),
                  to::N=t.nil) where {K, V,
                                      T <: AbstractBST{K, V},
                                      N <: AbstractNode{K, V}}
    Iterator{K, V, T, N}(t, from, to)
end

function Iterator(t::T, from::K, to::K) where {K, V, T <: AbstractBST{K, V}}
    to < from && error("Cannot initialize iterator where `from > to`.")
    n, d = _search(t, t.root, from)
    if d == 0
        nn = _predecessor(n, t)
        while nn.k == n.k && nn !== n
            n = nn
            nn = _predecessor(n, t)
        end
    elseif d > 0
        n = _successor(n, t)
        isnil(t, n) && return Iterator(t, t.nil, t.nil)
    end
    fromN = n
    to < fromN && return Iterator(t, t.nil, t.nil)
    n, d = _search(t, t.root, to)
    nn = n
    if d == 0
        nn = _successor(n, t)
        while nn.k == n.k && !isnil(t, nn)
            n = nn
            nn = _successor(n, t)
        end
        n = nn
    elseif d < 0
        n = _predecessor(n, t)
    else
        n = _successor(n, t)
    end
    toN = n
    return Iterator(t, fromN, toN)
end

Base.IteratorSize(it::Iterator) = Base.SizeUnknown()

Base.iterate(it::Iterator) = iterate(it, it.from)

function Base.iterate(it::Iterator{K, V, T, N},
                      n::N) where {K, V,
                                   T <: AbstractBST{K, V},
                                   N <: AbstractNode{K, V}}
    n === it.to && return nothing
    return ((n.k => n.v), _successor(n, it.tree))
end

#
# Based on -
# Random binary search tree with equal elements
#   by Tomi A. Pasanen
#
# Published in Theoretical Computer Science 411 (2010) 3867–3872

mutable struct RandomizedBSTNode{K, V} <: AbstractNode{K, V}
    k::K
    v::V
    l::RandomizedBSTNode{K, V}
    r::RandomizedBSTNode{K, V}
    p::RandomizedBSTNode{K, V}
    n::Int

    function RandomizedBSTNode{K, V}() where {K, V}
        self = new{K, V}()
        self.n = 0
        self.l = self.r = self.p = self
    end
end

mutable struct RandomizedBinarySearchTree{K, V}  <: AbstractBST{K, V}
    root::RandomizedBSTNode{K, V}
    nil::RandomizedBSTNode{K, V}
    unique::Bool

    function RandomizedBinarySearchTree{K, V}() where {K, V}
        self = new{K, V}()
        nil = RandomizedBSTNode{K, V}()
        @assert nil.l === nil && nil.r === nil && nil.n == 0
        self.unique = false
        self.root = self.nil = nil
        return self
    end
end

function RandomizedBSTNode(t::RandomizedBinarySearchTree{K, V}, k::K, v::V) where {K, V}
    self = RandomizedBSTNode{K, V}()
    self.k, self.v, self.l, self.r, self.p, self.n = k, v, t.nil, t.nil, t.nil, 0
    return self
end

function _inorder2(f::Function, n::AbstractNode, t::AbstractBST)
    if !isnil(t, n)
        proceed = true
        proceed = (proceed && _inorder2(f, n.l, t))
        res = f(n, t)
        if res isa Bool
            proceed = (proceed && res)
        end
        proceed = (proceed && _inorder2(f, n.r, t))
        return proceed
    else
        return true
    end
end

function print_tree(n::RandomizedBinarySearchTree{K, V}) where {K, V}
    if isempty(n)
        println("[]")
    else
        println("Root = ", n.root.k, " and size = ", n.root.n)
        _inorder2(print_node, n.root, n)
    end
end

function print_node(n::RandomizedBSTNode{K, V}, t::RandomizedBinarySearchTree{K, V}) where {K, V}
  str = "[$(n.k)], left = $(n.l == t.nil ? "NULL" : n.l.k), right = $(n.r == t.nil ? "NULL" : n.r.k), parent = $(n.p == t.nil ? "NIL" : n.p.k) and size=$(n.n)"
  println(str)
  return str
end

isnil(t::RandomizedBinarySearchTree, n::RandomizedBSTNode) = n === t.nil
Base.empty!(t::RandomizedBinarySearchTree) = (t.root = t.nil)
Base.length(t::RandomizedBinarySearchTree) = t.root.n

function Base.insert!(t::RandomizedBinarySearchTree, k::K, v::V) where {K, V}
    r = _insert!(t.root, k, v, t, t.unique)
end

@inline function _insert!(node::RandomizedBSTNode{K, V},
                          k::K, v::V,
                          t::RandomizedBinarySearchTree, unique::Bool=true) where {K, V}
    rn = rand(0:node.n)
    if rn == node.n
        r = _insert_at_root!(node, k, v, unique, t)
        _fix_node!(t, r)
        return r
    else
        if k < node.k
            node.l = _insert!(node.l, k, v, t, unique)
            _fix_node!(t, node.l)
            _fix_node!(t, node)
        elseif !unique || (node.k < k)
            node.r =  _insert!(node.r, k, v, t, unique)
            _fix_node!(t, node.r)
            _fix_node!(t, node)
        else
            unique && error("Key $k already exists.")
        end
    end
    return node
end

@inline function _fix_node!(t::RandomizedBinarySearchTree{K, V}, n::RandomizedBSTNode{K, V}) where {K, V}
    if n != t.nil
        n.n = n.l.n + n.r.n + 1
        n.l != t.nil ? n.l.p = n : nothing
        n.r != t.nil ? n.r.p = n : nothing
    end
end

@inline function _insert_at_root!(t::RandomizedBSTNode{K, V},
                                 k::K, v::V, unique::Bool,
                                 tree::RandomizedBinarySearchTree{K, V}) where {K, V}
    node = RandomizedBSTNode(tree, k, v)
    node.p = t.p
    l, r = _split(t, k, v, unique, tree)
    node.l, node.r = l, r
    _fix_node!(tree, node)
    if isempty(tree) || tree.root == t
        tree.root = node
    end
    return node
end

# Given a k and a subtree - node returns (l, r) where all elements in l are < than k and all elements in r are > k
@inline function _split(node::RandomizedBSTNode{K, V}, k::K, v::V, unique::Bool,
                        t::RandomizedBinarySearchTree) where {K, V}
    if node == t.nil
       return t.nil, t.nil
    end

    if k < node.k
        r = node
        l, r.l = _split(node.l, k, v, unique, t)
        _fix_node!(t, r.l)
        _fix_node!(t, r)
        _fix_node!(t, l)
        return l, r
    elseif !unique || node.k < k
        l = node
        l.r, r = _split(node.r, k, v, unique, t)
        _fix_node!(t, l.r)
        _fix_node!(t, l)
        _fix_node!(t, r)
        return l, r
    else
        unique && error("Key $k already exists.")
    end
end

@inline function Base.delete!(t::RandomizedBinarySearchTree{K, V}, k::K) where {K, V}
    isempty(t) && error("Cannot delete from empty tree")
    x, s = _delete!(t.root, k, t)
    return s != t.nil ? s.k => s.v : nothing
end

@inline function _delete!(node::RandomizedBSTNode{K, V}, k::K, t::RandomizedBinarySearchTree{K, V}) where {K, V}
    found = t.nil
    if node == t.nil
        return t.nil, t.nil
    end
    if k < node.k
        node.l, found = _delete!(node.l, k, t)
    elseif k > node.k
        node.r, found = _delete!(node.r, k, t)
    else
        found = node
        node = _join(node.l,node.r, t)
        if (found == t.root)
            t.root = node
            node.p = t.nil
        end
    end
    _fix_node!(t, node)
    return node, found
end

@inline function _join(l::RandomizedBSTNode{K, V}, r::RandomizedBSTNode{K, V}, t::RandomizedBinarySearchTree{K, V}) where {K, V}
    m, n, total = l.n, r.n, l.n + r.n

    if total == 0
        return t.nil
    end

    rn = rand(0:total-1)
    if (rn < m)
        l.r = _join(l.r, r, t)
        _fix_node!(t, l.r)
        _fix_node!(t, l)
        return l
    else
         r.l = _join(l, r.l, t)
         _fix_node!(t, r.l)
         _fix_node!(t, r)
         return r
    end
end


text:
###############################################################################
#
# Leaf parameters container
#
###############################################################################
"""
    mutable struct Leaf{FT}

Struct to store leaf information.

# Fields
$(TYPEDFIELDS)
"""
Base.@kwdef mutable struct Leaf{FT<:AbstractFloat}
    # Temperature related
    "Temperature `[K]`"
    T::FT = 298.15
    "Old Temperature `[K]`, if not T, run leaf_temperature_dependence!"
    T_old::FT = 0

    # Photosynthesis system
    "Rate constant for thermal dissipation"
    Kd::FT = 0.85
    "Rate constant for fluorescence (const)"
    Kf::FT = 0.05
    "Reversible NPQ rate constant (initially zero)"
    Kr::FT = 0
    "Sustained NPQ rate constant (for seasonal changes, default is zero)"
    Ks::FT = 0
    "Rate constant for photochemistry (all reaction centers open)"
    Kp::FT = 4
    "Maximal Kp"
    Kp_max::FT = 4
    "max PSII yield (Kr=0, all RC open)"
    maxPSII::FT = Kp/(Kp+Kf+Kd)
    "Fraction of absorbed light used by PSII ETR"
    PSII_frac::FT = 0.5

    # CO₂ pressures
    "Leaf internal CO₂ partial pressure `[Pa]`"
    p_i::FT = 10
    "Leaf surface CO₂ partial pressure `[Pa]`"
    p_s::FT = 40
    "Saturation H₂O vapor pressure `[Pa]`"
    p_sat::FT = saturation_vapor_pressure(T)
    "Leaf diffusive conductance to CO₂ `[mol m⁻² s⁻¹]`"
    g_bc::FT = 3 / 1.35
    "Leaf diffusive conductance to CO₂ `[mol m⁻² s⁻¹]`"
    g_lc::FT = 0.01

    # Photosynthesis related
    "RubisCO limited photosynthetic rate `[μmol m⁻² s⁻¹]`"
    Ac::FT = 0
    "Light limited photosynthetic rate `[μmol m⁻² s⁻¹]`"
    Aj::FT = 0
    "Gross photosynthetic rate `[μmol m⁻² s⁻¹]`"
    Ag::FT = 0
    "Net photosynthetic rate `[μmol m⁻² s⁻¹]`"
    An::FT = 0
    "Product limited photosynthetic rate `[μmol m⁻² s⁻¹]`"
    Ap::FT = 0
    "Electron transport `[μmol m⁻² s⁻¹]`"
    J::FT = 0
    "Potential Electron Transport Rate `[μmol m⁻² s⁻¹]`"
    J_pot::FT = 0
    "Maximal electron transport rate `[μmol m⁻² s⁻¹]`"
    Jmax::FT = 120
    "Maximal electron transport rate at 298.15 K `[μmol m⁻² s⁻¹]`"
    Jmax25::FT = 120
    "RubisCO coefficient Kc `[Pa]`"
    Kc::FT = 0
    "RubisCO coefficient Ko `[Pa]`"
    Ko::FT = 0
    "PEP coefficient Ko `[Pa]`"
    Kpep::FT = 0
    "Michaelis-Menten's coefficient `[Pa]`"
    Km::FT = 0
    "Respiration rate `[μmol m⁻² s⁻¹]`"
    Rd::FT = 1
    "Respiration rate at 298.15 K `[μmol m⁻² s⁻¹]`"
    Rd25::FT = 1
    "Maximal carboxylation rate `[μmol m⁻² s⁻¹]`"
    Vcmax::FT = 60
    "Maximal carboxylation rate at 298.15 K `[μmol m⁻² s⁻¹]`"
    Vcmax25::FT = 60
    "Maximal PEP carboxylation rate `[μmol m⁻² s⁻¹]`"
    Vpmax::FT = 120
    "Maximal PEP carboxylation rate at 298.15 K `[μmol m⁻² s⁻¹]`"
    Vpmax25::FT = 120
    "CO₂ compensation point with the absence of Rd `[Pa]`"
    Γ_star::FT = 0

    # Cytochrome related
    "Total concentration of Cytochrome b₆f `[μmol m⁻²]`"
    C_b₆f::FT = 350 / 300
    "Maximal turnover rate of Cytochrome b₆f `[e⁻ s⁻¹]`"
    k_q::FT = 300
    "Maximal Cytochrome b₆f activity `[μmol e⁻ m⁻² s⁻¹]`"
    Vqmax::FT = C_b₆f * k_q
    "Rate constant of consititutive heat loss from the antennae `[s⁻¹]`"
    K_D1::FT = 5.5e8
    "rate constant of fluorescence `[s⁻¹]`"
    K_F1::FT = 5e7
    "Rate constant of regulated heat loss for PS I `[s⁻¹]`"
    K_N1::FT = 1.45e10
    "Rate constant of photochemistry for PS I `[s⁻¹]`"
    K_P1::FT = 1.45e10
    "Rate constant of photochemistry for PS II `[s⁻¹]`"
    K_P2::FT = 4.5e9
    "Rate constant of excitation sharing for PS II `[s⁻¹]`"
    K_U2::FT = 0
    "PPFD absorbed by PS I per incident PPFD"
    α_1::FT = 0.41
    "PPFD absorbed by PS II per incident PPFD"
    α_2::FT = 0.44
    "Weighting factor for PS I"
    ϵ_1::FT = 0
    "Weighting factor for PS II"
    ϵ_2::FT = 1
    "Maximal PS I photochemical yield"
    φ_P1_max::FT = K_P1 / (K_P1 + K_D1 + K_F1)
    "Coupling efficiency of cyclic electron flow `[mol ATP mol⁻¹ e⁻]`"
    n_C::FT = 1
    "Coupling efficiency of linear electron flow `[mol ATP mol⁻¹ e⁻]`"
    n_L::FT = 0.75
    "ratio between J_P700 and J_P680"
    η::FT = 0

    # related to fluorescence
    "PS II electron transport rate `[μmol e⁻ m⁻² s⁻¹]`"
    J_P680_a::FT = 0
    "Rubisco limited PS II electron transport rate `[μmol e⁻ m⁻² s⁻¹]`"
    J_P680_c::FT = 0
    "Light limited PS II electron transport rate `[μmol e⁻ m⁻² s⁻¹]`"
    J_P680_j::FT = 0
    "Product limited PS II electron transport rate `[μmol e⁻ m⁻² s⁻¹]`"
    J_P680_p::FT = 0
    "PS I electron transport rate `[μmol e⁻ m⁻² s⁻¹]`"
    J_P700_a::FT = 0
    "Rubisco limited PS I electron transport rate `[μmol e⁻ m⁻² s⁻¹]`"
    J_P700_c::FT = 0
    "Light limited PS I electron transport rate `[μmol e⁻ m⁻² s⁻¹]`"
    J_P700_j::FT = 0
    "Product limited PS I electron transport rate `[μmol e⁻ m⁻² s⁻¹]`"
    J_P700_p::FT = 0

    # Well watered condition values (to use with β function over PS)
    "Well watered maximal electron transport rate at 298.15 K `[μmol m⁻² s⁻¹]`"
    Jmax25WW::FT = Jmax25
    "Well watered respiration rate at 298.15 K `[μmol m⁻² s⁻¹]`"
    Rd25WW::FT = Rd25
    "Well watered maximal carboxylation rate at 298.15 K `[μmol m⁻² s⁻¹]`"
    Vcmax25WW::FT = Vcmax25
    "Well watered maximal PEP carboxylation rate at 298.15 K `[μmol m⁻² s⁻¹]`"
    Vpmax25WW::FT = Vpmax25

    # Fluorescence related
    "Total efficiency, incl. photorespiration `[mol CO₂ mol⁻¹ e-]`"
    e2c::FT = 1 / 6
    "dark adapted yield (`Kp=0`)"
    Fm::FT = 0
    "light adapted yield (`Kp=0`)"
    Fm′::FT = 0
    "dark-adapted fluorescence yield (`Kp=max`)"
    Fo::FT = 0
    "light-adapted fluorescence yield in the dark (`Kp=max`)"
    Fo′::FT = 0
    "Actual electron transport rate `[μmol m⁻² s⁻¹]`"
    Ja::FT = 0
    "Non-Photochemical quenching "
    NPQ::FT = 0
    "Photochemical quenching"
    qQ::FT = 0
    "energy quenching"
    qE::FT = 0
    "PSII yield"
    φ::FT = 0
    "Steady-state (light-adapted) yield (aka Fs)"
    φs::FT = 0

    # Environment related
    "Absorbed photosynthetic active radiation `[μmol m⁻² s⁻¹]`"
    APAR::FT = 100
end

max_stars_repo_name: climate-machine/Land


text:
const SUPPORTED_ENVS = (:section, :columns)

# NOTE: \begin{section} or \begin{:section} considered equiv.

function lx_begin(com, _)
    content = strip(proc(com))
    args, kwargs = lxargs(content, "begin")
    length(args) == 1 ||
        error("Expected a single argument such as `:section`.")
    env = args[1]
    env isa Symbol ||
        error("Expected a symbol as first argument such  as `:section`.")
    env in SUPPORTED_ENVS ||
        error("Unknown environment: $env.")
    starter = Symbol("_begin_$env")
    return eval(starter)(; kwargs...)
end

function lx_end(com, _)
    env = Symbol(strip(strip(proc(com)), ':'))
    env in SUPPORTED_ENVS ||
        error("Unknown environment: $env.")
    closer = Symbol("_end_$env")
    return eval(closer)()
end

#
# COLUMNS
#

_begin_columns() = html("""

    """)

_end_columns() = html("""

    """)

lx_column(com, _) = "@@col $(proc(com)) @@"

#
# SECTION
#

function _begin_section(; title="", name=title)
    id = F.refstring(name)
    pair = (id => name)
    if pair in F.locvar("sections")
        empty!(F.LOCAL_VARS["sections"].first)
    end
    push!(F.LOCAL_VARS["sections"].first, pair)
    counter = F.globvar("section_counter")
    F.set_var!(F.GLOBAL_VARS, "section_counter", counter+1)
    class = ""
    if iseven(counter)
        class = "class=\"section-bg-color\""
    end
    return html(
        """

                  $title

        """)
end

function _end_section()
    return html(
        """

        """)
end


text:
# from Pkg.jl
Base.@kwdef mutable struct MiniProgressBar
    max::Int = 1
    header::String = ""
    color::Symbol = :white
    width::Int = 32
    current::Int = 0
    prev::Int = 0
    has_shown::Bool = false
    time_shown::Float64 = 0.0
    percentage::Bool = true
    indent::Int = 4
end

const NONINTERACTIVE_TIME_GRANULARITY = Ref(2.0)
const PROGRESS_BAR_PERCENTAGE_GRANULARITY = Ref(0.1)

function showprogress(io::IO, p::MiniProgressBar)
    if p.max == 0
        perc = 0.0
        prev_perc = 0.0
    else
        perc = p.current / p.max * 100
        prev_perc = p.prev / p.max * 100
    end
    if !isinteractive()
        t = time()
        if p.has_shown && (t - p.time_shown) < NONINTERACTIVE_TIME_GRANULARITY[]
            return
        end
        p.time_shown = t
    end
    p.prev = p.current
    p.has_shown = true
    n_filled = ceil(Int, p.width * perc / 100)
    n_left = p.width - n_filled
    print(io, " "^p.indent)
    printstyled(io, p.header, color=p.color, bold=true)
    print(io, " [")
    print(io, "="^n_filled, ">")
    print(io, " "^n_left, "]  ", )
    if p.percentage
        @printf io "%2.1f %%" perc
    else
        print(io, p.current, "/",  p.max)
    end
    print(io, "\r")
end

function end_progress(io, p::MiniProgressBar)
    ansi_enablecursor = "\e[?25h"
    ansi_clearline = "\e[2K"
    print(io, ansi_enablecursor)
    print(io, ansi_clearline)
end

max_issues_repo_name: matthieugomez/FixedEffects

max_forks_repo_name: matthieugomez/FixedEffects


text:
# LinearAlgebra's definition of `det`:
#=
function det(A::AbstractMatrix{T}) where T
    if istriu(A) || istril(A) # is triangular upper or is triangular lower
        S = typeof((one(T)*zero(T) + zero(T))/one(T))
        return convert(S, det(UpperTriangular(A)))
    end
    return det(lu(A; check = false)) # LU factorisation: https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/lu.jl
end

det(x::Number) = x

# https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/lu.jl#L442
function det(F::LU{T}) where T
    n = checksquare(F)
    issuccess(F) || return zero(T)
    P = one(T)
    c = 0
    @inbounds for i = 1:n
        P *= F.factors[i,i]
        if F.ipiv[i] != i
            c += 1
        end
    end
    s = (isodd(c) ? -one(T) : one(T))
    return P * s
end
=#

function _first_cayley_hyperdet(A::AbstractArray{T, N}) where {T <: Real, N} # usually denoted "det₀"

end

function _second_cayley_hyperdet(A::AbstractArray{T, N}) where {T <: Real, N} # usually denoted "Det."

end

function LinearAlgebra.det(A::AbstractArray{T, N}) where {T <: Real, N}
    N < 3 && return det(A)
    iseven(ndims(A)) && return _first_cayley_hyperdet(A)
    # size(A) == ntuple(_ -> 2, ndims(A)) && return _first_cayley_hyperdet(A)

    return "hello"
end

# LinearAlgebra's definition of `det`:
#=
# https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/dense.jl
function inv(A::StridedMatrix{T}) where T
    checksquare(A)
    S = typeof((one(T)*zero(T) + one(T)*zero(T))/one(T))
    AA = convert(AbstractArray{S}, A)
    if istriu(AA)
        Ai = triu!(parent(inv(UpperTriangular(AA))))
    elseif istril(AA)
        Ai = tril!(parent(inv(LowerTriangular(AA))))
    else
        Ai = inv!(lu(AA))
        Ai = convert(typeof(parent(Ai)), Ai)
    end
    return Ai
end

=#
# code lowered:
#=
julia> @code_typed inv(A)
CodeInfo(
1 ── %1   = Base.arraysize(A, 1)::Int64
│    %2   = Base.arraysize(A, 2)::Int64
│    %3   = (%1 === %2)::Bool
└───        goto #3 if not %3
2 ──        goto #4
3 ── %6   = Base.arraysize(A, 1)::Int64
│    %7   = Base.arraysize(A, 2)::Int64
│    %8   = Core.tuple(%6, %7)::Tuple{Int64,Int64}
│    %9   = invoke Base.print_to_string("matrix is not square: dimensions are "::String, %8::Vararg{Any,N} where N)::String
│    %10  = %new(Base.DimensionMismatch, %9)::DimensionMismatch
│           LinearAlgebra.throw(%10)::Union{}
└───        $(Expr(:unreachable))::Union{}
4 ┄─ %13  = invoke AbstractArray{Float64,2}(_2::Array{Int8,2})::Array{Float64,2}
│    %14  = LinearAlgebra.istriu::Core.Compiler.Const(LinearAlgebra.istriu, false)
│    %15  = invoke %14(%13::Array{Float64,2}, 0::Int64)::Bool
└───        goto #23 if not %15
5 ──        Base.arraysize(%13, 1)::Int64
│           Base.arraysize(%13, 2)::Int64
│    %19  = Base.or_int(false, false)::Bool
│    %20  = Base.or_int(%19, false)::Bool
│    %21  = Base.not_int(%20)::Bool
└───        goto #7 if not %21
6 ──        goto #8
7 ── %24  = %new(Core.ArgumentError, "offset arrays are not supported but got an array with index other than 1")::ArgumentError
│           Base.throw(%24)::Union{}
└───        $(Expr(:unreachable))::Union{}
8 ┄─ %27  = Base.arraysize(%13, 1)::Int64
│    %28  = Base.arraysize(%13, 2)::Int64
│    %29  = (%27 === %28)::Bool
└───        goto #10 if not %29
9 ──        goto #11
10 ─ %32  = Base.arraysize(%13, 1)::Int64
│    %33  = Base.arraysize(%13, 2)::Int64
│    %34  = Core.tuple(%32, %33)::Tuple{Int64,Int64}
│    %35  = invoke Base.print_to_string("matrix is not square: dimensions are "::String, %34::Vararg{Any,N} where N)::String
│    %36  = %new(Base.DimensionMismatch, %35)::DimensionMismatch
│           LinearAlgebra.throw(%36)::Union{}
└───        $(Expr(:unreachable))::Union{}
11 ┄        goto #12
12 ─        goto #13
13 ─ %41  = Base.arraysize(%13, 1)::Int64
│    %42  = Base.arraysize(%13, 1)::Int64
│    %43  = LinearAlgebra.I::UniformScaling{Bool}
│    %44  = Core.tuple(%41, %42)::Tuple{Int64,Int64}
│    %45  = invoke Array{Float64,2}(%43::UniformScaling{Bool}, %44::Tuple{Int64,Int64})::Array{Float64,2}
│    %46  = LinearAlgebra.LAPACK.trtrs!::typeof(LinearAlgebra.LAPACK.trtrs!)
│    %47  = invoke %46('U'::Char, 'N'::Char, 'N'::Char, %13::Array{Float64,2}, %45::Array{Float64,2})::Array{Float64,2}
│           Base.arraysize(%47, 1)::Int64
│           Base.arraysize(%47, 2)::Int64
│    %50  = Base.or_int(false, false)::Bool
│    %51  = Base.or_int(%50, false)::Bool
│    %52  = Base.not_int(%51)::Bool
└───        goto #15 if not %52
14 ─        goto #16
15 ─ %55  = %new(Core.ArgumentError, "offset arrays are not supported but got an array with index other than 1")::ArgumentError
│           Base.throw(%55)::Union{}
└───        $(Expr(:unreachable))::Union{}
16 ┄ %58  = Base.arraysize(%47, 1)::Int64
│    %59  = Base.arraysize(%47, 2)::Int64
│    %60  = (%58 === %59)::Bool
└───        goto #18 if not %60
17 ─        goto #19
18 ─ %63  = Base.arraysize(%47, 1)::Int64
│    %64  = Base.arraysize(%47, 2)::Int64
│    %65  = Core.tuple(%63, %64)::Tuple{Int64,Int64}
│    %66  = invoke Base.print_to_string("matrix is not square: dimensions are "::String, %65::Vararg{Any,N} where N)::String
│    %67  = %new(Base.DimensionMismatch, %66)::DimensionMismatch
│           LinearAlgebra.throw(%67)::Union{}
└───        $(Expr(:unreachable))::Union{}
19 ┄        goto #20
20 ─        goto #21
21 ─        goto #22
22 ─ %73  = invoke LinearAlgebra.triu!(%47::Array{Float64,2}, 0::Int64)::Array{Float64,2}
└───        goto #56
23 ─ %75  = LinearAlgebra.istril::Core.Compiler.Const(LinearAlgebra.istril, false)
│    %76  = invoke %75(%13::Array{Float64,2}, 0::Int64)::Bool
└───        goto #42 if not %76
24 ─        Base.arraysize(%13, 1)::Int64
│           Base.arraysize(%13, 2)::Int64
│    %80  = Base.or_int(false, false)::Bool
│    %81  = Base.or_int(%80, false)::Bool
│    %82  = Base.not_int(%81)::Bool
└───        goto #26 if not %82
25 ─        goto #27
26 ─ %85  = %new(Core.ArgumentError, "offset arrays are not supported but got an array with index other than 1")::ArgumentError
│           Base.throw(%85)::Union{}
└───        $(Expr(:unreachable))::Union{}
27 ┄ %88  = Base.arraysize(%13, 1)::Int64
│    %89  = Base.arraysize(%13, 2)::Int64
│    %90  = (%88 === %89)::Bool
└───        goto #29 if not %90
28 ─        goto #30
29 ─ %93  = Base.arraysize(%13, 1)::Int64
│    %94  = Base.arraysize(%13, 2)::Int64
│    %95  = Core.tuple(%93, %94)::Tuple{Int64,Int64}
│    %96  = invoke Base.print_to_string("matrix is not square: dimensions are "::String, %95::Vararg{Any,N} where N)::String
│    %97  = %new(Base.DimensionMismatch, %96)::DimensionMismatch
│           LinearAlgebra.throw(%97)::Union{}
└───        $(Expr(:unreachable))::Union{}
30 ┄        goto #31
31 ─        goto #32
32 ─ %102 = Base.arraysize(%13, 1)::Int64
│    %103 = Base.arraysize(%13, 1)::Int64
│    %104 = LinearAlgebra.I::UniformScaling{Bool}
│    %105 = Core.tuple(%102, %103)::Tuple{Int64,Int64}
│    %106 = invoke Array{Float64,2}(%104::UniformScaling{Bool}, %105::Tuple{Int64,Int64})::Array{Float64,2}
│    %107 = LinearAlgebra.LAPACK.trtrs!::typeof(LinearAlgebra.LAPACK.trtrs!)
│    %108 = invoke %107('L'::Char, 'N'::Char, 'N'::Char, %13::Array{Float64,2}, %106::Array{Float64,2})::Array{Float64,2}
│           Base.arraysize(%108, 1)::Int64
│           Base.arraysize(%108, 2)::Int64
│    %111 = Base.or_int(false, false)::Bool
│    %112 = Base.or_int(%111, false)::Bool
│    %113 = Base.not_int(%112)::Bool
└───        goto #34 if not %113
33 ─        goto #35
34 ─ %116 = %new(Core.ArgumentError, "offset arrays are not supported but got an array with index other than 1")::ArgumentError
│           Base.throw(%116)::Union{}
└───        $(Expr(:unreachable))::Union{}
35 ┄ %119 = Base.arraysize(%108, 1)::Int64
│    %120 = Base.arraysize(%108, 2)::Int64
│    %121 = (%119 === %120)::Bool
└───        goto #37 if not %121
36 ─        goto #38
37 ─ %124 = Base.arraysize(%108, 1)::Int64
│    %125 = Base.arraysize(%108, 2)::Int64
│    %126 = Core.tuple(%124, %125)::Tuple{Int64,Int64}
│    %127 = invoke Base.print_to_string("matrix is not square: dimensions are "::String, %126::Vararg{Any,N} where N)::String
│    %128 = %new(Base.DimensionMismatch, %127)::DimensionMismatch
│           LinearAlgebra.throw(%128)::Union{}
└───        $(Expr(:unreachable))::Union{}
38 ┄        goto #39
39 ─        goto #40
40 ─        goto #41
41 ─ %134 = invoke LinearAlgebra.tril!(%108::Array{Float64,2}, 0::Int64)::Array{Float64,2}
└───        goto #56
42 ─ %136 = $(Expr(:foreigncall, :(:jl_array_copy), Ref{Array{Float64,2}}, svec(Any), 0, :(:ccall), :(%13)))::Array{Float64,2}
│    %137 = LinearAlgebra.lu!::typeof(lu!)
│    %138 = Base.sle_int(1, 1)::Bool
└───        goto #44 if not %138
43 ─ %140 = Base.sle_int(1, 0)::Bool
└───        goto #45
44 ─        nothing::Nothing
45 ┄ %143 = φ (#43 => %140, #44 => false)::Bool
└───        goto #47 if not %143
46 ─        invoke Base.getindex(()::Tuple, 1::Int64)::Union{}
└───        $(Expr(:unreachable))::Union{}
47 ┄        goto #48
48 ─        goto #49
49 ─        goto #50
50 ─        goto #51
51 ─ %151 = invoke LinearAlgebra.:(var"#lu!#132")(true::Bool, %137::typeof(lu!), %136::Array{Float64,2}, $(QuoteNode(Val{true}()))::Val{true})::LU{Float64,Array{Float64,2}}
└───        goto #52
52 ─        goto #53
53 ─        goto #54
54 ─        goto #55
55 ─ %156 = LinearAlgebra.LAPACK.getri!::typeof(LinearAlgebra.LAPACK.getri!)
│    %157 = LinearAlgebra.getfield(%151, :factors)::Array{Float64,2}
│           Base.arraysize(%157, 1)::Int64
│           Base.arraysize(%157, 2)::Int64
│    %160 = LinearAlgebra.getfield(%151, :factors)::Array{Float64,2}
│    %161 = LinearAlgebra.getfield(%151, :factors)::Array{Float64,2}
│           Base.arraysize(%161, 1)::Int64
│           Base.arraysize(%161, 2)::Int64
│    %164 = LinearAlgebra.getfield(%151, :ipiv)::Array{Int64,1}
└─── %165 = invoke %156(%160::Array{Float64,2}, %164::Array{Int64,1})::Array{Float64,2}
56 ┄ %166 = φ (#22 => %73, #41 => %134, #55 => %165)::Array{Float64,2}
└───        return %166
) => Array{Float64,2}
=#

#=
# *(A::AbstractArray{T,2} where T, B::AbstractArray{T,2} where T) in LinearAlgebra at /Applications/Julia-1.5.app/Contents/Resources/julia/share/julia/stdlib/v1.5/LinearAlgebra/src/matmul.jl:151
function (*)(A::AbstractMatrix, B::AbstractMatrix)
    TS = promote_op(matprod, eltype(A), eltype(B))
    mul!(similar(B, TS, (size(A,1), size(B,2))), A, B)
end
=#

function get_I(S::NTuple{N, T}) where {T <: Integer, N}
    A = spzeros(S...)
    for i in CartesianIndices(A)
        if allequal(i.I)
            A[i] = 1
        end
    end

    return A
end

function _promote_I(S::NTuple{N, T}) where {T  <: Integer, N}
    return T.(reshape(repeat(I(first(S)), inner = ntuple(i -> i == 2 ? first(S) : 1, first(S))), ntuple(_ -> first(S), N)))
end

function LinearAlgebra.I(S::NTuple{N, T}; dims::Int = 2) where {T <: Integer, N}
    allequal(S) || (s = rand(S); throw(DimensionMismatch("""
        The identity array should be square.  Size received: $(size(A)).  Example size: $(ntuple(_ -> s, length(A))).
    """)))

    dims == 2 && return _promote_I(S)
    return I(S, dims = (dims - 1))

    # return reshape(repeat(I(first(S)), inner = ntuple(i -> i == 2 ? first(S) : 1, first(S))), ntuple(_ -> first(S), N))

    # reshape(repeat(I(3), inner = (1, 3, 1)), (3, 3, 3))

    # cat(Iterators.repeated(I(3), 3)..., dims = 1)
    #
    # sz = ntuple(d -> d == dims ? n : size(M, d), max(N, dims))
    # return cat(M, fill(fill_elem, sz); dims = dims)
end
LinearAlgebra.I(A::Integer...) = LinearAlgebra.I(A)

LinearAlgebra.I(N::Integer...) =

function _inv_nD(A::AbstractArray{T, N}) where {T <: Real, N}
    issquare(A) || (s = rand(size(A)); throw(DimensionMismatch("""
        The inverse is defined as the matrix required to multiply with the input matrix to produce an identity matrix.  The identity matrix is square, so you need to ensure every dimension of the input array has the same size.  Size received: $(size(A)).  Example size: $(ntuple(_ -> s, ndims(A))).
    """)))
end

function LinearAlgebra.inv(A::AbstractArray{T, N}) where {T <: Real, N}
    N < 3 && return inv(A)
    return _inv_nD(A)
end

import Base: *
@generated function *(A::AbstractArray{T, N}, B::AbstractArray{T, N}) where {T, N}
    N < 3 && return :(A * B)
    C = transpose(B)
    return :(@matmul M[i,j] := sum(k) A[i, k] * C[j, k])
end


text:
module JToric

# use other Julia packages
using Pkg
using GAP
using CapAndHomalg
using JToric

"""
Initializing function for 'JToric'.
"""
# initialization
function __init__()
    # load necessary gap packages
    GAP.Globals.LoadPackage( GapObj( "JConvex" ) )
    GAP.Globals.LoadPackage( GapObj( "ToricVarieties" ) )

    # inform that JToric has been loaded
    print("Welcome to JToric ")
    printstyled("$version\n", color = :green)
    println("Enjoy this software to perform computations on toric geometry!")
    #println("Type: ?JToric for more information")
end

"""
    JToric.version

The version number of the loaded `CapAndHomalg`.
Please mention this number in any bug report.
"""
global version = 0

if VERSION >= v"1.4"
    deps = Pkg.dependencies()
    if Base.haskey(deps, Base.UUID("9bfa4af0-0a41-40a3-86ae-7e8f6cc9e7ab"))
        ver = Pkg.dependencies()[Base.UUID("9bfa4af0-0a41-40a3-86ae-7e8f6cc9e7ab")]
        if occursin("/dev/", ver.source)
            version = VersionNumber("$(ver.version)-dev")
        else
            version = VersionNumber("$(ver.version)")
        end
    else
        version = "not installed"
    end
end

# include files
include("AttributesAndMethods.jl")

end

max_stars_repo_name: fingolfin/JToric

max_issues_repo_name: fingolfin/JToric

max_forks_repo_name: fingolfin/JToric


text:
macro load_flux_fmi()
    @info "Loading differentiable compatibility to \u001b[3m\u001b[34;1mFlux.jl\u001b[0m\n\n          \u001b[34;1mhttps://github.com/FluxML/Flux.jl\n "
    expr = quote
        using Zygote
        using Flux
        using Flux: Chain, Dense, update!

        macro track(expr)
        end
    end

    expr = MacroTools.prewalk(unblock ∘ rmlines, expr)
    esc(expr)
end


text:
using Test
using WordTokenizers

@testset "Poorman's Tokenizer" begin
	@test poormans_tokenize("a b c") == ["a", "b", "c"]
	@test poormans_tokenize("a's b c") == ["as", "b", "c"]
	@test poormans_tokenize("a b- c") == ["a", "b", "c"]
end

@testset "Punctuation Space Tokenizer" begin
	@test punctuation_space_tokenize("a b c") == ["a", "b", "c"]
	@test punctuation_space_tokenize("a's b c") == ["a's", "b", "c"]
	@test punctuation_space_tokenize("a b- c") == ["a", "b", "c"]
end


text:
using Test
using JSExpr

@testset "juliaisms" begin
    @testset "range expression" begin
        expr_str = string(@js(
            for i in 2:11
                console.log(i)
            end
        ))
        @test occursin("new Array(10)", expr_str)
        @test occursin("console.log(i);", expr_str)

        @test_throws ErrorException JSExpr.crawl(:(1.5:3.5))
        @test_throws ErrorException JSExpr.crawl(:(1:2:10))
    end
end


text:
# SPDX-License-Identifier: X11
# 2020-10-18 (Timed)
# Back and Forth (300pt)

sx, sy, tx, ty = parse.(Int, split(readline()))

xd = tx - sx
yd = ty - sy

println('R'^xd *
        'U'^yd *
        'L'^xd *
        'D'^yd *
        'D' *
        'R'^(xd + 1) *
        'U'^(yd + 1) *
        'L' *
        'U' *
        'L'^(xd + 1) *
        'D'^(yd + 1) *
        'R')

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
mutable struct REPLWindow <: GtkWindow

    handle::Ptr{Gtk.GObject}
    console_manager
    statusBar

    function REPLWindow()

        w = GtkWindow("GtkREPL.jl - v$(VERSION)", 800, 600)

        signal_connect(main_window_key_press_cb, w, "key-press-event", Cint, (Ptr{Gtk.GdkEvent}, ), false)
        signal_connect(main_window_quit_cb, w, "delete-event", Cint, (Ptr{Gtk.GdkEvent}, ), false)

        style_css(w, "window, view, textview, buffer, text {
            font-family: Monaco, Consolas, Courier, monospace;
            margin:0px;
            font-size:$(fontsize)px;
          }"
        )

        n = new(w.handle)
        Gtk.gobject_move_ref(n, w)
    end

end

function init!(w::REPLWindow, console_mng)

    main_window.console_manager = console_mng
    main_window.statusBar = GtkStatusbar()
    set_gtk_property!(main_window.statusBar, :margin, 2)
    push!(main_window.statusBar, "main", "Julia $VERSION")

    main_window |>
        ((mainVbox = GtkBox(:v)) |>
            console_mng |>
            main_window.statusBar
        )
end

console_manager(main_window::T) where T<:GtkWindow = main_window.console_manager

# hooks for external process
on_command_done(main_window::T, console) where T<:GtkWindow = nothing
on_console_mng_switch_page(cm, c) = nothing

@guarded (PROPAGATE) function main_window_key_press_cb(widgetptr::Ptr, eventptr::Ptr, user_data)
    event = convert(Gtk.GdkEvent, eventptr)
    main_window = convert(GtkWindow, widgetptr)

    return PROPAGATE
end

@guarded (PROPAGATE) function main_window_quit_cb(widgetptr::Ptr, eventptr::Ptr, user_data)
    main_window = convert(GtkWindow, widgetptr)
    return PROPAGATE
end


text:
using JudiLing # our package
using CSV # read csv files
using DataFrames # parse data into dataframes

mkpath(joinpath(@__DIR__, "data"))
download("https://osf.io/2ejfu/download", joinpath(@__DIR__, "data", "latin.csv"))
# load latin file
latin = DataFrame(CSV.File(joinpath(@__DIR__, "data", "latin.csv")))
display(latin)

# create C matrixes for training datasets
cue_obj = JudiLing.make_cue_matrix(
    latin,
    grams = 3,
    target_col = :Word,
    tokenized = false,
    keep_sep = false,
)

# retrieve dim of C
# we set the S matrixes as the same dimensions
n_features = size(cue_obj.C, 2)
S = JudiLing.make_S_matrix(
    latin,
    ["Lexeme"],
    ["Person", "Number", "Tense", "Voice", "Mood"],
    ncol = n_features,
)

# we use cholesky function to calculate mapping G from S to C
G = JudiLing.make_transform_matrix(S, cue_obj.C)

# we calculate Chat matrixes by multiplying S and G
Chat = S * G
@show JudiLing.eval_SC(Chat, cue_obj.C)

# we calculate F as we did for G
F = JudiLing.make_transform_matrix(cue_obj.C, S)

# we calculate Shat as we did for Chat
Shat = cue_obj.C * F
@show JudiLing.eval_SC(Shat, S)

# here we only use a adjacency matrix as we got it from training dataset
A = cue_obj.A

# we calculate how many timestep we need for learn_paths and huo function
max_t = JudiLing.cal_max_timestep(latin, :Word)

# we calculate learn_paths and build_paths function
res_learn, gpi_learn = JudiLing.learn_paths(
    latin,
    latin,
    cue_obj.C,
    S,
    F,
    Chat,
    A,
    cue_obj.i2f,
    cue_obj.f2i, # api changed in 0.3.1
    check_gold_path = true,
    gold_ind = cue_obj.gold_ind,
    Shat_val = Shat,
    max_t = max_t,
    max_can = 10,
    grams = 3,
    threshold = 0.05,
    tokenized = false,
    keep_sep = false,
    target_col = :Word,
    verbose = true,
)

acc_learn = JudiLing.eval_acc(res_learn, cue_obj.gold_ind, verbose = false)

println("Acc for learn: $acc_learn")

res_build = JudiLing.build_paths(
    latin,
    cue_obj.C,
    S,
    F,
    Chat,
    A,
    cue_obj.i2f,
    cue_obj.gold_ind,
    max_t = max_t,
    n_neighbors = 3,
    verbose = true,
)

acc_build = JudiLing.eval_acc(res_build, cue_obj.gold_ind, verbose = false)

println("Acc for build: $acc_build")

# you can save results into csv files or dfs
JudiLing.write2csv(
    res_learn,
    latin,
    cue_obj,
    cue_obj,
    "latin_learn_res.csv",
    grams = 3,
    tokenized = false,
    sep_token = nothing,
    start_end_token = "#",
    output_sep_token = "",
    path_sep_token = ":",
    target_col = :Word,
    root_dir = @__DIR__,
    output_dir = "latin_out",
)

df_learn = JudiLing.write2df(
    res_learn,
    latin,
    cue_obj,
    cue_obj,
    grams = 3,
    tokenized = false,
    sep_token = nothing,
    start_end_token = "#",
    output_sep_token = "",
    path_sep_token = ":",
    target_col = :Word,
)

JudiLing.write2csv(
    res_build,
    latin,
    cue_obj,
    cue_obj,
    "latin_build_res.csv",
    grams = 3,
    tokenized = false,
    sep_token = nothing,
    start_end_token = "#",
    output_sep_token = "",
    path_sep_token = ":",
    target_col = :Word,
    root_dir = @__DIR__,
    output_dir = "latin_out",
)

df_build = JudiLing.write2df(
    res_build,
    latin,
    cue_obj,
    cue_obj,
    grams = 3,
    tokenized = false,
    sep_token = nothing,
    start_end_token = "#",
    output_sep_token = "",
    path_sep_token = ":",
    target_col = :Word,
)

display(df_learn)
display(df_build)

# cross-validation
download("https://osf.io/2ejfu/download", joinpath(@__DIR__, "data", "latin_train.csv"))
download("https://osf.io/bm7y6/download", joinpath(@__DIR__, "data", "latin_val.csv"))

latin_train =
    DataFrame(CSV.File(joinpath(@__DIR__, "data", "latin_train.csv")))
latin_val =
    DataFrame(CSV.File(joinpath(@__DIR__, "data", "latin_val.csv")))

# create C matrices for both training and validation datasets
cue_obj_train, cue_obj_val = JudiLing.make_cue_matrix(
    latin_train,
    latin_val,
    grams = 3,
    target_col = :Word,
    tokenized = false,
    keep_sep = false,
)

# create S matrices
n_features = size(cue_obj_train.C, 2)
S_train, S_val = JudiLing.make_S_matrix(
    latin_train,
    latin_val,
    ["Lexeme"],
    ["Person", "Number", "Tense", "Voice", "Mood"],
    ncol = n_features,
)

# here we learning mapping only from training dataset
G_train = JudiLing.make_transform_matrix(S_train, cue_obj_train.C)
F_train = JudiLing.make_transform_matrix(cue_obj_train.C, S_train)

# we predict S and C for both training and validation datasets
Chat_train = S_train * G_train
Chat_val = S_val * G_train
Shat_train = cue_obj_train.C * F_train
Shat_val = cue_obj_val.C * F_train

# we evaluate them
@show JudiLing.eval_SC(Chat_train, cue_obj_train.C)
@show JudiLing.eval_SC(Chat_val, cue_obj_val.C)
@show JudiLing.eval_SC(Shat_train, S_train)
@show JudiLing.eval_SC(Shat_val, S_val)

# we can use build path and learn path
A = cue_obj_train.A
max_t = JudiLing.cal_max_timestep(latin_train, latin_val, :Word)

res_learn_train, gpi_learn_train = JudiLing.learn_paths(
    latin_train,
    latin_train,
    cue_obj_train.C,
    S_train,
    F_train,
    Chat_train,
    A,
    cue_obj_train.i2f,
    cue_obj_train.f2i, # api changed in 0.3.1
    gold_ind = cue_obj_train.gold_ind,
    Shat_val = Shat_train,
    check_gold_path = true,
    max_t = max_t,
    max_can = 10,
    grams = 3,
    threshold = 0.05,
    tokenized = false,
    sep_token = "_",
    keep_sep = false,
    target_col = :Word,
    issparse = :dense,
    verbose = true,
)

res_learn_val, gpi_learn_val = JudiLing.learn_paths(
    latin_train,
    latin_val,
    cue_obj_train.C,
    S_val,
    F_train,
    Chat_val,
    A,
    cue_obj_train.i2f,
    cue_obj_train.f2i, # api changed in 0.3.1
    gold_ind = cue_obj_val.gold_ind,
    Shat_val = Shat_val,
    check_gold_path = true,
    max_t = max_t,
    max_can = 10,
    grams = 3,
    threshold = 0.05,
    is_tolerant = true,
    tolerance = -0.1,
    max_tolerance = 2,
    tokenized = false,
    sep_token = "-",
    keep_sep = false,
    target_col = :Word,
    issparse = :dense,
    verbose = true,
)

acc_learn_train =
    JudiLing.eval_acc(res_learn_train, cue_obj_train.gold_ind, verbose = false)
acc_learn_val = JudiLing.eval_acc(res_learn_val, cue_obj_val.gold_ind, verbose = false)

res_build_train = JudiLing.build_paths(
    latin_train,
    cue_obj_train.C,
    S_train,
    F_train,
    Chat_train,
    A,
    cue_obj_train.i2f,
    cue_obj_train.gold_ind,
    max_t = max_t,
    n_neighbors = 3,
    verbose = true,
)

res_build_val = JudiLing.build_paths(
    latin_val,
    cue_obj_train.C,
    S_val,
    F_train,
    Chat_val,
    A,
    cue_obj_train.i2f,
    cue_obj_train.gold_ind,
    max_t = max_t,
    n_neighbors = 20,
    verbose = true,
)

acc_build_train =
    JudiLing.eval_acc(res_build_train, cue_obj_train.gold_ind, verbose = false)
acc_build_val = JudiLing.eval_acc(res_build_val, cue_obj_val.gold_ind, verbose = false)

@show acc_learn_train
@show acc_learn_val
@show acc_build_train
@show acc_build_val

# Once you are done, you may want to clean up the workspace
rm(joinpath(@__DIR__, "data"), force = true, recursive = true)
rm(joinpath(@__DIR__, "latin_out"), force = true, recursive = true)


text:
using Distances

if VERSION >= v"0.5-"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

tests = ["dists"]

println("Running tests ...")

for t in tests
    fn = "test_$t.jl"
    println("* $fn ...")
    include(fn)
end


text:
ENV["GKSwstype"] = "100"  # set 'GR environment' to 'no output' (for Travis CI)
using Documenter, Reachability

makedocs(
    doctest=true, # use this flag to skip doctests (saves time!)
    modules = [Reachability],
    format = :html,
    assets = ["assets/juliareach.css"],
    sitename = "Reachability.jl",
    pages = [
        "Home" => "index.md",
        "Library" => Any[
        "User interface" => "lib/interface.md",
        "Systems" => "lib/systems.md",
        "Transformations" => "lib/transformations.md",
        "Discretization" => "lib/discretize.md",
        "Distributed computations" => "lib/distributed.md"],
        "About" => "about.md"
    ]
)

deploydocs(
    repo = "github.com/JuliaReach/Reachability.jl.git",
    target = "build",
    osname = "linux",
    julia  = "0.6",
    deps = nothing,
    make = nothing
)


text:
function bfs_test()
    A = load_matrix_network("bfs_example")
    bfs(MatrixNetwork(A),1)

    A = [1 0 1 1; 0 1 1 1; 1 1 0 1; 1 1 1 0]
    (d,dt,pred) = bfs(MatrixNetwork(sparse(A)),1)
    d_ref = [0,2,1,1]
    dt_ref = [0,3,1,2]
    pred_ref = [1,3,1,1]
    if d_ref != d_ref || dt_ref !=dt || pred_ref != pred
        error("bfs failed")
    end
    return true
end


text:
export dtw, dtwdt, computeErrorFunction, accumulateErrorFunction, backtrackDistanceFunction, computeDTWerror

"""
    dtw(ref, cur, t, window, fs; dtwnorm='L2', maxlag=80, b=1, direction=1)

Returns minimum distance time lag and index in dist array, and dtw error between traces.

# Arguments
- `ref::AbstractArray`: Reference time series
- `cur::AbstractArray`: Current time series
- `t::AbstractArray`: Time axis common to both signals
- `window::AbstractArray`: Indices for a subset of t over which to measure phase lags
- `fs::Float64`: Sampling frequency
- `dtwnorm::String`: Norm used to calculate distance; effect on the unit of dtw error. (L2 or L1)
- `maxlag::Int64`: Number of maxLag id to search the distance.
- `b::Int64t`: b value to controll in distance calculation algorithm (see Mikesell et al. 2015).
- `direction::Int64`: Direction of error accumulation (1=forward, -1=backward, 0=double to smooth)

# Returns
- `dvv::Float64`: dv/v for current correlation
- `dvv_err::Float64`: Error for calculation of dv/v
- `int::Float64`: Intercept for regression calculation
- `int_err::Float64`: Error for intercept
- `dvv0::Float64`: dv/v for current correlation forced through the origin
- `dvv0_err::Float64`: Error for calculation of dvv0
"""
function dtw(ref::AbstractArray, cur::AbstractArray, t::AbstractArray, window::AbstractArray, fs::Float64;
    dtwnorm::String="L2",
    maxLag::Int64=80,
    b::Int64=1,
    direction::Int64=1,
    norm::Bool=true)

    # measure phase shifts
    stbarTime, stbar, dist, error = dtwdt(ref, cur, t, window, fs, maxLag=maxLag, b=b, direction=direction)
    # perform linear regression of dt/t=-dv/v
    dvv, dvv_err, int, int_err, dvv0, dvv0_err = dvv_lstsq(t[window], stbarTime)

    return dvv, dvv_err, int, int_err, dvv0, dvv0_err
end

"""

    dtw(ref, cur, t, window, fs, freqbands; dtwnorm='L2', maxlag=80, b=1, direction=1, norm=true)

returns minimum distance time lag and index in dist array, and dtw error between traces.

# Arguments
- `ref::AbstractArray`: Reference time series
- `cur::AbstractArray`: Current time series
- `t::AbstractArray`: Time axis common to both signals
- `window::AbstractArray`: Indices for a subset of t over which to measure phase lags
- `fs::Float64`: Sampling frequency
- `freqbands::AbstractArray`: Frequency bands over which to compute dv/v
- `dtwnorm::String`: Norm used to calculate distance; effect on the unit of dtw error. (L2 or L1)
- `maxlag::Int64`: Number of maxLag id to search the distance.
- `b::Int64t`: b value to controll in distance calculation algorithm (see Mikesell et al. 2015).
- `direction::Int64`: Direction of error accumulation (1=forward, -1=backward, 0=double to smooth)
- `norm::Bool`: Whether or not to normalize signals before dv/v

# Returns
- `freqbands::AbstractArray`: Array of frequencies where dv/v was measured
- `dvv::Float64`: dv/v for current correlation for a range of frequencies
- `dvv_err::Float64`: Error for calculation of dv/v for a range of frequencies
- `int::Float64`: Intercept for regression calculation
- `int_err::Float64`: Error for intercept
- `dvv0::Float64`: dv/v for current correlation forced through the origin
- `dvv0_err::Float64`: Error for calculation of dvv0 for a range of frequencies
"""
function dtw(ref::AbstractArray, cur::AbstractArray, t::AbstractArray, window::AbstractArray, fs::Float64, freqbands::AbstractArray;
    dtwnorm::String="L2",
    maxLag::Int64=80,
    b::Int64=1,
    direction::Int64=1,
    norm::Bool=true)

    # define sample frequency
    dt = 1/fs

    # calculate the CWT of the time series, using identical parameters for both calculations
    cwt1, sj, freqs, coi = cwt(ref, dt, minimum(freqbands), maximum(freqbands))
    cwt2, sj, freqs, coi = cwt(cur, dt, minimum(freqbands), maximum(freqbands))

    # if a frequency window is given (instead of a set of frequency bands), we assume
    # dv/v should be calculated for each frequency. We construct a 2D array of the
    # form [f1 f1; f2 f2; ...], which can be treated the same as a 2D array of frequency bands
    if ndims(freqbands)==1
        freqbands = hcat(freqs, freqs)
    end
    # number of frequency bands
    (nbands,_) = size(freqbands)

    # initialize arrays
    dvv = zeros(nbands)
    dvv_err = similar(dvv)
    int = similar(dvv)
    int_err = similar(dvv)
    dvv0 = similar(dvv)
    dvv0_err = similar(dvv)

    # loop over frequency bands
    for iband=1:nbands
        (fmin, fmax) = freqbands[iband, :]

        # get current frequencies over which we apply icwt
        # frequency checks
        if fmax < fmin
            println("Error: please ensure columns 1 and 2 are right frequency limits in freqbands!")
        else
            freq_ind = findall(f->(f>=fmin && f<=fmax), freqs)
        end

        # perform icwt
        icwt1 = icwt(cwt1[:,freq_ind], sj[freq_ind], dt)
        icwt2 = icwt(cwt2[:,freq_ind], sj[freq_ind], dt)
        wcwt1 = real.(icwt1)
        wcwt2 = real.(icwt2)

        # normalize both signals, if appropriate
        if norm
            ncwt1 = ((wcwt1 .- mean(wcwt1)) ./ std(wcwt1))
            ncwt2 = ((wcwt2 .- mean(wcwt2)) ./ std(wcwt2))
        else
            ncwt1 = wcwt1
            ncwt2 = wcwt2
        end

        # perform dv/v
        stbarTime, stbar, dist, error = dtwdt(ncwt1, ncwt2, t, window, fs, maxLag=maxLag, b=b, direction=direction)
        # perform linear regression
        dvv[iband], dvv_err[iband], int[iband], int_err[iband], dvv0[iband], dvv0_err[iband] = dvv_lstsq(t[window], stbarTime)
    end

    return freqbands, dvv, dvv_err, int, int_err, dvv0, dvv0_err
end

"""

    dtwdt(u0, u1, t, window, fs; dtwnorm='L2', maxlag=80, b=1, direction=1)

returns minimum distance time lag and index in dist array, and dtw error between traces.

# Arguments
- `u0::AbstractArray`: Time series #1
- `u1::AbstractArray`: Time series #2
- `t::AbstractArray`: Time axis common to both signals
- `window::AbstractArray`: Indices for a subset of t over which to measure phase lags
- `fs::Float64`: Sampling frequency
- `dtwnorm::String`: Norm used to calculate distance; effect on the unit of dtw error. (L2 or L1)
- `maxlag::Int64`: Maximum distance in samples to search for best phase lag
- `b::Int64t`: Value to control in distance calculation algorithm (see Mikesell et al. 2015).
- `direction::Int64`: Direction of error accumulation (1=forward, -1=backward, 0=double to smooth)

# Returns
- `stbarTime::Array{Float64,1}`: Series of time shift at t.
- `stbar::Array{Int64,1}`: Series of minimum distance index in distance array.
- `dist::Array{Int64,2}`: Distance array.
- `error::Float64`: DTW error (distance) between two time series.
"""
function dtwdt(u0::AbstractArray, u1::AbstractArray, t::AbstractArray, window::AbstractArray, fs::Float64;
    dtwnorm::String="L2",
    maxLag::Int64=80,
    b::Int64=1,
    direction::Int64=1)

    npts = length(window) # number of samples in the signal subset we consider

    # check signal length before computing distances
    if length(u0) != length(u1) error("u0 and u1 must be same length.") end

    #compute distance between traces
    err = computeErrorFunction(u0[window], u1[window], npts, maxLag, norm=dtwnorm)

    #compute distance array and backtrack index
    if direction == 1 || direction == -1
        dist  = accumulateErrorFunction(direction, err, npts, maxLag, b) # forward accumulation to make distance function
        stbar = backtrackDistanceFunction(-direction, dist, err, -maxLag, b) # find shifts
    elseif direction == 0
        #calculate double time to smooth distance array
        dist1 = accumulateErrorFunction(-1, err, npts, maxLag, b) # forward accumulation to make distance function
        dist2 = accumulateErrorFunction(1, err, npts, maxLag, b) # backwward accumulation to make distance function
        dist  = dist1 .+ dist2 .- err # add them and remove 'err' to not count twice (see Hale's paper)
        stbar = backtrackDistanceFunction(-1, dist, err, -maxLag, b)
    else
        error("direction must be +1, -1 or 0(smoothing).")
    end

    stbarTime = stbar ./ fs # convert units from samples to time
    t2 = t[window] + stbarTime # warp the time axis

    #accumulate distance in distance array to calculate dtw error
    error = computeDTWerror(err, stbar, maxLag)

    return stbarTime, stbar, dist, error
end

"""

    computerErrorFunction(u1, u0, nSample, lag, norm='L2')

Compute error function for each sample and lag, see Hale, 2013.
Dave Hale, (2013), "Dynamic warping of seismic images," GEOPHYSICS 78: S105-S115.
https://doi.org/10.1190/geo2012-0327.1

# Arguments
- `u1::AbstractArray`: Trace we intend to warp
- `u0::AbstractArray`: Reference trace to compare with
- `nSample::Int64`: Number of points to compare in the traces
- `lag::Int64`: Maximum lag in sample number to search
- `norm::String`: 'L2' or 'L1'

# Returns
- `err::AbstractArray`: 2D error function in (samples, lags)
"""
function computeErrorFunction(u1::AbstractArray, u0::AbstractArray, nSample::Int64, lag::Int64; norm::String="L2")
    if lag >= nSample
        error("computeErrorFunction:lagProblem ","lag must be smaller than nSample")
    end

    # Allocate error function variable
    err = zeros(Float64, nSample, 2 * lag + 1 )

    # initial error calculation
    for ll = -lag:lag # loop over lags
        thisLag = ll + lag + 1
        for ii = 1:nSample # loop over samples
            if ( ii + ll >= 1 && ii + ll <= nSample ) # skip corners for now, we will come back to these
                diff = u1[ii] - u0[ii + ll] # sample difference
                if norm == "L2"
                        err[ii, thisLag] = diff^2 # difference squared error
                elseif norm == "L1"
                        err[ii, thisLag] = abs(diff) # absolute value errors
                else
                    error("norm type is not defined.")
                end
            end
        end
    end

    # Now fix corners with constant extrapolation
    for ll = -lag:lag # loop over lags
        thisLag = ll + lag + 1
        for ii = 1:nSample # loop over samples
            if ( ii + ll < 1 ) # lower left corner (negative lag, early time)
                err[ii, thisLag] = err[-ll + 1, thisLag]
            elseif ( ii + ll > nSample ) # upper right corner (positive lag, late time)
                err[ii, thisLag] = err[nSample - ll, thisLag]
            end
        end
    end

    return err
end

"""

    accumulateErrorFunction(dir, err, nSample, lag, b)

# Arguments
- `dir::Int64`: accumulation direction (dir > 0 ⟹ forward in time, dir ≦ 0 ⟹ backward in time)
- `err::Array{Float64,2}`: 2D error function
- `nSample::Int64`: number of points to compare in the traces
- `lag::Int64`: maximum lag in sample number to search
- `b::Int64`: strain limit (b∈Ints, b≥1)

# Returns
- `d::AbstractArray`: 2D distance function

The function is equation 6 in Hale, 2013.
Original by Di Yang Last modified by Dylan Mikesell (25 Feb. 2015)
"""
function accumulateErrorFunction(dir::Int64, err::Array{Float64,2}, nSample::Int64, lag::Int64, b::Int64)
    nLag = (2 * lag ) + 1 # number of lags from [ -lag : +lag ]

    # allocate distance matrix
    d = zeros(Float64, nSample, nLag)

    #--------------------------------------------------------------------------
    # Setup indices based on forward or backward accumulation direction
    #--------------------------------------------------------------------------
    if dir > 0            # FORWARD
        iBegin = 1        # start index
        iEnd   = nSample  # end index
        iInc   = 1        # increment
    else                  # BACKWARD
        iBegin = nSample  # start index
        iEnd   = 1        # stop index
        iInc   = -1       # increment
    end

    #--------------------------------------------------------------------------
    # Loop through all times ii in forward or backward direction
    for ii = iBegin:iInc:iEnd
        # min/max to account for the edges/boundaries
        ji = max(1, min(nSample, ii - iInc ))     # i-1 index
        jb = max(1, min(nSample, ii - iInc * b )) # i-b index

        # loop through all lags l
        for ll = 1:nLag
            # -----------------------------------------------------------------
            # check limits on lag indices
            lMinus1 = ll - 1 # lag at l-1

            if lMinus1 < 1  # check lag index is greater than 1
                lMinus1 = 1 # make lag = first lag
            end

            lPlus1 = ll + 1 # lag at l+1

            if lPlus1 > nLag # check lag index less than max lag
                lPlus1 = nLag # D.M. version
            end
            # -----------------------------------------------------------------

            # get distance at lags (ll-1, ll, ll+1)
            distLminus1 = d[jb, lMinus1] # minus:  d( i-b, j-1 )
            distL       = d[ji, ll]      # actual: d( i-1, j   )
            distLplus1  = d[jb, lPlus1]  # plus:   d( i-b, j+1 )

            if ji != jb # equation 10 in Hale (2013)
                for kb = ji:-iInc:jb+iInc # sum errors over i-1:i-b+1
                    distLminus1 = distLminus1 + err[kb, lMinus1]
                    distLplus1  = distLplus1  + err[kb, lPlus1]
                end
            end

            # equation 6 (if b=1) or 10 (if b>1) in Hale (2013) after treating boundaries
            d[ii, ll] = err[ii, ll] + min(distLminus1, distL, distLplus1)
        end
    end
    return d
end

"""

    backtrackDistanceFunction(dir, d, err, lmin, b)

# Arguments
- `dir::Int64`: side to start minimization (dir > 0 ⟹ forward in time, dir ≦ 0 ⟹ backward in time)
- `d::Array{Float64,2}`: 2D distance function
- `err:Array{Float64,2}`: 2D error function
- `lmin::Int64`: minimum lag to search over
- `b::Int64`: strain limit (b∈Ints, b≥1)

# Returns
- `stbar::AbstractArray`: vector of integer shifts subject to |u(i)-u(i-1)| <= 1/b

See equation 2 in Hale, 2013.
Original by Di Yang and modified from DTWDT.jl by Kurama Okubo
"""
function backtrackDistanceFunction(dir::Int64, d::Array{Float64,2}, err::Array{Float64,2}, lmin::Int64, b::Int64)
    nSample = size(d,1) # number of samples
    nLag    = size(d,2) # number of lags
    stbar   = zeros(Int64, nSample) # allocate

    #--------------------------------------------------------------------------
    # Setup indices based on forward or backward accumulation direction
    #--------------------------------------------------------------------------
    if dir > 0            # FORWARD
        iBegin = 1        # start index
        iEnd   = nSample  # end index
        iInc   = 1        # increment
    else                  # BACKWARD
        iBegin = nSample  # start index
        iEnd   = 1        # stop index
        iInc   = -1       # increment
    end
    #--------------------------------------------------------------------------
    # start from the end (front or back)
    ll0 = argmin(d[iBegin,:]) # find minimum accumulated distance at front or back depending on 'dir'
    stbar[iBegin] = ll0 + lmin - 1 # absolute value of integer shift
    #--------------------------------------------------------------------------
    # move through all time samples in forward or backward direction
    ii = iBegin

    while ii != iEnd
        if ii == iBegin
            ll = ll0
        else
            ll = ll_next
        end

        # min/max for edges/boundaries
        ji = max( 1, min(nSample, ii + iInc) )
        jb = max( 1, min(nSample, ii + iInc * b) )

        # -----------------------------------------------------------------
        # check limits on lag indices

        lMinus1 = ll - 1 # lag at l-1

        if lMinus1 < 1 # check lag index is greater than 1
            lMinus1 = 1 # make lag = first lag
        end

        lPlus1 = ll + 1 # lag at l+1

        if lPlus1 > nLag # check lag index less than max lag
            lPlus1 = nLag # D.M. and D.Y. version
        end
        # -----------------------------------------------------------------
        # get distance at lags (ll-1, ll, ll+1)
        distLminus1 = d[jb, lMinus1] # minus:  d( i-b, j-1 )
        distL       = d[ji, ll]      # actual: d( i-1, j   )
        distLplus1  = d[jb, lPlus1]  # plus:   d( i-b, j+1 )

        if ji != jb # equation 10 in Hale (2013)
            for kb = ji:iInc:jb-iInc # sum errors over i-1:i-b+1
                distLminus1 = distLminus1 + err[kb, lMinus1]
                distLplus1  = distLplus1  + err[kb, lPlus1]
            end
        end

        dl = min(distLminus1, distL, distLplus1) # update minimum distance to previous sample

        if ( dl != distL ) # then ll != ll and we check forward and backward
            if ( dl == distLminus1 )
                global ll_next = lMinus1
            else # ( dl == lPlus1 )
                global ll_next = lPlus1
            end
        else
            ll_next = ll
        end

        # assume ii = ii - 1
        ii += iInc # previous time sample

        stbar[ii] = ll_next + lmin - 1 # absolute integer of lag
        # now move to correct time index, if smoothing difference over many
        # time samples using 'b'

        if ( ll_next == lMinus1 || ll_next == lPlus1 ) # check edges to see about b values
            if ( ji != jb ) # if b>1 then need to move more steps
                for kb = ji:iInc:jb - iInc
                    ii = ii + iInc # move from i-1:i-b-1
                    stbar[ii] = ll_next + lmin - 1 # constant lag over that time
                end
            end
        end
    end

    return stbar
end

"""

    computeDTWerror(Aerr, u, lag0)

Compute the accumulated error along the warping path for DTW

# Arguments
- `Aerr::Array{Float64,2}`: error matrix (see eqn. 13 in Hale, 2013)
- `u::Array{Int64,1}`: warping function
- `lag0::Int64`: value of minimum lag

# Returns
- `error::AbstractArray`: accumulatd error along warping path

Original by Dylan Mikesell and modified from DTWDT.jl by Kurama Okubo
"""
function computeDTWerror(Aerr::Array{Float64,2}, u::Array{Int64,1}, lag0::Int)
    npts = length(u)

    if size(Aerr,1) != npts
        error("Funny things with dimensions of error matrix: check inputs.")
    end

    error = 0 # initialize

    # accumulate error
    for ii = 1:npts
        idx = lag0 + 1 + u[ii] # index of lag
        error += Aerr[ii,idx] # sum error
    end

    return error
end


text:
type NullOptionletVolatilityStructure <: OptionletVolatilityStructure end

type ShiftedLognormalVolType <: VolatilityType end
type NormalVolType <: VolatilityType end

type FlatSmileSection{DC <: DayCount, V <: VolatilityType} <: AbstractSmileSection
  exerciseDate::Date
  exerciseTime::Float64
  vol::Float64
  dc::DC
  referenceDate::Date
  atmLevel::Float64
  shift::Float64
  volType::V
  isFloating::Bool
end

function FlatSmileSection{DC <: DayCount}(d::Date, vol::Float64, dc::DC, referenceDate::Date = Date(), atmLevel::Float64 = -1.0, shift::Float64 = 0.0)
  d >= referenceDate || error("Exercise Date must be greater than reference date")
  exerciseTime = year_fraction(dc, referenceDate, d)
  isFloating = referenceDate == Date()
  return FlatSmileSection{DC, ShiftedLognormalVolType}(d, exerciseTime, vol, dc, referenceDate, atmLevel, shift, ShiftedLognormalVolType(), isFloating)
end

volatility_impl(smile::FlatSmileSection, ::Float64) = smile.vol
volatility(smile::AbstractSmileSection, rate::Float64) = volatility_impl(smile, rate)

type ConstantOptionVolatility{B <: BusinessCalendar, C <: BusinessDayConvention, DC <: DayCount} <: OptionletVolatilityStructure
  settlementDays::Int
  referenceDate::Date
  calendar::B
  bdc::C
  volatility::Float64
  dc::DC
end

function ConstantOptionVolatility{B <: BusinessCalendar, C <: BusinessDayConvention, DC <: DayCount}(settlementDays::Int, calendar::B, bdc::C, volatility::Float64, dc::DC)
  today = settings.evaluation_date
  ref_date = advance(Dates.Day(settlementDays), calendar, today, bdc)
  ConstantOptionVolatility{B, C, DC}(settlementDays, ref_date, calendar, bdc, volatility, dc)
end

# Swaption Volatility structures
type ConstantSwaptionVolatility{B <: BusinessCalendar, C <: BusinessDayConvention, DC <: DayCount} <: SwaptionVolatilityStructure
  settlementDays::Int
  referenceDate::Date
  calendar::B
  bdc::C
  volatility::Quote
  dc::DC
end

# floating reference date, floating market data
ConstantSwaptionVolatility{B <: BusinessCalendar, C <: BusinessDayConvention, DC <: DayCount}(settlementDays::Int, cal::B, bdc::C, volatility::Quote, dc::DC) =
                          ConstantSwaptionVolatility{B, C, DC}(settlementDays, Date(), cal, bdc, volatility, dc)

# Local Vol Term Structure #
type LocalConstantVol{DC <: DayCount} <: LocalVolTermStructure
  referenceDate::Date
  settlementDays::Int
  volatility::Quote
  dc::DC
end

LocalConstantVol(refDate::Date, volatility::Float64, dc::DayCount) = LocalConstantVol(refDate, 0, Quote(volatility), dc)

local_vol_impl(volTS::LocalConstantVol, ::Float64, ::Float64) = volTS.volatility.value

function local_vol(volTS::LocalVolTermStructure, t::Float64, underlyingLevel::Float64)
  # TODO check stuff
  return local_vol_impl(volTS, t, underlyingLevel)
end

function black_variance(ovs::OptionletVolatilityStructure, option_date::Date, strike::Float64)
  v = calc_volatility(ovs, option_date, strike)
  t = time_from_reference(ovs, option_date)
  return v * v * t
end

function calc_volatility(ovs::OptionletVolatilityStructure, option_date::Date, strike::Float64)
  # TODO checks - see optionletvolatilitystructure.hpp, volatility
  return volatility_impl(ovs, option_date, strike)
end

volatility_impl(ovs::OptionletVolatilityStructure, option_date::Date, strike::Float64) =
              volatility_impl(ovs, time_from_reference(ovs, option_date), strike)

volatility_impl(const_opt_vol::ConstantOptionVolatility, ::Float64, ::Float64) = const_opt_vol.volatility

# Swaption Volatility methods
function swap_length{V <: SwaptionVolatilityStructure}(swapVol::V, _start::Date, _end::Date)
  result = round(float(_end - _start) / 365.25 * 12.0) # month unit
  result /= 12.0 # year unit
  return result
end

function calc_volatility{V <: SwaptionVolatilityStructure}(swapVol::V, optionDate::Date, swapLength::Float64, strike::Float64)
  # TODO add checks
  optionTime = time_from_reference(swapVol, optionDate)
  return volatility_impl(swapVol, optionTime, swapLength, strike)
end

volatility_impl(swapVol::ConstantSwaptionVolatility, ::Float64, ::Float64, ::Float64) = swapVol.volatility.value

function black_variance{V <: SwaptionVolatilityStructure}(swapVol::V, optionDate::Date, swapLength::Float64, strike::Float64)
  v = calc_volatility(swapVol, optionDate, swapLength, strike)
  optionTime = time_from_reference(swapVol, optionDate)
  return v * v * optionTime
end

function smile_section_impl(swapVol::ConstantSwaptionVolatility, optionDate::Date, ::Dates.Period)
  atmVol = swapVol.volatility.value
  return FlatSmileSection(optionDate, atmVol, swapVol.dc, reference_date(swapVol))
end

smile_section(swapVol::SwaptionVolatilityStructure, optionDate::Date, swapTenor::Dates.Period) = smile_section_impl(swapVol, optionDate, swapTenor)


text:
load("utils/req.jl")
req("dag/transforms.jl")

# -- toexpr -------------------------------------------------------------------

toexpr(ex::LiteralEx) = ex.value
toexpr(ex::SymbolEx)  = ex.name

toexpr(ex::CallEx,     args...) = expr(:call,  args...)
toexpr(ex::RefEx,      args...) = expr(:ref,   args...)
toexpr(ex::TupleEx,    args...) = expr(:tuple, args...)
toexpr(ex::EllipsisEx, args...) = expr(:(...), args...)

toexpr(ex::AssignEx,   args...) = expr(:(=), args[1:2]...)# remove dependencies

# -- UntangleContext ----------------------------------------------------------

type Untangle
    code::Vector
    Untangle() = new({})
end
typealias UntangleContext Context{Untangle}

emit(c::UntangleContext, ex) = (push(c.c.code, ex); nothing)

# -- untangle -----------------------------------------------------------------

function untangle(sink::Node)
    sink = name_fanout_nodes(sink)
    sink = name_nodes_uniquely(sink)
    raw_untangle(sink)
end

function raw_untangle(sink::Node)
    c = UntangleContext()
    value = evaluate(c, sink)
    value, c.c.code
end

untangle(c::UntangleContext, node::KnotNode) = evaluate(c, node.args)[end]
untangle(c::UntangleContext, n::Node) = toexpr(n.val, evaluate(c, n.args)...)

function evaluate(c::UntangleContext, node::AssignNode)
    emit(c, untangle(c, node))
    return (@cached evaluate(c, get_rhs(node)))
end
function evaluate(c::UntangleContext, node::Node)
    ex = untangle(c, node)
    if !is(node.name, nothing)
        emit(c, expr(:(=), node.name, ex))
        return node.name
    end
    return ex
end

# -- wrap_kernel_body ---------------------------------------------------------

function wrap_kernel_body(flat_code::Vector, indvars)
    prologue = { :(indvars=$(quoted_tuple(indvars))) }

    body = expr(:block, append(prologue, flat_code))
    for k = 1:length(indvars)
        indvar = indvars[k]
        body = expr(:for, :(($indvar) = 1:shape[$k]), body)
    end
    body
end

# -- printing -----------------------------------------------------------------

print_untangled(sink::Node) = (print_list(untangle(sink)[2]))

max_stars_repo_name: toivoh/julia-kernels

max_issues_repo_name: toivoh/julia-kernels

max_forks_repo_name: toivoh/julia-kernels


text:
module Losses

export batchlogitbinarycrossentropy, batchmse

include("basic.jl")
end

max_stars_repo_name: dhonza/FluxGoodies

max_issues_repo_name: dhonza/FluxGoodies

max_forks_repo_name: dhonza/FluxGoodies