text:
## -------------------------------------------------------------------
# Data from
# Nanchen, Annik, Alexander Schicker, and Uwe Sauer. “Nonlinear Dependency of
# Intracellular Fluxes on Growth Rate in Miniaturized Continuous Cultures of
# Escherichia Coli.” Applied and Environmental Microbiology 72, no. 2 (February 2006):
# 1164–72. https://doi.org/10.1128/AEM.72.2.1164-1172.2006.

## -------------------------------------------------------------------
# From FIG 3: Multidimensional diagram for mini-scale chemostat cultures of E. coli.
# This data was take by hand from the graph
const FIG3 = Dict()

function _populate_FIG3()

    empty!(FIG3)

    FIG3["D"] = Dict(
        "name" => "Dilution rate",
        "unit" => "1/ hr",
        "val" => [
            0.05, 0.04, 0.05, 0.06, 0.09, 0.09, 0.11, 0.11, 0.11, 0.19, 0.19,
            0.19, 0.29, 0.31, 0.31, 0.41, 0.40, 0.41, 0.51, 0.55, 0.58, 0.69, 0.75
        ]
    )
    FIG3["X"] = Dict(
        "name" => "Biomass concentration",
        "unit" => "gDW/ L",
        "val" => [
            0.23, 0.27, 0.33, 0.28, 0.28, 0.31, 0.31, 0.35, 0.35, 0.36, 0.37,
            0.37, 0.37, 0.37, 0.37, 0.39, 0.39, 0.38, 0.41, 0.41, 0.41, 0.37, 0.20
        ]
    )

end
## -------------------------------------------------------------------
# From supplementary materials TABLE 2
const TABLE2 = Dict()

function _populate_TABLE2()

    empty!(TABLE2)

    TABLE2["Biomass"] = Dict(
        "unit" => "mmol/ gDW hr",
        "val" => [0.04, 0.05, 0.06, 0.08, 0.09, 0.09, 0.09, 0.11, 0.11, 0.18, 0.18, 0.28, 0.28, 0.31, 0.31, 0.4, 0.41],
        "err" => [0.0, 0.01, 0.0, 0.01, 0.01, 0.01, 0.01, 0.0, 0.0, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.02]
    )

    TABLE2["O2 + 2NADH -> 2P/O x ATP"] = Dict(
        "unit" => "mmol/ gDW hr",
        "val" => [3.84, 5.42, 4.68, 6.21, 6.22, 6.68, 6.68, 7.83, 7.87, 9.11, 9.09, 13.58, 13.6, 14.11, 14.14, 17.19, 17.47],
        "err" => [0.84, 1.17, 1.07, 1.34, 1.34, 1.23, 1.23, 1.78, 1.77, 1.01, 1.01, 1.59, 1.59, 1.77, 1.77, 1.66, 2.24]
    )

    TABLE2["NADPH -> NADH"] = Dict(
        "unit" => "mmol/ gDW hr",
        "val" => [0.84, 1.45, 0.53, 0.29, 0.48, 0.56, 0.71, 0.55, 0.76, 0.58, 0.61, 1.82, 1.93, 1.34, 1.45, 2.58, 2.51],
        "err" => [0.78, 0.87, 0.77, 0.91, 0.92, 0.9, 0.86, 0.9, 0.92, 0.8, 0.8, 1.01, 1.01, 1.06, 1.07, 1.0, 1.23]
    )

    TABLE2["ICT + AcCoA -> MAL + FUM + NADH"] = Dict(
        "unit" => "mmol/ gDW hr",
        "val" => [0.12, 0.11, 0.61, 0.98, 0.94, 0.8, 0.81, 0.97, 0.96, 1.33, 1.32, 0.54, 0.54, 0.94, 0.94, 0.0, 0.0],
        "err" => [0.06, 0.09, 0.13, 0.18, 0.17, 0.06, 0.04, 0.19, 0.19, 0.14, 0.14, 0.17, 0.17, 0.17, 0.17, 0.0, 0.0]
    )

    TABLE2["AcCoA -> Acetate + ATP"] = Dict(
        "unit" => "mmol/ gDW hr",
        "val" => [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
        "err" => [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    )

    TABLE2["PEP + CO2 -> OAA"] = Dict(
        "unit" => "mmol/ gDW hr",
        "val" => [0.42, 0.64, 0.35, 0.5, 0.53, 0.64, 0.63, 0.66, 0.65, 0.84, 0.83, 1.68, 1.67, 1.9, 1.89, 2.59, 2.61],
        "err" => [0.08, 0.13, 0.08, 0.1, 0.1, 0.11, 0.11, 0.14, 0.14, 0.1, 0.1, 0.19, 0.19, 0.22, 0.22, 0.22, 0.34]
    )

    TABLE2["OAA + ATP -> PEP + CO2"] = Dict(
        "unit" => "mmol/ gDW hr",
        "val" => [0.35, 0.66, 0.61, 1.19, 1.17, 0.86, 0.85, 0.72, 0.71, 1.16, 1.16, 0.64, 0.64, 1.83, 1.82, 0.88, 0.9],
        "err" => [0.06, 0.12, 0.13, 0.17, 0.17, 0.12, 0.12, 0.15, 0.15, 0.14, 0.13, 0.12, 0.12, 0.21, 0.21, 0.15, 0.16]
    )

    TABLE2["MAL -> PYR + CO2 + NADH"] = Dict(
        "unit" => "mmol/ gDW hr",
        "val" => [0.1, 0.0, 0.21, 0.0, 0.0, 0.26, 0.28, 0.56, 0.55, 0.41, 0.41, 0.67, 0.67, 0.0, 0.0, 0.46, 0.42],
        "err" => [0.04, 0.0, 0.07, 0.07, 0.06, 0.13, 0.16, 0.14, 0.14, 0.12, 0.12, 0.13, 0.13, 0.01, 0.01, 0.13, 0.34]
    )

    TABLE2["MAL -> OAA + NADH"] = Dict(
        "unit" => "mmol/ gDW hr",
        "val" => [1.24, 1.84, 1.87, 2.89, 2.82, 2.59, 2.55, 2.84, 2.82, 3.52, 3.51, 3.81, 3.8, 4.99, 4.98, 4.45, 4.58],
        "err" => [0.25, 0.37, 0.39, 0.48, 0.47, 0.38, 0.38, 0.6, 0.6, 0.35, 0.35, 0.43, 0.43, 0.56, 0.56, 0.47, 0.53]
    )

    TABLE2["FUM -> MAL"] = Dict(
        "unit" => "mmol/ gDW hr",
        "val" => [1.21, 1.73, 1.46, 1.91, 1.89, 2.05, 2.02, 2.43, 2.41, 2.61, 2.6, 3.94, 3.92, 4.05, 4.04, 4.91, 5.0],
        "err" => [0.26, 0.38, 0.32, 0.45, 0.45, 0.42, 0.42, 0.57, 0.57, 0.32, 0.32, 0.51, 0.51, 0.57, 0.57, 0.53, 0.74]
    )

    TABLE2["OGA -> FUM + CO2 + 1.5*ATP + 2NADH"] = Dict(
        "unit" => "mmol/ gDW hr",
        "val" => [1.09, 1.61, 0.85, 0.93, 0.94, 1.25, 1.21, 1.46, 1.45, 1.28, 1.28, 3.4, 3.38, 3.11, 3.1, 4.91, 5.0],
        "err" => [0.27, 0.41, 0.21, 0.49, 0.49, 0.48, 0.39, 0.46, 0.46, 0.25, 0.25, 0.57, 0.57, 0.63, 0.63, 0.53, 0.74]
    )

    TABLE2["ICT -> OGA + CO2 + NADPH"] = Dict(
        "unit" => "mmol/ gDW hr",
        "val" => [1.14, 1.66, 0.92, 1.06, 1.08, 1.4, 1.36, 1.63, 1.61, 1.56, 1.56, 3.81, 3.8, 3.58, 3.56, 5.48, 5.58],
        "err" => [0.27, 0.4, 0.22, 0.46, 0.46, 0.46, 0.37, 0.46, 0.46, 0.24, 0.24, 0.56, 0.56, 0.62, 0.62, 0.53, 0.72]
    )

    TABLE2["OAA + AcCoA -> ICT"] = Dict(
        "unit" => "mmol/ gDW hr",
        "val" => [1.26, 1.77, 1.53, 2.05, 2.02, 2.2, 2.17, 2.59, 2.57, 2.88, 2.88, 4.35, 4.34, 4.52, 4.5, 5.48, 5.58],
        "err" => [0.26, 0.37, 0.33, 0.44, 0.43, 0.4, 0.4, 0.57, 0.57, 0.31, 0.31, 0.51, 0.51, 0.56, 0.56, 0.53, 0.72]
    )

    TABLE2["PYR -> AcCoA + CO2 + NADH"] = Dict(
        "unit" => "mmol/ gDW hr",
        "val" => [1.47, 1.98, 2.26, 3.24, 3.17, 3.23, 3.21, 3.82, 3.79, 4.64, 4.63, 5.57, 5.55, 6.2, 6.19, 6.42, 6.55],
        "err" => [0.27, 0.36, 0.46, 0.47, 0.46, 0.33, 0.41, 0.72, 0.71, 0.41, 0.41, 0.5, 0.5, 0.56, 0.56, 0.53, 0.7]
    )

    TABLE2["PEP -> PYR + ATP"] = Dict(
        "unit" => "mmol/ gDW hr",
        "val" => [1.45, 2.04, 2.15, 3.44, 3.37, 3.19, 3.15, 3.51, 3.47, 4.58, 4.53, 5.48, 5.47, 6.85, 6.84, 6.86, 7.03],
        "err" => [0.23, 0.34, 0.4, 0.45, 0.44, 0.36, 0.35, 0.61, 0.6, 0.35, 0.35, 0.43, 0.43, 0.57, 0.56, 0.49, 0.51]
    )

    TABLE2["PGA -> PEP"] = Dict(
        "unit" => "mmol/ gDW hr",
        "val" => [1.56, 2.07, 1.95, 2.83, 2.8, 3.06, 3.02, 3.56, 3.52, 4.42, 4.38, 6.77, 6.76, 7.2, 7.19, 8.91, 9.09],
        "err" => [0.24, 0.33, 0.35, 0.36, 0.35, 0.33, 0.32, 0.57, 0.56, 0.3, 0.3, 0.49, 0.49, 0.55, 0.55, 0.54, 0.64]
    )

    TABLE2["T3P -> PGA + ATP + NADH"] = Dict(
        "unit" => "mmol/ gDW hr",
        "val" => [1.63, 2.14, 2.04, 2.97, 2.94, 3.21, 3.17, 3.74, 3.7, 4.72, 4.68, 7.24, 7.23, 7.73, 7.71, 9.59, 9.79],
        "err" => [0.23, 0.32, 0.35, 0.35, 0.35, 0.32, 0.31, 0.57, 0.56, 0.3, 0.3, 0.49, 0.49, 0.55, 0.55, 0.54, 0.63]
    )

    TABLE2["S7P + T3P -> E4P + F6P"] = Dict(
        "unit" => "mmol/ gDW hr",
        "val" => [0.07, 0.09, 0.09, 0.1, 0.13, 0.11, 0.14, 0.13, 0.16, 0.31, 0.31, 0.39, 0.41, 0.41, 0.43, 0.5, 0.49],
        "err" => [0.02, 0.02, 0.02, 0.03, 0.03, 0.03, 0.03, 0.04, 0.04, 0.05, 0.05, 0.06, 0.06, 0.07, 0.07, 0.06, 0.06]
    )

    TABLE2["P5P + E4P -> F6P + T3P"] = Dict(
        "unit" => "mmol/ gDW hr",
        "val" => [0.05, 0.07, 0.06, 0.06, 0.09, 0.07, 0.1, 0.08, 0.11, 0.23, 0.23, 0.27, 0.29, 0.28, 0.3, 0.34, 0.33],
        "err" => [0.02, 0.02, 0.02, 0.03, 0.03, 0.03, 0.03, 0.04, 0.04, 0.05, 0.05, 0.06, 0.06, 0.07, 0.07, 0.06, 0.06]
    )

    TABLE2["2P5P -> S7P + T3P"] = Dict(
        "unit" => "mmol/ gDW hr",
        "val" => [0.07, 0.09, 0.09, 0.1, 0.13, 0.11, 0.14, 0.13, 0.16, 0.31, 0.31, 0.39, 0.41, 0.41, 0.43, 0.5, 0.49],
        "err" => [0.02, 0.02, 0.02, 0.03, 0.03, 0.03, 0.03, 0.04, 0.04, 0.05, 0.05, 0.06, 0.06, 0.07, 0.07, 0.06, 0.06]
    )

    TABLE2["F6P + ATP -> 2T3P"] = Dict(
        "unit" => "mmol/ gDW hr",
        "val" => [0.77, 1.0, 0.95, 1.4, 1.37, 1.52, 1.48, 1.77, 1.73, 2.12, 2.08, 3.31, 3.3, 3.53, 3.52, 4.44, 4.53],
        "err" => [0.11, 0.16, 0.17, 0.18, 0.18, 0.17, 0.16, 0.29, 0.28, 0.17, 0.16, 0.26, 0.26, 0.29, 0.29, 0.29, 0.33]
    )

    TABLE2["6PG -> T3P + PYR"] = Dict(
        "unit" => "mmol/ gDW hr",
        "val" => [0.05, 0.07, 0.08, 0.11, 0.12, 0.11, 0.12, 0.14, 0.15, 0.28, 0.31, 0.38, 0.38, 0.42, 0.42, 0.41, 0.46],
        "err" => [0.05, 0.06, 0.07, 0.11, 0.11, 0.12, 0.11, 0.14, 0.13, 0.14, 0.14, 0.19, 0.19, 0.21, 0.2, 0.21, 0.22]
    )

    TABLE2["G6P -> F6P"] = Dict(
        "unit" => "mmol/ gDW hr",
        "val" => [0.65, 0.85, 0.81, 1.25, 1.16, 1.35, 1.25, 1.57, 1.46, 1.58, 1.55, 2.68, 2.63, 2.86, 2.8, 3.63, 3.73],
        "err" => [0.09, 0.13, 0.14, 0.16, 0.14, 0.14, 0.13, 0.25, 0.23, 0.13, 0.12, 0.21, 0.21, 0.24, 0.23, 0.25, 0.28]
    )

    TABLE2["6PG -> P5P + CO2 + NADPH"] = Dict(
        "unit" => "mmol/ gDW hr",
        "val" => [0.21, 0.27, 0.26, 0.31, 0.39, 0.34, 0.42, 0.39, 0.5, 0.97, 0.97, 1.23, 1.29, 1.33, 1.39, 1.65, 1.63],
        "err" => [0.05, 0.07, 0.07, 0.09, 0.09, 0.09, 0.09, 0.11, 0.12, 0.14, 0.13, 0.18, 0.18, 0.2, 0.2, 0.18, 0.19]
    )

    TABLE2["G6P -> 6PG + NADPH"] = Dict(
        "unit" => "mmol/ gDW hr",
        "val" => [0.26, 0.34, 0.34, 0.42, 0.51, 0.45, 0.54, 0.53, 0.65, 1.24, 1.28, 1.61, 1.67, 1.74, 1.8, 2.06, 2.09],
        "err" => [0.05, 0.06, 0.07, 0.08, 0.08, 0.08, 0.08, 0.1, 0.11, 0.11, 0.11, 0.16, 0.16, 0.17, 0.17, 0.19, 0.2]
    )

    TABLE2["GLC + ATP -> G6P"] = Dict(
        "unit" => "mmol/ gDW hr",
        "val" => [0.91, 1.19, 1.15, 1.67, 1.67, 1.8, 1.8, 2.11, 2.12, 2.84, 2.84, 4.31, 4.31, 4.62, 4.62, 5.71, 5.84],
        "err" => [0.12, 0.17, 0.18, 0.18, 0.18, 0.16, 0.16, 0.3, 0.3, 0.15, 0.15, 0.25, 0.25, 0.29, 0.29, 0.28, 0.33]
    )

    TABLE2["uAC"] = Dict(
        "unit" => "mmol/ gDW hr",
        "val" => [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
        "err" => [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    )

    TABLE2["uGLC"] = Dict(
        "unit" => "mmol/ gDW hr",
        "val" => -[0.9, 1.2, 1.1, 1.7, 1.7, 1.8, 1.8, 2.0, 2.0, 2.8, 2.8, 4.2, 4.2, 4.6, 4.6, 5.7, 5.8],
        "err" => [0.1, 0.2, 0.2, 0.2, 0.2, 0.2, 1.2, 0.3, 0.3, 0.2, 0.2, 0.3, 0.4, 0.3, 0.3, 0.3, 0.3]
    )

    TABLE2["D"] = Dict(
        "unit" => "1/ hr",
        "val" => [0.04, 0.05, 0.06, 0.09, 0.09, 0.09, 0.09, 0.11, 0.11, 0.19, 0.19, 0.29, 0.29, 0.31, 0.31, 0.4, 0.41],
        "err" => [0.01, 0.01, 0.0, 0.01, 0.01, 0.01, 0.01, 0.0, 0.0, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.02]
    )
end

## -------------------------------------------------------------------
const DAT_IDERS = [
    "GLC + ATP -> G6P", "Biomass", "PEP -> PYR + ATP", "FUM -> MAL",
    "AcCoA -> Acetate + ATP", "O2 + 2NADH -> 2P/O x ATP", "G6P -> 6PG + NADPH",
    "T3P -> PGA + ATP + NADH", "S7P + T3P -> E4P + F6P", "MAL -> OAA + NADH",
    "PYR -> AcCoA + CO2 + NADH", "xi", "PGA -> PEP", "P5P + E4P -> F6P + T3P",
    "F6P + ATP -> 2T3P", "ICT + AcCoA -> MAL + FUM + NADH", "X", "MAL -> PYR + CO2 + NADH",
    "uGLC", "PEP + CO2 -> OAA", "D", "NADPH -> NADH", "uAC", "OGA -> FUM + CO2 + 1.5*ATP + 2NADH",
    "cGLC", "ICT -> OGA + CO2 + NADPH", "G6P -> F6P", "OAA + AcCoA -> ICT", "6PG -> T3P + PYR",
    "2P5P -> S7P + T3P", "6PG -> P5P + CO2 + NADPH", "OAA + ATP -> PEP + CO2"
]

## -------------------------------------------------------------------
# BUNDLE
const EXPS = 1:17
const BUNDLE = Dict()
function _populate_bundle()

    empty!(BUNDLE)
    _populate_TABLE2()
    _populate_FIG3()

    BUNDLE["D"] = TABLE2["D"]

    # cell concentration (X)
    # I will average values extracted from FIG 3
    # estimated error from the handed procedure of copying from fig3
    rewrite_err = 0.03 # gDW/ L
    fig3Xs = Dict()
    for (fig3D, fig3X) in zip(FIG3["D"]["val"], FIG3["X"]["val"])
        get!(fig3Xs, fig3D, [])
        push!(fig3Xs[fig3D], fig3X)
    end

    BUNDLE["X"] = Dict()
    BUNDLE["X"]["name"] = "cell concentration"
    BUNDLE["X"]["unit"] = "gDW / L"
    BUNDLE["X"]["val"] = []
    BUNDLE["X"]["err"] = []
    for D in BUNDLE["D"]["val"]
        av = mean(fig3Xs[D])
        err = length(fig3Xs[D]) > 1 ? max(std(fig3Xs[D]), rewrite_err) : rewrite_err
        push!(BUNDLE["X"]["val"], av)
        push!(BUNDLE["X"]["err"], err)
    end

    # xi
    BUNDLE["xi"] = Dict()
    BUNDLE["xi"]["name"] = "dilution specific cell concentration"
    BUNDLE["xi"]["units"] = "gCDW/ L h"
    BUNDLE["xi"]["val"] = BUNDLE["X"]["val"] ./ BUNDLE["D"]["val"]
    BUNDLE["xi"]["err"] = BUNDLE["X"]["err"]

    # cGLC
    # "M9 minimal medium with 1 g per liter glucose as
    # the limiting substrate at a pH of 7.0"
    # MM_glc = 180.156 g/mol
    # c(g/L) / MM_glc(g/mol) = (mol/L) * 1e3 = (mmol/L)
    BUNDLE["cGLC"] = Dict()
    BUNDLE["cGLC"]["name"] = "feed glc concentration"
    BUNDLE["cGLC"]["units"] = "mM"
    BUNDLE["cGLC"]["val"] = fill((1.0 / 180.156) * 1e3, length(BUNDLE["D"]["val"]))
    BUNDLE["cGLC"]["err"] = fill(0.0, length(BUNDLE["D"]["val"]))

    # rxns
    for (rxn, dat) in TABLE2
        rxn == "D" && continue
        BUNDLE[rxn] = TABLE2[rxn]
    end
end

## ------------------------------------------------------------------
# API
const msd_mets = ["GLC", "AC"]

_get_val(id, dk) = BUNDLE[string(id)][dk]
_get_val(id, dk, exp::Int) = BUNDLE[string(id)][dk][exp]
_get_val(id, dk, ref, dflt) = try; _get_val(id, dk, ref); catch err; dflt end

for fun in [:val, :err]
    dk = string(fun)
    @eval $(fun)(id) = _get_val(id, $(dk))
    @eval $(fun)(id, ref) = _get_val(id, $(dk), ref)
    @eval $(fun)(id, ref, dflt) = _get_val(id, $(dk), ref, dflt)

    for p in [:u, :c]
        pstr = string(p)
        pfun = Symbol(p, fun)
        @eval $(pfun)(id, args...) = $(fun)(string($(pstr), id), args...)
    end
end

unit(id) = BUNDLE[string(id)]["unit"]

ciD_X(id) = [cval(id, exp, 0.0) * val(:D, exp) / val(:X, exp) for exp in EXPS]
ciD_X(id, exp) = cval(id, exp, 0.0) * val(:D, exp) / val(:X, exp)


text:
using GeometricIntegrators

set_config(:tab_compensated_summation, true)

const out_dir = "../../docs/src/integrators/runge_kutta/guiding_center_4d/"
const run_id  = "fast_barely_passing_firk_comp"

include("guiding_center_4d_fast_barely_passing.jl")
include("../guiding_center_4d_firk.jl")

ext: jl


text:
# this script splits the yearlong seisio files into smaller files to
# optimize memory use per core
# For the test data, using cc_len = 2^15 seconds and cc_step = 26788 seconds,
# there are 1,177 time windows per year
# using Npersplit = 25 gives 1,177 \ 25 = 48 splits
using SeisIO, Glob, Dates

files = glob("*SAC",expanduser("~/SeisNoise_paper/DATA/YEARSAC/"))
OUTDIR = expanduser("~/SeisNoise_paper/DATA/SPLITSEIS")
cc_len, cc_step = 2^15, 26768
Npersplit = 25  # number of time-windows per split
starttime = DateTime(2019,1,1)

function splitsta(file::String,cc_len::Int,cc_step::Int,starttime::DateTime,Npersplit::Int,OUTDIR::String)

    println("Reading file ",basename(file))
    S = read_data("sac",file)
    demean!(S)
    OUTDIR = joinpath(OUTDIR,string(Npersplit))
    _,endtime  = u2d.(SeisIO.t_win(S[1].t,S[1].fs) .* 1e-6)

    if !isdir(OUTDIR)
        mkpath(OUTDIR)
    end

    # find number of splits1:
    Nsplits = convert(Int,ceil(((endtime - starttime).value / 1000 - cc_len) / cc_step / Npersplit))
    t1 = starttime
    splitNum = 1

    while t1 < endtime
        t2 = t1 + Second((Npersplit - 1) * cc_step + cc_len - 1)
        t2 = min(t2,endtime)
        Ssplit = sync(S,s=t1,t=t2,pad=true)

        splitstr = lpad(string(splitNum), Int(ceil(log10(Nsplits))), '0')
        SPLITDIR = joinpath(OUTDIR,"SPLIT" * splitstr)
        if !isdir(SPLITDIR)
            mkpath(SPLITDIR)
        end
        fileout = joinpath(SPLITDIR,Ssplit[1].id * ".sac")
        writesac(Ssplit,fname=fileout)

        # update
        t1 += Second(Npersplit * cc_step)
        splitNum += 1
    end
    return nothing
end

for file in files
    splitsta(file,cc_len,cc_step,starttime,Npersplit,OUTDIR)
end


text:
module TestShow

include("preamble.jl")

xforms = [
    Cat(),
    Count(),
    Count(2),
    Count(3, 4),
    GetIndex([0]),
    # GetIndex{true}([0]),
    Iterated(sqrt, 0.1),
    # Iterated(sqrt, 0.1, Number),
    Map(sin),
    Filter(isfinite),
    Scan(*),
    ZipSource(Filter(isfinite)),
    ZipSource(Filter(isfinite) |> Map(sin)),
    ZipSource(Filter(isfinite) |> Map(sin)) |> Map(sum),
    let xf = ZipSource(Filter(isfinite) |> Map(sin))
        xf |> Map(sum) |> xf
    end,
    let xf = Map(first) |> Map(last)
        xf = ZipSource(ZipSource(xf) |> Map(identity)) |> xf
        xf |> ZipSource(xf) |> xf
    end,
    # Zip(Map(sin), Map(cos), Map(tan)) |> Map(prod),
]

@testset "smoke test summary(xf)" begin
    for xf in xforms
        @test sprint(summary, xf) isa String
    end
end

# Note: Use `ENV["JULIA_DEBUG"] = "Main"` to enable debugging
@testset "eval(show(xf))" begin
    @testset "$(summary(xf))" for xf in xforms
        code = sprint(show, xf)
        @debug """
        show(xf) =
        $code
        """
        xf2 = include_string(@__MODULE__, code)
        @test xf == xf2
    end
end
@testset "eval(show(text/plain, xf))" begin
    @testset "$(summary(xf))" for xf in xforms
        code = sprint(show, "text/plain", xf)
        @debug """
        show("text/plain", xf) =
        $code
        """
        xf2 = include_string(@__MODULE__, code)
        @test xf == xf2
    end
end

@testset "Reduced" begin
    @test sprint(show, Reduced(1)) == "Transducers.Reduced(1)"
    @test sprint(show, Reduced(1); context = :limit => true) == "Reduced(1)"
end

@testset "Completing(+)" begin
    rf = Completing(+)
    kw = (; context = :module => Base)
    @test repr(rf; kw...) == "Transducers.Completing(+)"
    @test sprint(print, rf; kw...) == "Transducers.Completing(+)"
    @test sprint(show, rf; kw...) == "Transducers.Completing(+)"
    @test sprint(show, "text/plain", rf; kw...) == "Transducers.Completing(+)"
end

@testset "Completing{Any}(+)" begin
    rf = Completing{Any}(+)
    kw = (; context = :module => Base)
    @test repr(rf; kw...) == "Transducers.Completing{Any}(+)"
    @test sprint(print, rf; kw...) == "Transducers.Completing{Any}(+)"
    @test sprint(show, rf; kw...) == "Transducers.Completing{Any}(+)"
    @test sprint(show, "text/plain", rf; kw...) == "Transducers.Completing{Any}(+)"
end

@testset "TeeRF(min, max)" begin
    rf = TeeRF(min, max)
    kw = (; context = :module => Base)
    @test repr(rf; kw...) == "Transducers.TeeRF(min, max)"
    @test sprint(print, rf; kw...) == "Transducers.TeeRF(min, max)"
    @test sprint(show, rf; kw...) == "Transducers.TeeRF(min, max)"
    @test sprint(show, "text/plain", rf; kw...) == "Transducers.TeeRF(min, max)"
end

@testset "TeeRF{Tuple{Any,Any}}((min, max))" begin
    rf = TeeRF{Tuple{Any,Any}}((min, max))
    kw = (; context = :module => Base)
    @test repr(rf; kw...) == "Transducers.TeeRF{Tuple{Any,Any}}((min, max))"
    @test sprint(print, rf; kw...) == "Transducers.TeeRF{Tuple{Any,Any}}((min, max))"
    @test sprint(show, rf; kw...) == "Transducers.TeeRF{Tuple{Any,Any}}((min, max))"
    @test sprint(show, "text/plain", rf; kw...) ==
          "Transducers.TeeRF{Tuple{Any,Any}}((min, max))"
end

@testset "ProductRF(min, max)" begin
    rf = ProductRF(min, max)
    kw = (; context = :module => Base)
    @test repr(rf; kw...) == "Transducers.ProductRF(min, max)"
    @test sprint(print, rf; kw...) == "Transducers.ProductRF(min, max)"
    @test sprint(show, rf; kw...) == "Transducers.ProductRF(min, max)"
    @test sprint(show, "text/plain", rf; kw...) == "Transducers.ProductRF(min, max)"
end

@testset "ProductRF{Tuple{Any,Any}}((min, max))" begin
    rf = ProductRF{Tuple{Any,Any}}((min, max))
    kw = (; context = :module => Base)
    @test repr(rf; kw...) == "Transducers.ProductRF{Tuple{Any,Any}}((min, max))"
    @test sprint(print, rf; kw...) == "Transducers.ProductRF{Tuple{Any,Any}}((min, max))"
    @test sprint(show, rf; kw...) == "Transducers.ProductRF{Tuple{Any,Any}}((min, max))"
    @test sprint(show, "text/plain", rf; kw...) ==
          "Transducers.ProductRF{Tuple{Any,Any}}((min, max))"
end

end  # module


text:
module DataDrop

using JSON
using SparseArrays
using HDF5

const _sep = "-"

"""
    clean_file_name(name)

Construct a clean file name.

Replace '.', ':', ' ' (dot, colon, space) with underscores.

Example:
```
s =  clean_file_name("something or other.1=0.5_5:7")
s == "something_or_other_1=0_5_5_7"
```
"""
function clean_file_name(name)
    cleanname = replace(replace(replace(name,'.'=>'_'),':'=>'_'),' '=>'_')
    return cleanname
end

"""
    file_extension(filename)

Extract extension from the string representing a file name.

Expected result:
```
p = "matrices.dir"
file_extension(p) == ".dir"
p = "matrices.h5"
file_extension(p) == ".h5"
p = "matri.ces.dir"
file_extension(p) == ".dir"
p = "matrices."
file_extension(p) == "."
p = "matrices"
file_extension(p) == ""
```
"""
function file_extension(filename)
    f, e = splitext(filename)
    return e
end

"""
    with_extension(filename, ext)

Make sure the file name has the desired extension.

`ext` can be either with or without the leading dot.

```
with_extension("fo_1=0_5_5_7", ".ext") == "fo_1=0_5_5_7.ext"
with_extension("fo_1=0_5_5_7", "ext") == "fo_1=0_5_5_7.ext"
with_extension("fo_1=0_5_5_7.dat", "ext") == "fo_1=0_5_5_7.ext"
```
"""
function with_extension(filename, ext)
    if ext == ""
        # Remove the extension completely
        f, e = splitext(filename)
        return f
    else
        # If the file has the wrong extension, remove it
        f, e = splitext(filename)
        # Replace or fix up the extension
        if ext[1] == '.'
            ext = ext[2:end]
        end
        if match(Regex(".*\\." * ext * "\$"), f) == nothing
            f = f * "." * ext
        end
        return f
    end
end

"""
    retrieve_json(j)

Retrieve a dictionary from a JSON file named `j`.

Example:
```
d = retrieve_json("myfile.json")
```
"""
function retrieve_json(j)
    j = with_extension(j, "json")
    return open(j, "r") do file
        JSON.parse(file)
    end
end

"""
    store_json(j, d)

Store a dictionary `d` into a JSON file named `j`.

The dictionary is stored into a JSON file. If the file name `fname` is supplied
without an extension, it is appended the ".json" extension.
"""
function store_json(j, d)
    j = with_extension(j, "json")
    open(j, "w") do file
        JSON.print(file, d, 4)
    end
end

const MATRIX_UNKNOWN = -1
const MATRIX_DENSE = 0
const MATRIX_SPARSE = 1

"""
    retrieve_matrix(fname)

Retrieve a matrix from file `fname`.
"""
function retrieve_matrix(fname)
    return retrieve_matrix(fname, "")
end

"""
    empty_hdf5_file(fname)

Empty an HDF5 file. All contents will be erased.
"""
function empty_hdf5_file(fname)
    # If the file exists, delete all contents
    h5open(fname, "w") do fid
    end
end

"""
    retrieve_matrix(fname, mpath)

Retrieve a matrix from file `fname`.

The matrix data is stored under the path `mpath`.

Example:
Storing a sparse matrix as "/data" with
`store_matrix("test/matrices/matrix_SInt64.h5f", "/data", d)` gives:
```
julia> f = h5open("test/matrices/matrix_SInt64.h5f", "r")
�
└─ � data
   ├─ � I
   ├─ � J
   ├─ � V
   ├─ � matrix_type
   ├─ � ncols
   └─ � nrows
```

Example:
Storing a dense matrix under the default path ("/") with
`store_matrix("matrix_d_Float32.h5", d)` gives:
```
julia> f = h5open("matrix_d_Float32.h5", "r")
�
├─ � matrix
└─ � matrix_type
```
"""
function retrieve_matrix(fname, mpath)
    if file_extension(fname) == ""
        fname = with_extension(fname, "h5")
    end
    if mpath != "" && mpath[end:end] != "/"
        mpath = mpath * "/"
    end
    typ = MATRIX_UNKNOWN
    return h5open(fname, "r") do f
        typ = read(f, mpath * "matrix_type")
        if typ == MATRIX_DENSE
            v = read(f, mpath * "matrix")
        elseif typ == MATRIX_SPARSE
            I = read(f, mpath * "I")
            J = read(f, mpath * "J")
            V = read(f, mpath * "V")
            nrows = read(f, mpath * "nrows")
            ncols = read(f, mpath * "ncols")
            v = sparse(I, J, V, nrows, ncols)
        end
        v
    end
end

"""
    store_matrix(fname, matrix)

Store a dense matrix into the file named `fname`.
"""
function store_matrix(fname, matrix)
    store_matrix(fname, "", matrix)
end

"""
    store_matrix(fname, mpath, matrix)

Store a dense matrix under the path `mpath` into the file `fname`.
"""
function store_matrix(fname, mpath, matrix)
    if file_extension(fname) == ""
        fname = with_extension(fname, "h5")
    end
    if mpath != "" && mpath[end:end] != "/"
        mpath = mpath * "/"
    end
    h5open(fname, "cw") do file
        write(file, mpath * "matrix_type", MATRIX_DENSE)
        write(file, mpath * "matrix", matrix)
    end
end

"""
    store_matrix(fname, mpath, matrix::SparseArrays.SparseMatrixCSC{T, Int64}) where {T}

Store a sparse matrix under the path `mpath` into the file named `fname`.

The matrix is stored into an HDF5 file. If the file name `fname` is supplied
without an extension, it is appended the ".h5" extension.
"""
function store_matrix(fname, mpath, matrix::SparseArrays.SparseMatrixCSC{T, Int64}) where {T}
    I, J, V = findnz(matrix)
    if file_extension(fname) == ""
        fname = with_extension(fname, "h5")
    end
    if mpath != "" && mpath[end:end] != "/"
        mpath = mpath * "/"
    end
    h5open(fname, "cw") do file
        write(file, mpath * "matrix_type", MATRIX_SPARSE)
        write(file, mpath * "I", I)
        write(file, mpath * "J", J)
        write(file, mpath * "V", V)
        write(file, mpath * "nrows", size(matrix, 1))
        write(file, mpath * "ncols", size(matrix, 2))
    end
end

"""
    store_matrix(fname, matrix::SparseArrays.SparseMatrixCSC{T, Int64}) where {T}

Store a sparse matrix into the file named `fname`.

The matrix is stored into an HDF5 file. If the file name `fname` is supplied
without an extension, it is appended the ".h5" extension.
"""
function store_matrix(fname, matrix::SparseArrays.SparseMatrixCSC{T, Int64}) where {T}
    store_matrix(fname, "", matrix)
end

"""
    store_value(fname, vpath, num)

Store a number or a string into the file named `fname` under the name `vpath`.
"""
function store_value(fname, vpath, num)
    if file_extension(fname) == ""
        fname = with_extension(fname, "h5")
    end
    if vpath == ""
        vpath = "number"
    end
    if vpath != "" && vpath[end:end] != "/"
        vpath = vpath * "/"
    end

    h5open(fname, "cw") do file
        write(file, vpath, num)
    end
end

"""
    retrieve_value(fname, vpath)

Retrieve a number or a string known as `vpath` from the file named `fname`.
"""
function retrieve_value(fname, vpath)
    if file_extension(fname) == ""
        fname = with_extension(fname, "h5")
    end
    if vpath != "" && vpath[end:end] != "/"
        vpath = vpath * "/"
    end
    return h5open(fname, "r") do f
        v = read(f, vpath)
        v
    end
end

end # module


text:
"""
REINFORCE learner (Williams 1992)
with a suprise-modulated learning rate.

See also REINFORCElearner.jl and sarsaLmodulatedlearner.jl

Specifically for 2 actions and reward only at the end of episode!

Note
----
The policy parameters θ as well as the feature vector ϕ are
column vectors. Here we use a (states x actions) matrix for simpler access.
also, the feature vector ϕ is a one-hot vector in our case. therefore we
do not need it to keep a vector but access the entries by the
state/action pair [s_t, a_t]
"""

mutable struct REINFORCEModulated
    nr_states::Int
    nr_actions::Int
    α::Float64
    α_eff::Float64
    γ::Float64
    γ_k::Float64
    action_selection_prob_in_episode::Array{Any,1}
    states_in_episode::Array{Any,1}
    actions_in_episode::Array{Any,1}
    θ::Array{Float64, 2}
    target::Array{Float64, 1}
    modulatingfactor_in_episode::Array{Any,1}
    do_debug_print::Bool
end
function REINFORCEModulated(; nr_states = 10, nr_actions = 4, α = 0.1, γ = .9,
    do_debug_print = false)
    θ = zeros(Float64, nr_states, nr_actions)

    REINFORCEModulated(nr_states, nr_actions, α, α, γ, 1., [], [], [], θ, [0.],
                        [], do_debug_print)
end
export REINFORCEModulated


text:
filename() = "day22.txt"

function read_data()
    blocks = split(String(read(filename())), "\n\n")
    parse_block(b) = parse.(Int, split(b, "\n")[2:end])
    return parse_block.(blocks)
end

function play(p1, p2)
    while length(p1) > 0 && length(p2) > 0
        c1, c2 = popfirst!(p1), popfirst!(p2)
        if c1 > c2
            append!(p1, [c1, c2])
        else
            append!(p2, [c2, c1])
        end
    end
    return length(p1) > 0 ? p1 : p2
end

function part1()
    p1, p2 = read_data()
    hands = play(p1, p2)
    multipliers = range(length(hands), 1, step = -1)
    return sum(hands .* multipliers)
end

function play2(p1, p2)
    mem = Set()
    while length(p1) > 0 && length(p2) > 0

        # check for instant win, else remember this round
        hash((p1,p2)) ∈ mem && return (1, p1)
        push!(mem, hash((p1,p2)))

        # draw cards
        c1, c2 = popfirst!(p1), popfirst!(p2)

        # play sub game?
        if length(p1) >= c1 && length(p2) >= c2
            winner, _ = play2(copy(p1)[1:c1], copy(p2)[1:c2])
        else
            winner = c1 > c2 ? 1 : 2
        end

        if winner == 1
            append!(p1, [c1, c2])
        else
            append!(p2, [c2, c1])
        end
    end
    return length(p1) > 0 ? (1, p1) : (2, p2)
end

function part2()
    p1, p2 = read_data()
    winner, hands = play2(p1, p2)
    multipliers = range(length(hands), 1, step = -1)
    return sum(hands .* multipliers)
end


text:
# - Julia online documentation
# - Collections and Data Structures
# - General Collections
using Base.Test
@testset "General Collections" begin
c::Vector{Int64} = []
# --------------------------------------------- #
@test Base.isempty(c) == true
c = 5:-1:1
@test Base.isempty(c) == false
# --------------------------------------------- #
c = 1:5
c = Base.empty!(c)
@test Base.isempty(c) == true
# --------------------------------------------- #
c = 10:-1:1
@test Base.length(c) == 10
@test Base.endof(c) == 10
end

# ============================================= #

# ============================================= #
# --------------------------------------------- #
# ============================================= #
# --------------------------------------------- #
# ============================================= #
# --------------------------------------------- #
# ============================================= #
# --------------------------------------------- #
# ============================================= #
# --------------------------------------------- #
# ============================================= #
# --------------------------------------------- #
# ============================================= #
# --------------------------------------------- #


text:
#=
These functions are taken
  from https://github.com/mlubin/jlSimplex/blob/master/jlSimplex.jl
  and modified as needed.
=#

using GLPK

function GLPK_solve(mpsfile::String)
    lp = GLPK.Prob()
    ret = GLPK.read_mps(lp, GLPK.MPS_FILE, mpsfile)
    @assert ret == 0

    lps_opts = GLPK.SimplexParam()
    # lps_opts.presolve = GLPK.ON
    lps_opts.msg_lev = GLPK.MSG_ALL
    lps_opts.meth = GLPK.PRIMAL

    GLPK.simplex(lp, lps_opts)
end

function GLPK_read_mps(mpsfile::String)

    # Read MPS file from GLPK
    lp = GLPK.glp_create_prob()
    ret = GLPK.glp_read_mps(lp, GLPK.GLP_MPS_FILE, C_NULL, mpsfile)
    @assert ret == 0

    # Numbers of rows and columns
    nrow::Int = GLPK.glp_get_num_rows(lp)
    ncol::Int = GLPK.glp_get_num_cols(lp)
    nnz::Int = GLPK.glp_get_num_nz(lp);
    # @show nrow
    # @show ncol
    # @show nnz

    index1 = Array{Int32,1}(undef, nrow)
    coef1 = Array{Float64,1}(undef, nrow)

    c = Vector{Float64}(undef, ncol)
    xlb = Vector{Float64}(undef, ncol)
    xub = Vector{Float64}(undef, ncol)
    l = Vector{Float64}(undef, nrow)
    u = Vector{Float64}(undef, nrow)

    for i in 1:ncol
        c[i] = GLPK.glp_get_obj_coef(lp,i)
        # @show i, c[i]
        t = GLPK.glp_get_col_type(lp,i)
        # @show t
        if t == GLPK.GLP_FR
            xlb[i] = typemin(Float64)
            xub[i] = typemax(Float64)
        elseif t == GLPK.GLP_UP
            xlb[i] = typemin(Float64)
            xub[i] = GLPK.glp_get_col_ub(lp,i)
        elseif t == GLPK.GLP_LO
            xlb[i] = GLPK.glp_get_col_lb(lp,i)
            xub[i] = typemax(Float64)
        elseif t == GLPK.GLP_DB || t == GLPK.GLP_FX
            xlb[i] = GLPK.glp_get_col_lb(lp,i)
            xub[i] = GLPK.glp_get_col_ub(lp,i)
        else
            @error "Unexpected column type $t"
        end
    end

    objname = GLPK.glp_get_obj_name(lp)
    GLPK.glp_create_index(lp)
    objrow = GLPK.glp_find_row(lp,objname)

    for i in 1:nrow
        reali = i
        # if (i >= objrow)
        #     reali += 1
        # end
        t = GLPK.glp_get_row_type(lp,reali)
        if t == GLPK.GLP_UP
            l[i] = typemin(Float64)
            u[i] = GLPK.glp_get_row_ub(lp,reali)
        elseif t == GLPK.GLP_LO
            l[i] = GLPK.glp_get_row_lb(lp,reali)
            u[i] = typemax(Float64)
        elseif t == GLPK.GLP_DB || t == GLPK.GLP_FX
            l[i] = GLPK.glp_get_row_lb(lp,reali)
            u[i] = GLPK.glp_get_row_ub(lp,reali)
        end
    end

    Arows = Int[]
    Acols = Int[]
    Avals = Float64[]
    sizehint!(Arows, nnz)
    sizehint!(Acols, nnz)
    sizehint!(Avals, nnz)
    for j in 1:ncol
        nnz1 = GLPK.glp_get_mat_col(lp, j, index1, coef1)
        for k = 1:nnz1
            if index1[k+1] == objrow
                continue
            end
            if objrow > 0 && index1[k+1] > objrow
                push!(Arows, index1[k+1]-1)
            else
                push!(Arows, index1[k+1])
            end
            push!(Acols, j)
            push!(Avals, coef1[k+1])
        end
    end

    A = sparse(Arows, Acols, Avals, nrow, ncol)

    return c, xlb, xub, l, u, A
end


text:
# testing output of this implementation with the one
# of Enzo de Sena: http://www.desena.org/sweep/
Fs = 4e4            # Sampling frequency
c  = 343    # create new acoustic env with default values
Nt = round(Int64,4E4/4)   # Number of time samples
xs = (2.,1.5,1.)          # Source position
xr = (1.,2.,2.)           # Receiver position
L  = (4.,4.,4.)           # Room dimensions
beta =  0.93.*(ones(6)...,)  # Reflection coefficient
Rd = 0.00          # random displacement

Tw = 40            # samples of Low pass filter
Fc = 0.9           # cut-off frequency

# generate IR with randomization
@time h, = rim(xs,xr,L,beta,Nt,Fs; Tw = Tw, Fc = Fc, Rd = Rd)

hm = readdlm("h.txt")
@test norm(hm-h)<1e-8
println("matlab reference test passed")


text:
export lorenz63, lorenz96

"""

    lorenz63(du, u, p, t)

Lorenz-63 dynamical model ``u = [x, y, z]`` and ``p = [\\sigma, \\rho, \\mu]``:
```math
\\frac{dx}{dt} = σ(y-x) \\\\
\\frac{dy}{dt} = x(ρ-z) - y \\\\
\\frac{dz}{dt} = xy - βz \\\\
```

- [Example Catalog](@ref)
- [Lorenz system on wikipedia](https://en.wikipedia.org/wiki/Lorenz_system)
"""
function lorenz63(du, u, p, t)

    du[1] = p[1] * (u[2] - u[1])
    du[2] = u[1] * (p[2] - u[3]) - u[2]
    du[3] = u[1] * u[2] - p[3] * u[3]

end

"""
    lorenz96(S, t, F, J)

Lorenz-96 dynamical model. For ``i=1,...,N``:

```math
\\frac{dx_i}{dt} = (x_{i+1}-x_{i-2})x_{i-1} - x_i + F
```

where it is assumed that ``x_{-1}=x_{N-1},x_0=x_N`` and ``x_{N+1}=x_1``.
Here ``x_i`` is the state of the system and ``F`` is a forcing constant.

- [Lorenz 96 model on wikipedia](https://en.wikipedia.org/wiki/Lorenz_96_model)
"""
function lorenz96(dx, x, p, t)
    F = p[1]
    N = Int64(p[2])
    # 3 edge cases
    dx[1] = (x[2] - x[N-1]) * x[N] - x[1] + F
    dx[2] = (x[3] - x[N]) * x[1] - x[2] + F
    dx[N] = (x[1] - x[N-2]) * x[N-1] - x[N] + F
    # then the general case
    for n = 3:(N-1)
        dx[n] = (x[n+1] - x[n-2]) * x[n-1] - x[n] + F
    end
end


text:
module MultiInterface

export Interface,
       Policy,
       @interface,
       @impl,
       enable_multiinterface_validation,
       get_policy,
       set_policy,
       with_policy,
       @unpack,
       list_impls,
       get_imp,
       get_method

using MacroTools
import Parameters: with_kw, @unpack

SE = Union{Symbol,Expr}

__policy = nothing
__num_imps = Dict()
__modules = Dict{Symbol,Module}()
# __interface = Dict{Symbol, Signature}()

enable_display() = global display = true
disable_display() = global display = false
display = false

gen_argname() = Symbol(replace(string(gensym("arg")), "#"=>""))

mi_types(e::Expr) = length(e.args) == 1 ? e.args[1] : e.args[2]
mi_types(::Symbol) = :Any

mi_argname(arg::Expr) = length(arg.args) == 1 ? gen_argname() : arg.args[1]
mi_argname(arg::Symbol) = arg

mi_getname(f::Symbol) = f
mi_getname(f::Expr) = f.args[2].value

function register_module(str, mod)
    global __modules
    push!(__modules, str=>mod)
end

function get_module(mod, str)
    global __modules
    get(__modules, str, mod)
end

# TODO do we want to allow the same interface name in different modules?
function register_interface(source, str, mod, args, retval, wparams)
    global __interface
    haskey(__interface, str) && @warn "reregistering interface for $str at $source" was___=get_interface(str) willbe=Signature(mod, args, retval, wparams)
    push!(__interface, str=>Signature(mod, args, retval, wparams))
end

function get_interface(str)
    get(__interface, str, nothing)
end

function getvars(t::Type)
    v = []

    getvars(t) = nothing
    getvars(t::DataType) = append!(v, collect(t.parameters))
    getvars(t::UnionAll) = begin
        push!(v, t.var)
        getvars(t.body)
    end

    unique(v)
end

function get_and_inc_imp_number(imp_type)
    global __num_imps
    if ~(imp_type in keys(__num_imps))
        __num_imps[imp_type] = 0
    end

    __num_imps[imp_type] += 1
    return __num_imps[imp_type]
end

function get_policy()
    global __policy
    return __policy
end

function set_policy(new_policy)
    global __policy
    old_policy = __policy
    __policy = new_policy
    return old_policy
end

abstract type Interface end

abstract type Policy end
get_imp(policy::Policy, args...) = nothing

OptionalPolicy = Union{Policy, Nothing}

function list_impls(imp_type, args_type=nothing)
    if isnothing(args_type)
        args_type = Tuple{<:Interface, Vararg{Any}}
    else
        args_type = Tuple{<:Interface, args_type.parameters...}
    end
    return [get_impl(m.sig) for m in methods(imp_type, args_type)]
end
function get_impl(t::UnionAll)
    return get_impl(t.body)
end
function get_impl(t::DataType)
    return t.parameters[2]
end

function strip_module(s)
    l = findlast('.', s)
    l === nothing || return s[(l+1):end]
    s
end

function firstcaps_to_lowerunder(s)
    s = strip_module(s)
    return lowercase(strip(replace(s,r"[A-Z]"=>s"_\0"),['_']))
end

function lowerunder_to_firstcaps(s)
    s = strip_module(s)
    return replace(titlecase(s),"_"=>"")
end

nothing2any(::Nothing) = :Any
nothing2any(x::Union{Symbol,Expr}) = x

"""
    macro interface(interface_exp)

Called
```
@interface a(x::I)::O where {I,O}
```

Generates the following code:
```
abstract type A <: Interface end

function (a(x::I; )::O) where {I,O}
    policy = get_policy()
    return a(policy, x::I)
end

function (a(policy::Policy, x::I; )::O) where {I,O}
    imp = get_imp(policy, A, x::I)
    return a(imp, x::I)
end
```
"""
macro interface(interface_exp)
    return esc(interface_macro(__module__, __source__, interface_exp))
end

function interface_macro(mod, source, interface_exp)
    f, args, R, W = interface_capture(interface_exp)
    args = add_missing_argnames(args)
    if display
        println("name: ", f, ": ", typeof(f))
        println("args: ", args, ": ", typeof(args))
        println("returntype: ", R, ": ", typeof(R))
        println("where: ", W, ": ", typeof(W))
    end

    W = isnothing(W) ? () : W

    interface_type_name = Symbol(lowerunder_to_firstcaps(String(f)))

    body = quote
        policy = get_policy()
        return $f(policy, $(get_arg_vars(args)...))
    end

    use_policy_body = quote
        imp = get_imp(policy, $interface_type_name, $(get_arg_vars(args)...))
        return $f(imp, $(get_arg_vars(args)...))
    end

    bare_func_sig_dict = Dict(:name => f,
                              :args => args,
                              :kwargs => Any[],
                              :body => body,
                              :rtype => R,
                              :whereparams => W,
                              )

    default_func_sig_dict = Dict(:name => f,
                                 :args => [:(policy::Policy); args],
                                 :kwargs => Any[],
                                 :body => use_policy_body,
                                 :rtype => R,
                                 :whereparams => W,
                                 )

    register_module(interface_type_name, mod)
    # WIP
    # register_interface(source, interface_type_name, mod, get_arg_types(args), nothing2any(R), W)

    result = quote
        abstract type $interface_type_name <: Interface end
        $(combinedef(bare_func_sig_dict))
        $(combinedef(default_func_sig_dict))
    end
    if display
        println(prettify(result))
    end

    return result
end

function get_arg_vars(args)
    tovar(s::Symbol) = s
    tovar(expr::Expr) = tovar(Val(expr.head), expr)
    tovar(::Val{:(::)}, expr) = expr.args[1]

    map(tovar, args)
end

function add_missing_argnames(args)
    map(args) do arg
        m = match(r"^::\s*(.*)", string(arg))
        if m !== nothing
            arg = copy(arg)
            insert!(arg.args, 1, gen_argname())
        end
        return arg
    end
end

function get_arg_type_sigs(args)
    map(args) do arg
        m = match(r".*?::\s*(.*)", string(arg))
        s = m === nothing ? :Type : Meta.parse(string("Type{<:", m.captures[1],"}"))
        return Expr(:(::), s)
    end
end

function get_arg_types(args)::Vector{Union{Symbol, Expr}}
    map(args) do arg
        m = match(r".*?::\s*(.*)", string(arg))
        return m === nothing ? :Any : Meta.parse(m.captures[1])
    end
end

function mi_join_args(func_sig_dicts)
    zippedargs = zip([[mi_types(a) for a in fsd[:args]] for fsd in func_sig_dicts]...)
    argnames = [mi_argname(a) for a in first(func_sig_dicts)[:args]]

    function mi_merge_into!(s::Expr, t)
        t in s.args && return s
        push!(s.args, t)
        return s
    end

    do_union(s,t) = :(Union{$s,$t})
    mi_merge(::Val, s::Expr, t::Symbol) = do_union(s,t)
    mi_merge(::Val, s::Expr, ::Val, t::Expr) = do_union(s,t)

    function mi_merge(::Val{:Union}, s::Expr, ::Val{:Union}, t::Expr)
        v = deepcopy(s)
        append!(v.args, t.args)
        unique!(v.args)
        return v
    end

    mi_merge(::Val{:Union}, s::Expr, t::Symbol) = mi_merge_into!(deepcopy(s), t)
    mi_merge(::Val{:Union}, s::Expr, ::Val, t::Expr) = mi_merge_into!(deepcopy(s), t)
    mi_merge(::Val, s::Expr, ::Val{:Union}, t::Expr) = mi_merge(t, s)

    mi_merge(s::Symbol, t::Symbol) = s === t ? s : do_union(s,t)
    mi_merge(s::Expr, t::Symbol) = mi_merge(Val(s.args[1]), s, t)
    mi_merge(s::Symbol, t::Expr) = mi_merge(t,s)

    mi_merge(s::Expr, t::Expr) = mi_merge(Val(s.args[1]), s, Val(t.args[1]), t)

    args = [reduce(mi_merge, arg) for arg in zippedargs]
    exprs = map(argnames, args) do n,arg
        :($(n)::$(arg))
    end
    collect(exprs)
end

parameter_names(e::Symbol) = e
parameter_names(e::Expr) = length(e.args) == 1 ? gen_argname() : e.args[1]

"""
    macro impl(expr)

Called
```
@impl begin
    struct MyA
        precision::Int64
    end
    function a(x::Int)
      c = 1
      return x + c + precision
    end
end
```

The `MyA` can be considered an identifier for this particular method, as opposed to different implementations of
this function.

This macro generates the following code:
```
begin
    Base.@__doc__ struct MyA <: A
            precision::Int64
            MyA(; precision = error("Field '" * "precision" * "' has no default, supply it with keyword.")) = MyA(precision)
            MyA(precision) = new(precision)
        end
    ()
    ()
    MyA(pp::MyA; kws...) = (Parameters).reconstruct(pp, kws)
    MyA(pp::MyA, di::(Parameters).AbstractDict) = (Parameters).reconstruct(pp, di)
    MyA(pp::MyA, di::Vararg{Tuple{Symbol, Any}}) = (Parameters).reconstruct(pp, di)
    nothing
    macro unpack_MyA(ex)
        esc((Parameters)._unpack(ex, Any[:precision]))
    end
    macro pack_MyA()
        esc((Parameters)._pack_new(MyA, Any[:precision]))
    end
    MyA
    function a(impl::MyA, x::Int; )
        @unpack (precision,) = impl
        c = 1
        return x + c + precision
    end
    function a(policy::Nothing, x::Int; )
        return a(MyA(), x::Int)
    end
    @generated function get_method(imp::Type{MyA})
            m = ([m for m = methods(a, Tuple{MyA, Vararg})])[1]
            return m
        end
    (s::MyA)(args...) = a(s, args...)
end
```
"""
macro impl(expr)
    return esc(implement_macro(__module__, __source__, expr))
end

function implement_macro(mod, source, implement_expr)
    stripped = MacroTools.prewalk(rmlines, implement_expr)
    T, fields, fname, func_sig_dicts = struct_capture(stripped)

    interface_name = Symbol(lowerunder_to_firstcaps(string(fname)))
    impl_type_name = isnothing(T) ?
                         Symbol(string(interface_name, get_and_inc_imp_number(interface_name))) :
                         T

    first_func_sig_dict = first(func_sig_dicts)
    specific_imp_func_decs = Dict[]
    for func_sig_dict in func_sig_dicts
        specific_imp_func_dec = deepcopy(func_sig_dict)
        specific_imp_func_dec[:args] = [:(impl::$impl_type_name); specific_imp_func_dec[:args]]
        if length(fields) > 0
            fieldnames = [namify(f) for f in fields]
            # Have to do some wonky thing here to get this to parse right
            unpack_expr = :(@unpack $(fieldnames[1:(end-1)]...), $(fieldnames[end]) = impl)
        else
            unpack_expr = :(begin end)
        end
        specific_imp_func_dec[:body] = quote
            try
                $unpack_expr
                $(specific_imp_func_dec[:body])
            catch
                src = $(string(source))
                @error "error at $(src)"
                rethrow()
            end
        end
        push!(specific_imp_func_decs, specific_imp_func_dec)
    end

    joined_args = mi_join_args(func_sig_dicts)
    joined_params = [parameter_names(e) for e in joined_args]
    default_imp_func_dec = deepcopy(first_func_sig_dict)
    default_imp_func_dec[:args] = [:(impl::Nothing); joined_args]
    default_imp_func_dec[:whereparams] = unique(reduce((x,fsd)->append!(x,get(fsd,:whereparams,[])), func_sig_dicts;init=[]))
    default_imp_func_dec[:body] = quote
        try
            return $(default_imp_func_dec[:name])($impl_type_name(), $(joined_params...))
        catch
            src = $(string(source))
            @error "error at $(src)"
            rethrow()
        end
    end
    # TODO should do a typejoin on all retvals
    default_imp_func_dec[:rtype] = :Any

    if ~(:rtype in keys(first_func_sig_dict))
        first_func_sig_dict[:rtype] = :Any
    end

    interface_module = get_module(mod, interface_name)

    if length(fields) >= 1
        struct_exp = with_kw(:(struct $impl_type_name <: $(interface_module).$(interface_name)
                                   $(fields...)
                               end),
                             @__MODULE__,
                             false)
    else
        struct_exp = quote
            struct $impl_type_name <: $(interface_module).$(interface_name)
                $(fields...)
            end
        end
    end

    wparams = reduce(append!, [get(fsd, :whereparams, []) for fsd in func_sig_dicts]; init=[])
    wparams = unique(wparams)
    argtypes = [[get(fsd,:rtype,:Any);get_arg_types(fsd[:args])] for fsd in func_sig_dicts]
    source_str = string(source)

    result = quote
        $struct_exp
        $([combinedef(sifd) for sifd in specific_imp_func_decs]...)
        $(combinedef(default_imp_func_dec))
        @generated function get_method(imp::Type{$impl_type_name})
            m = [m for m in methods($(first_func_sig_dict[:name]), Tuple{$impl_type_name, Vararg{Any}})][1]
            return m
        end
        (s::$impl_type_name)(args...) = $(first_func_sig_dict[:name])(s, args...)
    end
    if display
        println(prettify(result))
    end
    return result
end

function interface_capture(interface_exp)
    @capture(interface_exp, (f_(args__)::R_ where W__) |
                            (f_(args__) where W__) |
                            (f_(args__)::R_) |
                            (f_(args__)))
    return f, args, R, W
end

function struct_capture(expr)
    if expr.args[1].head === :struct
        @capture(expr.args[1], (struct T_ fields__ end) | (mutable struct T_ fields__ end))
        fs = [splitdef(a) for a in expr.args[2:end]]
    else
        T = nothing
        fields = []
        fs = [splitdef(a) for a in expr.args]
    end
    fname = first(fs)[:name]
    T,fields,fname,fs
end

###
### USING POLICIES
###
function with_policy(f, new_policy)
    last_policy = set_policy(new_policy)
    result = nothing
    try
        result = f()
    catch err
        rethrow()
    finally
        _ = set_policy(last_policy)
    end
    return result
end

end # module


text:
# Run package tests
println("Testing Silo.jl in Julia version ", VERSION)

using Base.Test
include(joinpath("..", "src", "Silo.jl"))
# using Silo
# run(`cd $(dirname(@__FILE__))/files && make`)

include("test1dwriteInt.jl")
include("test1dreadwrite.jl")

# Silo.DBInqFile(dbfile.file_name)


text:
using Languages

@assert all(articles(EnglishLanguage) .== ["a", "an", "the"])
@assert all(indefinite_articles(EnglishLanguage) .== ["a", "an"])
@assert all(definite_articles(EnglishLanguage) .== ["the"])
@assert length(prepositions(EnglishLanguage)) == 33
@assert length(pronouns(EnglishLanguage)) == 36
@assert length(stopwords(EnglishLanguage)) == 488


text:
#Sigmoid function in Julia
#(presumably) accepts Float & Int arguments
function sigmoid(x)
    return (MathConstants.e^x)/(1 + MathConstants.e^x)
end

#Generates all primes under the limit
#only takes positive integer values, greater than 2
function genprimes(limit)
    list = fill(true, limit)
    primes = []
    for i in 2:isqrt(limit)
        if list[i] == true
            j = i^2
            while j <= limit
                list[j] = false
                j += i
            end
        end
    end
    for i in 2:limit
        if list[i] == true
            push!(primes, i)
        end
    end
    return primes
end

max_stars_repo_name: ian-double-u/NNfromScratch

max_issues_repo_name: ian-double-u/NNfromScratch

max_forks_repo_name: ian-double-u/NNfromScratch


text:
using PowerDynamics: AbstractLine, construct_edge
using LinearAlgebra: I
import PowerDynamics: PiModel
using NetworkDynamics: ODEEdge

begin
    @__doc__ struct RLLine <: AbstractLine
            from
            to
            R
            L
            ω0
            RLLine(; from, to, R, L, ω0) = new(from, to, R, L, ω0)
        end
    function PowerDynamics.construct_edge(par::RLLine)
        from = par.from
        to = par.to
        R = par.R
        L = par.L
        ω0 = par.ω0
        @assert R > 0 "resistance (R) should be positive"
        @assert L > 0 "inductance (L) should be positive"
        @assert ω0 > 0 "rated frequency (ω0) should be positive"
        Z = [R -ω0*L; ω0*L R]
        function rhs!(de, e, v_s, v_d, p, t)
            # the different minus signs are due to the PowerDynamics sign convention for currents
            i_left_right = [e[1]; e[2]]
            i_right_left = [e[3]; e[4]]
            v_left = [v_s[1]; v_s[2]]
            v_right = [v_d[1]; v_d[2]]
            v_left_right = (- Z * i_left_right .- v_left .+ v_right) ./ L
            v_right_left = (- Z * i_right_left .+ v_right .- v_left) ./ L
            de .= [v_left_right; v_right_left]
        end
        #TODO check if there's a bug in the macro missing the empty space
        return ODEEdge(f! = rhs!, dim=4, mass_matrix=I, sym=Symbol[:i_r, :i_i, :ir_r, :ir_i])
    end
    symbolsof(::RLLine) = begin
            [:i_r, :i_i, :ir_r, :ir_i]
    end
    dimension(::RLLine) = begin
            4
    end
end

PiModel(l::RLLine) = PiModel(1/(l.R + im *l.ω0*l.L), 0, 0, 1, 1)

max_stars_repo_name: luap-pik/PowerDynIC

max_issues_repo_name: luap-pik/PowerDynIC

max_forks_repo_name: luap-pik/PowerDynIC


text:
export LowRankFactorizedNEP

"
    struct LowRankFactorizedNEP <: AbstractSPMF
    function LowRankFactorizedNEP(L::Vector,U::Vector,f::Vector)
    function LowRankFactorizedNEP(L::Vector,U::Vector,A::Vector, f::Vector)

Representation of a `NEP` which has low rank in the sense that it is an `SPMF`
where each of the terms are factorized: `A[i]=L[i]*U[i]'`. The factorization
is provided in the `L` and `U` vectors and the full matrix `A[i]` can be
either provided (or is otherwise implicitly computed).

# Example:
```
julia> L=randn(5,1); U=randn(5,1);
julia> f=S->exp(S)
julia> nep=LowRankFactorizedNEP([L],[U],[f]);
julia> X=randn(5,2);
julia> norm(compute_Mlincomb(nep,0.0,X)-L*U'*X*ones(2),1)
6.661338147750939e-16
```

"
struct LowRankFactorizedNEP{S<:AbstractMatrix{<:Number}} <: AbstractSPMF{AbstractMatrix}
    spmf::SPMF_NEP
    r::Int          # Sum of ranks of matrices
    L::Vector{S}    # Low rank L factors of matrices
    U::Vector{S}    # Low rank U factors of matrices
end

function LowRankFactorizedNEP(L::AbstractVector{S}, U::AbstractVector{S},
                              A::AbstractVector{S}, f) where {S<:AbstractMatrix}
    rank = mapreduce(u -> size(u, 2), +, U)
    spmf=SPMF_NEP(A, f, align_sparsity_patterns=true, check_consistency=false);
    return LowRankFactorizedNEP{S}(spmf, rank, L, U)
end
function LowRankFactorizedNEP(L::AbstractVector{S}, U::AbstractVector{S},
                              f) where {S<:AbstractMatrix}
    A = L .* adjoint.(U);
    return LowRankFactorizedNEP(L,U,A,f)
end

LowRankFactorizedNEP(::Type{T}, n) where T<:Number =
    LowRankFactorizedNEP(SPMF_NEP(n), 0, Vector{Matrix{T}}(), Vector{Matrix{T}}())

# forward function calls to SPMF
compute_Mder(nep::LowRankFactorizedNEP, λ::T, i::Int = 0) where T<:Number =
    compute_Mder(nep.spmf, λ, i)

compute_Mlincomb(nep::LowRankFactorizedNEP, λ::Number, V::AbstractVecOrMat, a::Vector = ones(eltype(V), size(V, 2)))  =
    compute_Mlincomb(nep.spmf, λ, V, a)

compute_MM(nep::LowRankFactorizedNEP, par...)  =
    compute_MM(nep.spmf, par...)

size(nep::LowRankFactorizedNEP) = size(nep.spmf)
size(nep::LowRankFactorizedNEP, dim) = size(nep.spmf, dim)

get_Av(nep::LowRankFactorizedNEP) = get_Av(nep.spmf)
get_fv(nep::LowRankFactorizedNEP) = get_fv(nep.spmf)


text:
ps = PathSet(; symlink=true)

@testset "$(typeof(ps.root))" begin
    testsets = [
        test_constructor,
        test_registration,
        test_show,
        test_parse,
        test_convert,
        test_components,
        test_parents,
        test_join,
        test_basename,
        test_filename,
        test_extensions,
        test_isempty,
        test_norm,
        test_real,
        test_relative,
        test_abs,
        test_isdir,
        test_isfile,
        test_stat,
        test_size,
        test_modified,
        test_created,
        test_cd,
        test_readpath,
        test_walkpath,
        test_read,
        test_write,
        test_mkdir,
        test_cp,
        test_mv,
        test_sync,
        test_symlink,
        test_touch,
        test_tmpname,
        test_tmpdir,
        test_mktmp,
        test_mktmpdir,
        test_download,
    ]

    if isa(ps.root, PosixPath)
        append!(
            testsets,
            [
                test_issocket,
                test_isfifo,
                test_ischardev,
                test_isblockdev,
                test_ismount,
                test_isexecutable,
                test_isreadable,
                test_iswritable,
                test_chown,
                test_chmod,
            ]
        )
    end

    # Run all of the automated tests
    test(ps, testsets)

    # Test the system path specific macros
    cd(abs(parent(Path(@__FILE__)))) do
        @test @__PATH__() == Path(@__DIR__)
        @test @__FILEPATH__() == Path(@__FILE__)
        @test FilePathsBase.@LOCAL("foo.txt") == join(@__PATH__, "foo.txt")
        @test FilePathsBase.@LOCAL("foo.txt") == joinpath(@__PATH__, "foo.txt")
    end
end


text:
"""
    SFMDistribution{T <: Mesh}

Representation of the distribution of a DiscretisedFluidQueue.

The field `coeffs` is a row-vector and indexing a distribution indexes this field. Similarly,
arithmetic operations on a distribution act on this vector.

# Arguments:
- `coeffs::Array{Float64, 2}`: A row vector of coefficients which encode the distribution of a
    DiscretisedFluidQueue and can be used to reconstruct a solution.
- `dq::DiscretisedFluidQueue{T}`:
"""
struct SFMDistribution{T<:Mesh} <: AbstractVector{Float64}
    coeffs::Vector{Float64}
    dq::DiscretisedFluidQueue{T}
end

SFMDistribution(coeffs::Vector{Float64},dq::DiscretisedFluidQueue{T}) where T =
    SFMDistribution{T}(coeffs,dq)

SFMDistribution(coeffs::Matrix{Float64},dq::DiscretisedFluidQueue) =
    (size(coeffs,1)==1) ? SFMDistribution(coeffs[:],dq) : throw(DomainError("coeffs must be a row (1xn array) or column vector"))

"""
    SFMDistribution(dq::DiscretisedFluidQueue{T})

A blank initialiser for a fluid queue distribution a distribution with `coeffs=zeros`.
"""
SFMDistribution(dq::DiscretisedFluidQueue{T}) where T =
    SFMDistribution{T}(zeros(n_bases_per_phase(dq)*n_phases(dq)+N₊(dq))+N₋(dq))

size(d::SFMDistribution) = size(d.coeffs)
size(d::SFMDistribution,dim::Int) = size(d.coeffs,dim)
getindex(d::SFMDistribution,i) = d.coeffs[i]
setindex!(d::SFMDistribution,x,i) = (d.coeffs[i]=x)

show(io::IO, mime::MIME"text/plain", d::SFMDistribution) = show(io, mime, d.coeffs)

fast_mul(u::SFMDistribution,B::AbstractMatrix{Float64}) = SFMDistribution(transpose(u.coeffs)*B,u.dq)
fast_mul(B::AbstractMatrix{Float64},u::SFMDistribution) =
    throw(DomainError("you can only premultiply a SFMDistribution"))
fast_mul(u::SFMDistribution,x::Real) = SFMDistribution(u.coeffs*x,u.dq)
fast_mul(x::Real,u::SFMDistribution) = fast_mul(u,x)

include("approximation.jl")
include("reconstruction.jl")


text:
# ===================================
# Approximations in the infinity norm
# ===================================

"""
    box_approximation(S::LazySet{N})::Union{Hyperrectangle{N}, EmptySet{N}}
        where {N<:Real}

Overapproximate a convex set by a tight hyperrectangle.

### Input

- `S` -- convex set

### Output

A tight hyperrectangle.
"""
function box_approximation(S::LazySet{N}
                          )::Union{Hyperrectangle{N}, EmptySet{N}} where {N<:Real}
    return overapproximate(S, Hyperrectangle)
end

"""
    interval_hull

Alias for `box_approximation`.
"""
interval_hull = box_approximation

"""
    box_approximation_symmetric(S::LazySet{N}
                               )::Union{Hyperrectangle{N}, EmptySet{N}}
                                where {N<:Real}

Overapproximate a convex set by a tight hyperrectangle centered in the origin.

### Input

- `S` -- convex set

### Output

A tight hyperrectangle centered in the origin.

### Algorithm

The center of the box is the origin, and the radius is obtained by computing the
maximum value of the support function evaluated at the canonical directions.
"""
function box_approximation_symmetric(S::LazySet{N};
                                    )::Union{Hyperrectangle{N},
                                             EmptySet{N}} where {N<:Real}
    (c, r) = box_approximation_helper(S)
    if r[1] < 0
        return EmptySet{N}()
    end
    return Hyperrectangle(zeros(N, length(c)), abs.(c) .+ r)
end

# special case: empty set
box_approximation_symmetric(∅::EmptySet) = ∅

"""
    symmetric_interval_hull

Alias for `box_approximation_symmetric`.
"""
symmetric_interval_hull = box_approximation_symmetric

"""
    box_approximation_helper(S::LazySet{N}) where {N<:Real}

Common code of `box_approximation` and `box_approximation_symmetric`.

### Input

- `S` -- convex set

### Output

A tuple containing the data that is needed to construct a tightly
overapproximating hyperrectangle.

- `c` -- center
- `r` -- radius

### Algorithm

The center of the hyperrectangle is obtained by averaging the support function
of the given convex set in the canonical directions.
The lengths of the sides can be recovered from the distance among support
functions in the same directions.
"""
@inline function box_approximation_helper(S::LazySet{N}) where {N<:Real}
    zero_N = zero(N)
    one_N = one(N)
    n = dim(S)
    c = Vector{N}(undef, n)
    r = Vector{N}(undef, n)
    d = zeros(N, n)
    @inbounds for i in 1:n
        d[i] = one_N
        htop = ρ(d, S)
        d[i] = -one_N
        hbottom = -ρ(d, S)
        d[i] = zero_N
        c[i] = (htop + hbottom) / 2
        r[i] = (htop - hbottom) / 2
        if r[i] < 0
            # contradicting bounds => set is empty
            # terminate with first radius entry being negative
            r[1] = r[i]
            break
        end
    end
    return c, r
end

"""
    ballinf_approximation(S::LazySet{N}
                         )::Union{BallInf{N}, EmptySet{N}} where {N<:Real}

Overapproximate a convex set by a tight ball in the infinity norm.

### Input

- `S` -- convex set

### Output

A tight ball in the infinity norm.
"""
function ballinf_approximation(S::LazySet{N}
                              )::Union{BallInf{N}, EmptySet{N}} where {N<:Real}
    return overapproximate(S, BallInf)
end

# special case: empty set
ballinf_approximation(∅::EmptySet) = ∅


text:
using ROS
ROS.@include geometry_msgs: TransformStamped, PoseStamped

function init()
    ROS.init("testTF")
    nh = ROS.NodeHandle()

    b = ROS.tf.TransformBroadcaster()
    sb = ROS.tf.StaticTransformBroadcaster()
    bf = ROS.tf.Buffer()
    tl = ROS.tf.TransformListener(bf)
    ts = ROS.geometry_msgs_TransformStamped()
    ts.header.stamp = ROS.now()
    ts.header.frame_id = "base_link"
    ts.child_frame_id = "gps_link"
    ts.transform.translation.x = 2.3
    ts.transform.translation.y = 1.6
    ts.transform.rotation.w = 1
    ts2 = ROS.geometry_msgs_TransformStamped()
    ts2.header.stamp = ROS.now()
    ts2.header.frame_id = "base_link"
    ts2.child_frame_id = "camera_link"
    ts2.transform.translation.x = 0.3
    ts2.transform.translation.y = 1
    ts2.transform.rotation.w = 1

    ROS.tf.sendTransform(b, ts)

    # Send multiple TFs at once
    tss = [ts,ts2]
    ROS.tf.sendTransform(sb, tss)

    if ROS.tf.canTransform(bf, "base_link", "camera_link", ROS.Time(0), ROS.Duration(1))
        stamped = ROS.tf.lookupTransform(bf, "base_link", "camera_link", ROS.Time(0), ROS.Duration(1))
        println(stamped.transform.translation.x)
        println(stamped.header.frame_id)

        p = ROS.geometry_msgs_PoseStamped()
        p.header.frame_id = "camera_link"
        p.pose.position.x = 1
        p.pose.position.y = 0.6
        p.pose.orientation.w = 1
        try
            transformed = ROS.tf.transform(bf, p, "base_link", ROS.Time(0), "gps_link", ROS.Duration(1))
            println(transformed.pose.position.x)
        catch
            println("TF Exception :(")
        end
    end
end

init()


text:
@testset "API" begin

    #########
    # SETUP #
    #########

    using Evolutionary: value!, f_calls, BoxConstraints, PenaltyConstraints,
                        EvolutionaryObjective

    # objective function
    dimension = 5
    individual = ones(dimension)
    minval = 1.0
    func = x->sum(x)+1
    objfun = EvolutionaryObjective(func, individual)
    st = TestOptimizerState(individual, minval);
    mthd = TestOptimizer()
    cnstr= Evolutionary.NoConstraints()

    # options
    opts = Evolutionary.Options(store_trace=true; default_options(mthd)...)
    show(IOBuffer(), opts)

    # population
    ppl_size = 7
    ppl = [rand(dimension) for i in 1:ppl_size]

    #########
    # TRACE #
    #########
    @test value(st) == minval
    val = value!(objfun, rand(dimension))

    tr = Evolutionary.OptimizationTrace{typeof(value(objfun)), typeof(mthd)}()
    @test !Evolutionary.trace!(tr, 1, objfun, st, ppl, mthd, opts)
    @test tr[1].iteration == 1
    @test tr[1].value == value(st)
    @test tr[1].metadata["time"] <= time()
    @test !update_state!(objfun, cnstr, st, ppl, mthd, opts, 10)
    @test !Evolutionary.trace!(tr, 2, objfun, st, ppl, mthd, opts)
    @test tr[2].iteration == 2
    @test tr[2].value >= 1
    @test tr[2].metadata["time"] >= tr[1].metadata["time"]
    show(IOBuffer(), tr)

    # callback on state
    opts = Evolutionary.Options(store_trace=false, callback=(s->s.iteration == 2))
    tr = Evolutionary.OptimizationTrace{typeof(value(objfun)), typeof(mthd)}()
    @test Evolutionary.trace!(tr, 1, objfun, st, ppl, mthd, opts) == false
    @test Evolutionary.trace!(tr, 2, objfun, st, ppl, mthd, opts) == true

    # callback on trace
    opts = Evolutionary.Options(store_trace=true, callback=(tr->length(tr)>1))
    tr = Evolutionary.OptimizationTrace{typeof(value(objfun)), typeof(mthd)}()
    @test Evolutionary.trace!(tr, 1, objfun, st, ppl, mthd, opts) == false
    @test Evolutionary.trace!(tr, 2, objfun, st, ppl, mthd, opts) == true

    # custom trace
    function Evolutionary.trace!(record::Dict{String,Any}, objfun, state::TestOptimizerState, population, method, options)
        record["ppl_size"] = length(population)
    end
    tr = Evolutionary.OptimizationTrace{typeof(value(objfun)), typeof(mthd)}()
    @test !Evolutionary.trace!(tr, 1, objfun, st, ppl, mthd, opts)
    @test tr[1].iteration == 1
    @test tr[1].metadata["ppl_size"] == length(ppl)

    ##########
    # RESULT #
    ##########
    opts = Evolutionary.Options(store_trace=true)
    res = Evolutionary.EvolutionaryOptimizationResults(
        mthd, individual, value(st),
        opts.iterations, opts.show_trace, opts.store_trace,
        Evolutionary.metrics(mthd), tr, f_calls(objfun),
        1.0, 1.0, opts.show_trace,
    );
    @test summary(res) == summary(mthd)
    @test Evolutionary.minimizer(res) == individual
    @test minimum(res) == value(st)
    @test Evolutionary.iterations(res) == opts.iterations
    @test Evolutionary.iteration_limit_reached(res) == opts.show_trace
    @test Evolutionary.converged(res) == opts.store_trace
    @test f_calls(res) == f_calls(objfun)
    @test Evolutionary.time_run(res) == 1.0
    @test Evolutionary.time_limit(res) == 1.0
    @test Evolutionary.is_moo(res) == opts.show_trace
    @test length(Evolutionary.trace(res)) >= 1
    @test !Evolutionary.converged(res.metrics[1])
    show(IOBuffer(), res)

    ##############
    # POPULATION #
    ##############
    pop = Evolutionary.initial_population(mthd, BitVector(ones(dimension)))
    @test length(pop) == population_size(mthd)
    @test size(first(pop)) == (dimension,)

    @test_throws AssertionError Evolutionary.initial_population(mthd, fill(BitVector(ones(dimension)),4))
    pop = Evolutionary.initial_population(mthd, fill(BitVector(ones(dimension)),5))
    @test length(pop) == population_size(mthd)
    @test size(first(pop)) == (dimension,)

    pop = Evolutionary.initial_population(mthd, BitMatrix(ones(dimension,6)))
    @test length(pop) == population_size(mthd)
    @test size(first(pop)) == (dimension,6)

    pop = Evolutionary.initial_population(mthd, (()->rand(Bool,dimension)))
    @test length(pop) == population_size(mthd)
    @test size(first(pop)) == (dimension,)

    lb = [0,1,2,-1]
    ub = [0,3,2,1]
    cb = Evolutionary.ConstraintBounds(lb,ub,[],[])
    pop = Evolutionary.initial_population(mthd, cb)
    @test length(pop) == population_size(mthd)
    @test map(i->i[1],pop) == fill(ub[1], population_size(mthd))
    @test all(map(i->lb[2] <= i[2] <= ub[2],pop))
    @test map(i->i[3],pop) == fill(ub[3], population_size(mthd))
    @test all(map(i->lb[4] <= i[4] <= ub[4],pop))

    ###########
    # OPTIONS #
    ###########
    opts = Evolutionary.Options(;store_trace=true, iterations=7)
    @test opts.store_trace
    @test opts.iterations == 7

    ###############
    # CONVERGENCE #
    ###############
    cm = Evolutionary.AbsDiff()
    @test !Evolutionary.converged(cm)
    st.fitness = 1.0
    @test !Evolutionary.assess!(cm, st)
    @test !Evolutionary.converged(cm)
    @test Evolutionary.assess!(cm, st)
    @test Evolutionary.converged(cm)

    ###############
    # CONSTRAINTS #
    ###############
    c = Evolutionary.NoConstraints()
    @test isfeasible(c, 1)
    @test isa(c.bounds, Evolutionary.ConstraintBounds)
    @test value(c, 1) === nothing
    @test apply!(c, 1) == 1
    @test penalty(c, 1) == 0

    c = BoxConstraints(lb, ub)
    @test isfeasible(c, [0, 1, 2, 0])
    @test value(c, [1,4,2,-2]) === nothing
    @test apply!(c, [1,4,2,-2]) == [0,3,2,-1]
    @test penalty(c, [1,4,2,-2]) == 0

    con_c!(x) = [sum(x)]
    cb = Evolutionary.ConstraintBounds(fill(0, 3), fill(1, 3), [1], [1])
    c = PenaltyConstraints(1, cb, con_c!)
    objfun = EvolutionaryObjective(sum, zeros(3))
    x, y = [0, 1, 0], [0, -1, 3]
    @test isfeasible(c, x)    # feasible
    @test !isfeasible(c, y) # not feasible
    @test apply!(c, x) == x
    @test apply!(c, y) == y
    @test value(c, x) == [sum(x)]
    @test value(c, y) == [sum(y)]
    @test penalty(c, x) == 0.0 # c penalty
    @test penalty(c, y) == 1+2^2+1.0 # c penalty
    @test penalty!([1, 2], c, [x,y]) == [1, 2 + penalty(c, y)]

    ############
    # OPTIMIZE #
    ############
    res =  Evolutionary.optimize(sum, (()->BitVector(rand(Bool,dimension))), mthd, opts)
    @test Evolutionary.minimum(res) <= dimension
    @test length(Evolutionary.minimizer(res)) == dimension
    @test Evolutionary.iterations(res) == 7
    @test !Evolutionary.converged(res)
    @test length(Evolutionary.trace(res)) > 1
    @test length(res.metrics) == 1
    @test res.metrics[1].tol == 1e-10

    o = Evolutionary.Options(abstol=1e-3)
    res =  Evolutionary.optimize(sum, (()->BitVector(rand(Bool,dimension))), mthd, o)
    @test length(res.metrics) == 1
    @test res.metrics[1].tol == 1e-3

    o = Evolutionary.Options(abstol=1e-3, reltol=1e-4)
    res =  Evolutionary.optimize(sum, (()->BitVector(rand(Bool,dimension))), mthd, o)
    @test length(res.metrics) == 2
    @test res.metrics[1].tol == 1e-3
    @test res.metrics[2].tol == 1e-4

end


text:
module LintExample
using Lint
export @fancyfuncgen

macro fancyfuncgen(sym)
    s = string(sym)
    code = Expr(:function,
        :($sym()),
        quote
            println($s)
        end)
    esc(code)
end

function lint_helper(ex::Expr, ctx::LintContext)
    if ex.head == :macrocall
        if ex.args[1] == Symbol("@fancyfuncgen")
            if typeof(ex.args[2]) != Symbol
                Lint.msg(ctx, 2, "@fancyfuncgen must use a symbol argument")
            else
                Lint.addconst!(ctx.callstack[end], ex.args[2], Function, Lint.location(ctx))
            end
            return true
        end
    end
    return false
end

end


text:
# find the interval for integral for a pdf

function zero_finding(f, kmax, tol, b0)
    # step 0:
    # test if b0 is a good guess
    while f(b0) < tol
        b0 = b0/2
    end

    b = b0
    k = 0
    iter = 0

    intvl = [b, 2*b]

    # step 1:
        # find an interval [a, b] such that f(a) > tol, f(b) < tol
    while (f(intvl[1]) > tol && f(intvl[2]) < tol) == false
        intvl[1] = intvl[2]
        intvl[2] *= 2
        # print("Current interval: [$(intvl[1]), $(intvl[2])] \n")
        # print("Function values: [$(f(intvl[1])), $(f(intvl[2]))] \n")
        iter += 1
    end

    # step 2
        # shrink such interval using bisection
    while k < kmax
        midpt = (intvl[1] + intvl[2])/2
        midvalue = f(midpt)
        if midvalue < tol
            intvl[2] = midpt
            # print("Current interval: [$(intvl[1]), $(intvl[2])] \n")
            # print("Function values: [$(f(intvl[1])), $(f(intvl[2]))] \n")
        else
            intvl[1] = midpt
            # print("Current interval: [$(intvl[1]), $(intvl[2])] \n")
            # print("Function values: [$(f(intvl[1])), $(f(intvl[2]))] \n")
        end
        k += 1
    end
    ttl_iter = iter + k
    return intvl[2], ttl_iter
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
###############################################################################
#
#   Abstract types
#
###############################################################################

# broad mathematical domains
# these contain the type classes of parent objects

abstract type Set end

abstract type Group <: Set end

abstract type Ring <: Group end

abstract type Field <: Ring end

# elements of mathematical domains

abstract type SetElem end

abstract type GroupElem <: SetElem end

abstract type RingElem <: GroupElem end

abstract type FieldElem <: RingElem end

# parameterized domains

abstract type Module{T <: RingElem} <: Group end

# elements of parameterised domains

abstract type ModuleElem{T <: RingElem} <: GroupElem end

# rings, fields etc, parameterised by an element type
# these are the type classes of different kinds of
# mathematical rings/fields/etc, which have a base ring,
# and for which a generic implementation is possible
# over that base ring

abstract type PolyRing{T} <: Ring end

abstract type MPolyRing{T} <: Ring end

abstract type SeriesRing{T} <: Ring end

abstract type ResRing{T} <: Ring end

abstract type FracField{T} <: Field end

# not always really mathematical rings
# later we'll distinguish matrix algebras
# from the generic case
abstract type MatSpace{T} <: Ring end

# mathematical objects parameterised by an element type
# these are the type classes of mathematical objects
# that have some kind of base ring, and a generic
# implementation is meaningful over that base ring

abstract type PolyElem{T} <: RingElem end

abstract type MPolyElem{T} <: RingElem end

abstract type ResElem{T} <: RingElem end

abstract type FracElem{T} <: FieldElem end

abstract type SeriesElem{T} <: RingElem end

abstract type RelSeriesElem{T} <: SeriesElem{T} end

abstract type AbsSeriesElem{T} <: SeriesElem{T} end

   # not always mathematical ring elements
   # later we'll maybe distinguish MatAlgebraElem, MatModuleElem
abstract type MatElem{T} <: RingElem end

# additional abstract types for parents, added ad hoc to form
# collections of types as needed by applications

abstract type FinField <: Field end    # for fq, fq_nmod, etc

# additional abstract types for elements, added ad hoc to form
# collections of types as needed by applications

abstract type FinFieldElem <: FieldElem end # for fq, fq_nmod, etc

################################################################################
#
#   Promotion system
#
# The promote_rule functions are not extending Base.promote_rule. The Nemo
# promotion system is orthogonal to the built-in julia promotion system. The
# julia system assumes that whenever you have a method signature of the form
# Base.promote_rule(::Type{T}, ::Type{S}) = R, then there is also a
# corresponding Base.convert(::Type{R}, ::T) and similar for S. Since we
# cannot use the julia convert system (we need an instance of the type and not
# the type), we cannot use the julia promotion system.
#
# The Nemo promotion system is used to define catch all functions for
# arithmetic between arbitrary ring elements.
#
################################################################################

promote_rule(T, U) = Union{}


text:
include("../imports.jl")

@testset ExtendedTestSet "ImageSegmentation" begin
    @testset ExtendedTestSet "2D" begin
        method = ImageSegmentation((16, 16), 1:4)
        testencoding(method.encodings, method.blocks)
        DLPipelines.checkmethod_core(method)
        @test_nowarn methodlossfn(method)
        @test_nowarn methodmodel(method, Models.xresnet18())
        @testset "Show backends" begin
            @testset "ShowText" begin
                FastAI.test_method_show(method, ShowText(Base.DevNull()))
            end
            @testset "ShowMakie" begin
                FastAI.test_method_show(method, ShowMakie())
            end
        end
    end
    @testset ExtendedTestSet "3D" begin
        method = BlockMethod(
            (Image{3}(), Mask{3}(1:4)),
            (
                ProjectiveTransforms((16, 16, 16), inferencefactor=8),
                ImagePreprocessing(),
                FastAI.OneHot()
            )
        )
        testencoding(method.encodings, method.blocks)
        DLPipelines.checkmethod_core(method)
        @test_nowarn methodlossfn(method)
    end

end


text:
## ----------------------------------------------------------------------------
# utils
const _CLEAR_BIBTEX_ENTRY_STRIP_SET = ['{', '}', '\n', ' ', ',', ';']
const _alphanum_str = "a-zA-Z0-9"
const _punct_str = "\\,\\-\\;\\:\\."
const _reg1 = Regex("[$_punct_str]{2}")
const _reg2 = Regex("\\s{2}")
const _reg3 = Regex("[^$(_alphanum_str)$(_punct_str)/\\s]")

function _clear_bibtex_entry(str::AbstractString)

    str = strip(str, _CLEAR_BIBTEX_ENTRY_STRIP_SET)
    str = replace(str, "\n" => "")

    _till_unchanged() do
        str = replace(str, _reg1 => ",")
        str = replace(str, _reg2 => " ")
        for c in ",;:."
            str = replace(str, " $c" => c)
        end
        str = replace(str, _reg3 => "")
        return str
    end
    return str
end

## ----------------------------------------------------------------------------
function _till_unchanged(f!::Function; n = 100)
    ret = f!()
    ret0 = ret
    for i in 1:n
        ret = f!()
        (ret0 == ret) && return ret
    end
    return ret
end

## ----------------------------------------------------------------------------
function _push_csv!(col, str::AbstractString)
    vals = split(str, ",", keepempty = false)
    for vali in vals
        push!(col, strip(vali))
    end
    return col
end

function _push_csv!(col, str::AbstractString, strs::AbstractString...)
    _push_csv!(col, str)
    for stri in strs
        _push_csv!(col, stri)
    end
    return col
end

## ----------------------------------------------------------------------------
function _doi_to_url(doi::AbstractString)
    if startswith(doi, "http")
        return doi
    elseif startswith(doi, "doi.org")
        return string("https://", doi)
    elseif startswith(doi, r"\d{1,10}\.\d{3,20}")
        return string("https://doi.org/", doi)
    end
    return ""
end

function camel_alphanumeric(str)
    strs = split(str, r"[^a-zA-Z0-9]"; keepempty = false)
    join(uppercasefirst.(strs), " ")
end

function _if_nonempty(f::Function, str, dflt = "")
    isempty(str) ? dflt : f(str)
end
_if_nonempty(str, val, dflt = "") = _if_nonempty((x) -> val, str, dflt)


text:
using PyPlot
using PyCall
using BilliardModels

VERSION < v"0.4-" && using Docile

@pyimport matplotlib.patches as patches

#import PyPlot.draw
export bdraw

@doc """The various `bdraw` functions draw parts of the billiard(hence the `b`
at the start of the name, to avoid using `draw` which conflicts with PyPlot.)""" ->
function bdraw(d::Disc, subplot, offset)  # default is the current axis

    circ = patches.Circle((d.centre.x + offset.x,
                           d.centre.y + offset.y),
                          d.radius, alpha=0.5)
    subplot[:add_patch](circ)
end

function bdraw(p::AbstractPlane, subplot)  # don't draw anything for plane? This is too naive
end

function bdraw(table::BilliardTable, subplot, offset::Vector2D{Int}=Vector2D(0,0))

    # Calculate the minimum and maximum of the points used to define the planes:
    xmax, xmin = -Inf, Inf
    ymax, ymin = -Inf, Inf

    for obstacle in table.obstacles
        if isa(obstacle, AbstractPlane)
            xmin = min(xmin, obstacle.point_on_plane.x)
            xmax = max(xmax, obstacle.point_on_plane.x)
            ymin = min(ymin, obstacle.point_on_plane.y)
            ymax = max(ymax, obstacle.point_on_plane.y)
        else
            bdraw(obstacle, subplot, offset)  # draw the non-plane objects
        end
    end

    # Draw a square bounding box according to the points found:
    subplot[:plot]( [xmin, xmax, xmax, xmin, xmin],
                    [ymin, ymin, ymax, ymax, ymin],"-k",lw=2)
end

function bdraw(xs::Array, subplot)

    x = [pt.x for pt in xs]
    y = [pt.y for pt in xs]
    subplot[:plot](x, y, "-", alpha=0.5)
    subplot[:axis]("image")

end

function bdraw(p::AbstractParticle, subplot, draw_velocity=true)
    pos = p.x
    vel = p.v

    subplot[:plot]([pos.x], [pos.y], "o")

    dt = 0.05

    if draw_velocity
        subplot[:arrow](pos.x, pos.y, dt*vel.x, dt*vel.y, head_width=0.05, head_length=0.05)
    end

end

max_stars_repo_name: dpsanders/BilliardModel

max_issues_repo_name: dpsanders/BilliardModel

max_forks_repo_name: dpsanders/BilliardModel


text:
export eval_field
"Evaluate the field at a point x"
eval_field(A::BEMOperator, density::Expansion, x, quad = A.quad) =
    eval_field(integraloperator(A), density, x, quad)

function eval_field(intop::BoundaryIntegralOperator, density::Expansion, x, quad)
    integrand = t -> eval_kernel(intop, x, t) * unsafe_weightfun(measure(intop), t) * density(t)
    DomainIntegrals.integral(quad, integrand, support(density))
end

eval_field(intop::BoundaryIntegralOperator, density::Expansion, x, quad::QuadQBF) =
    _eval_field(intop, dictionary(density), coefficients(density), x, quad)

function _eval_field(intop, dict, coef, x, quad::QuadQBF)
    z = zero(eltype(coef))
    for i in 1:length(coef)
        z = z + coef[i] * collocation_BEM_entry(x, (dict, i), intop, quad)
    end
    z
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 3bd0847c-f2b4-11ea-35d8-c7657a170cf9
using DrWatson

# ╔═╡ 3bd0c52c-f2b4-11ea-09e6-05dbd556433f
begin
	@quickactivate "StatisticalRethinkingTuring"
	using StatisticalRethinking
	using PlutoUI
end

# ╔═╡ 28342aca-f2b1-11ea-342d-95590e306ff4
md"## Clip-04-01-06s.jl"

# ╔═╡ 3bd18872-f2b4-11ea-2271-470ff5d51737
md"### snippet 4.1"

# ╔═╡ 3be20a12-f2b4-11ea-2179-2995ca45302f
md"###### No attempt has been made to condense this too fewer lines of code."

# ╔═╡ 4114e7b0-f85f-11ea-248a-e1c6723d2f1a
@bind noofwalks Slider(5:200, default=9)

# ╔═╡ 3be2d1ea-f2b4-11ea-09ce-1db161f69c7e
begin
	noofsteps = 20;
	pos = Array{Float64, 2}(rand(Uniform(-1, 1), noofsteps, noofwalks));
	pos[1, :] = zeros(noofwalks);
	csum = cumsum(pos, dims=1);
	mx = minimum(csum) * 0.9
end;

# ╔═╡ 3bf08696-f2b4-11ea-0fae-815209e20f46
md"###### Plot and annotate the random walks."

# ╔═╡ 3bfe421a-f2b4-11ea-12c7-f7f775c210a9
md"###### Generate 3 plots of densities at 3 different step numbers (4, 8 and 16)."

# ╔═╡ 3bff112a-f2b4-11ea-338a-791bc65b719f
begin
	f = Plots.font("DejaVu Sans", 6)
	p1 = plot(csum, leg=false, title="Random walks ($(noofwalks))")
	plot!(p1, csum[:, Int(floor(noofwalks/2))], leg=false, title="Random walks ($(noofwalks))", 				color=:black)
	plot!(p1, [5], seriestype="vline")
	annotate!(5, mx, text("step 4", f, :left))
	plot!(p1, [9], seriestype="vline")
	annotate!(9, mx, text("step 8", f, :left))
	plot!(p1, [17], seriestype="vline")
	annotate!(17, mx, text("step 16", f, :left))

	p2 = Vector{Plots.Plot{Plots.GRBackend}}(undef, 3);
	plt = 1
	for step in [4, 8, 16]
		indx = step + 1 								# We added the first line of zeros
		global plt
	  	fitl = fit_mle(Normal, csum[indx, :])
	  	lx = (fitl.μ-4*fitl.σ):0.01:(fitl.μ+4*fitl.σ)
	  	p2[plt] = density(csum[indx, :], legend=false, title="$(step) steps")
	 	plot!( p2[plt], lx, pdf.(Normal( fitl.μ , fitl.σ ) , lx ), fill=(0, .5,:orange))
	  	plt += 1
	end
	p3 = plot(p2..., layout=(1, 3))
	plot(p1, p3, layout=(2,1))
end

# ╔═╡ 3c10cb0e-f2b4-11ea-210f-1f6986fbb7d8
md"## snippet 4.2"

# ╔═╡ 3c1b683e-f2b4-11ea-206b-b9178296ec19
prod(1 .+ rand(Uniform(0, 0.1), 12))

# ╔═╡ 3c248840-f2b4-11ea-1397-9f2c313a2676
md"## snippet 4.3"

# ╔═╡ 3c26427c-f2b4-11ea-23f9-6303c62de7f6
begin
	growth = [prod(1 .+ rand(Uniform(0, 0.1), 12)) for i in 1:10000];
	fit2 = fit_mle(Normal, growth)
	plot(Normal(fit2.μ , fit2.σ ), fill=(0, .5,:orange), lab="Normal distribution")
	density!(growth, lab="'sample' distribution")
end

# ╔═╡ 3c37aa58-f2b4-11ea-2e87-a11117582011
md"## snippet 4.4"

# ╔═╡ 3c390506-f2b4-11ea-3232-296f6952aaa7
begin
	big = [prod(1 .+ rand(Uniform(0, 0.5), 12)) for i in 1:10000];
	small = [prod(1 .+ rand(Uniform(0, 0.01), 12)) for i in 1:10000];
	fitb = fit_mle(Normal, big)
	fits = fit_mle(Normal, small)
	p5 = plot(Normal(fitb.μ , fitb.σ ), lab="Big normal distribution", fill=(0, .5,:orange))
	p4 = plot(Normal(fits.μ , fits.σ ), lab="Small normal distribution", fill=(0, .5,:orange))
	density!(p5, big, lab="'big' distribution")
	density!(p4, small, lab="'small' distribution")
	plot(p5, p4, layout=(1, 2))
end

# ╔═╡ 3c4245ee-f2b4-11ea-2add-750014544e83
md"## snippet 4.5"

# ╔═╡ 3c4be516-f2b4-11ea-086b-f1461fa2e55e
begin
	log_big = [log(prod(1 .+ rand(Uniform(0, 0.5), 12))) for i in 1:10000];
	fit3 = fit_mle(Normal, log_big)
	plot(Normal(fit3.μ , fit3.σ ), fill=(0, .5,:orange), lab="Normal distribution")
	density!(log_big, lab="'sample' distribution")
end

# ╔═╡ 3c567b9a-f2b4-11ea-2798-b5b3bffc3f84
md"## snippet 4.6"

# ╔═╡ 3c60487a-f2b4-11ea-32f7-abda1ad883cf
md"###### Grid of 1001 steps."

# ╔═╡ 3c6a0c78-f2b4-11ea-12a9-a506e7a7cacb
p_grid = range(0, step=0.001, stop=1);

# ╔═╡ 3c7435ce-f2b4-11ea-36a4-e72a5024cfa3
md"###### All priors = 1.0."

# ╔═╡ 3c7f4458-f2b4-11ea-3006-fb69892e25da
prior = ones(length(p_grid));

# ╔═╡ 3c89c428-f2b4-11ea-3216-7dab773dfe79
md"###### Binomial pdf."

# ╔═╡ 3c9436d8-f2b4-11ea-3b46-fbfefa79d535
likelihood = [pdf(Binomial(9, p), 6) for p in p_grid];

# ╔═╡ 3ca753ee-f2b4-11ea-21f5-93f4fee09437
md"###### A Uniform prior has been used, unstandardized posterior is equal to likelihood."

# ╔═╡ 3cab5208-f2b4-11ea-0cff-ddab2410c607
posterior = likelihood .* prior;

# ╔═╡ 3cb6dac6-f2b4-11ea-0c94-a109c74c2d22
md"###### Scale posterior such that they become probabilities."

# ╔═╡ 3cc23a06-f2b4-11ea-3dd3-c5ca6ebe71c3
posterior2 = posterior / sum(posterior);

# ╔═╡ 3cd03f70-f2b4-11ea-1877-11b0102993c8
md"###### Sample using the computed posterior values as weights. In this example we keep the number of samples equal to the length of p_grid, but that is not required."

# ╔═╡ 3cdd29a6-f2b4-11ea-071d-fb6328c01d5f
begin
	samples = sample(p_grid, Weights(posterior), length(p_grid));
	p = Vector{Plots.Plot{Plots.GRBackend}}(undef, 2)
	p[1] = scatter(1:length(p_grid), samples, markersize = 2, ylim=(0.0, 1.3), lab="Draws")
end

# ╔═╡ 3ce9eac4-f2b4-11ea-0bf6-9f130255ec5a
md"###### Analytical calculation."

# ╔═╡ 3cf665a6-f2b4-11ea-29ad-19ecaf1f423a
begin
	w = 6
	n = 9
	x = 0:0.01:1
	p[2] = density(samples, ylim=(0.0, 5.0), lab="Sample density")
	p[2] = plot!( x, pdf.(Beta( w+1 , n-w+1 ) , x ), lab="Conjugate solution")
end

# ╔═╡ 3d03bb5c-f2b4-11ea-315e-6fac666895eb
md"###### Quadratic approximation."

# ╔═╡ 3d11472c-f2b4-11ea-3bf6-9d517b13b291
begin
	plot!( p[2], x, pdf.(Normal( 0.67 , 0.16 ) , x ), lab="Normal approximation", fill=(0, .5,:orange))
	plot(p..., layout=(1, 2))
end

# ╔═╡ 3d1db2f0-f2b4-11ea-3bde-29b8f0be3087
md"## End of clip-04-01-06s.jl"

# ╔═╡ Cell order:
# ╟─28342aca-f2b1-11ea-342d-95590e306ff4
# ╠═3bd0847c-f2b4-11ea-35d8-c7657a170cf9
# ╠═3bd0c52c-f2b4-11ea-09e6-05dbd556433f
# ╟─3bd18872-f2b4-11ea-2271-470ff5d51737
# ╟─3be20a12-f2b4-11ea-2179-2995ca45302f
# ╠═4114e7b0-f85f-11ea-248a-e1c6723d2f1a
# ╠═3be2d1ea-f2b4-11ea-09ce-1db161f69c7e
# ╟─3bf08696-f2b4-11ea-0fae-815209e20f46
# ╟─3bfe421a-f2b4-11ea-12c7-f7f775c210a9
# ╠═3bff112a-f2b4-11ea-338a-791bc65b719f
# ╟─3c10cb0e-f2b4-11ea-210f-1f6986fbb7d8
# ╠═3c1b683e-f2b4-11ea-206b-b9178296ec19
# ╠═3c248840-f2b4-11ea-1397-9f2c313a2676
# ╠═3c26427c-f2b4-11ea-23f9-6303c62de7f6
# ╠═3c37aa58-f2b4-11ea-2e87-a11117582011
# ╠═3c390506-f2b4-11ea-3232-296f6952aaa7
# ╠═3c4245ee-f2b4-11ea-2add-750014544e83
# ╠═3c4be516-f2b4-11ea-086b-f1461fa2e55e
# ╟─3c567b9a-f2b4-11ea-2798-b5b3bffc3f84
# ╟─3c60487a-f2b4-11ea-32f7-abda1ad883cf
# ╠═3c6a0c78-f2b4-11ea-12a9-a506e7a7cacb
# ╟─3c7435ce-f2b4-11ea-36a4-e72a5024cfa3
# ╠═3c7f4458-f2b4-11ea-3006-fb69892e25da
# ╟─3c89c428-f2b4-11ea-3216-7dab773dfe79
# ╠═3c9436d8-f2b4-11ea-3b46-fbfefa79d535
# ╟─3ca753ee-f2b4-11ea-21f5-93f4fee09437
# ╠═3cab5208-f2b4-11ea-0cff-ddab2410c607
# ╠═3cb6dac6-f2b4-11ea-0c94-a109c74c2d22
# ╠═3cc23a06-f2b4-11ea-3dd3-c5ca6ebe71c3
# ╠═3cd03f70-f2b4-11ea-1877-11b0102993c8
# ╠═3cdd29a6-f2b4-11ea-071d-fb6328c01d5f
# ╟─3ce9eac4-f2b4-11ea-0bf6-9f130255ec5a
# ╠═3cf665a6-f2b4-11ea-29ad-19ecaf1f423a
# ╟─3d03bb5c-f2b4-11ea-315e-6fac666895eb
# ╠═3d11472c-f2b4-11ea-3bf6-9d517b13b291
# ╟─3d1db2f0-f2b4-11ea-3bde-29b8f0be3087


text:
using GCMAES, MPI

ENV["SLURM_ARRAY_TASK_ID"] = 1
ENV["SLURM_ARRAY_TASK_COUNT"] = 2

@mpirun GCMAES.make_virtual_jobarray(2)

rank = GCMAES.myrank(MPI.COMM_WORLD)
@show rank, GCMAES.myrank(), GCMAES.partjob(1:8), GCMAES.part(GCMAES.partjob(1:8))


text:
# Trotter factorized solution for systems with few modes.

"""
Trotter factorized solution for a small system.
"""
struct Trotter <: Solution
    "Partition function."
    Z::Float64
    "Energy."
    E::Float64
    "Heat capacity."
    Cv::Float64
end

"""
    Trotter(sys::System, beta::Float64, basis_size::Int, P::Int)

Calculate the solution for `sys` at `beta` with `basis_size` basis functions
and `P` links.
"""
function Trotter(sys::System, beta::Float64, basis_size::Int, P::Int)
    basis = Basis(sys, basis_size)
    h0, V = operators(basis, sys)

    # Single Trotter product.
    tau = beta / P
    rho = exp(-tau * h0) * exp(-tau * V)

    # Full path.
    path = rho^P

    Z = tr(path)
    E = tr(path * (h0 + V)) / Z
    Cv = (tr(path * (h0 + V)^2) / Z - E^2) * beta^2

    Trotter(Z, E, Cv)
end


text:
# ---
# title: Colorant conversion and channel layout
# cover: layout.gif
# description: a set of commonly used basic operations that wrapped by Augmentor
# ---

# Augmentor has warpped some commonly used basic operations so that you can use to build the augmentation
# pipeline. The `internal` column is what you'd probably do outside of `Augmentor`.

# | Category             | internal                | Augmentor          |
# | ---                  | ---                     | ---                |
# | Conversion           | `T.(img)`               | `ConvertEltype(T)` |
# | Information Layout   | `ImageCore.channelview` | `SplitChannels`    |
# | Information Layout   | `ImageCore.colorview`   | `CombineChannels`  |
# | Information Layout   | `Base.permutedims`      | `PermuteDims`      |
# | Information Layout   | `Base.reshape`          | `Reshape`          |

# It is not uncommon that machine learning frameworks require the data in a specific form and layout.
# For example many deep learning frameworks expect the colorchannel of the images to be encoded in
# the third dimension of a 4-dimensional array. Augmentor allows to convert from (and to) these
# different layouts using special operations that are mainly useful in the beginning or end of a
# augmentation pipeline.

using Augmentor
using ImageCore

## 300×400 Matrix{RGB{N0f8}, 2} => 300×400×3 Array{Float32, 3}
img = testpattern(RGB, ratio=0.5)
img_in = augment(img, SplitChannels() |> PermuteDims(2, 3, 1) |> ConvertEltype(Float32))

## 300×400×3 Array{Float32, 3} => 300×400 Matrix{RGB{N0f8}, 2}
img_out = augment(img_in, ConvertEltype(N0f8) |> PermuteDims(3, 1, 2) |> CombineChannels(RGB))

img_out == img

# ## References

#md # ```@docs
#md # ConvertEltype
#md # SplitChannels
#md # CombineChannels
#md # PermuteDims
#md # Reshape
#md # ```

## save covers #src
using ImageMagick #src
include(joinpath("..", "assets", "utilities.jl")) #src
cover = make_gif(testpattern(RGB, ratio=0.5), ConvertEltype(Gray{N0f8}), 2) #src
ImageMagick.save("layout.gif", cover; fps=1) #src


text:
pdd = ProbabilityDensityDistribution(["hej","farvel"],[x->x+1,x -> x+1])

@test ["hej","farvel"] == states(pdd)
@test_throws pdd = ProbabilityDensityDistribution(["hej","farvel","måske"],[x->x+1,x -> x+1])
@test_throws pdd = ProbabilityDensityDistribution(["hej","farvel","måske"],[x->x+1,x -> x+1, x -> "hej"])