text:
# In this file, we figure out how to link to Python (surprisingly complicated)
# and generate a deps/deps.jl file with the libpython name and other information
# needed for static compilation of PyReCall.

# As a result, if you switch to a different version or path of Python, you
# will probably need to re-run Pkg.build("PyReCall").

using Compat
import Conda

immutable UseCondaPython <: Exception end

try # make sure deps.jl file is removed on error

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

# Fix the environment for running `python`, and setts IO encoding to UTF-8.
# If cmd is the Conda python, then additionally removes all REPYTHON* and
# CONDA* environment variables.
function pythonenv(cmd::Cmd)
    env = copy(ENV)
    if dirname(cmd.exec[1]) == abspath(Conda.PYTHONDIR)
        pythonvars = String[]
        for var in keys(env)
            if startswith(var, "CONDA") || startswith(var, "REPYTHON")
                push!(pythonvars, var)
            end
        end
        for var in pythonvars
            pop!(env, var)
        end
    end
    # set PYTHONIOENCODING when running python executable, so that
    # we get UTF-8 encoded text as output (this is not the default on Windows).
    env["PYTHONIOENCODING"] = "UTF-8"
    setenv(cmd, env)
end

pyvar(python::AbstractString, mod::AbstractString, var::AbstractString) = chomp(readstring(pythonenv(`$python -c "import $mod; print($mod.$var)"`)))

pyconfigvar(python::AbstractString, var::AbstractString) = pyvar(python, "distutils.sysconfig", "get_config_var('$var')")
pyconfigvar(python, var, default) = let v = pyconfigvar(python, var)
    v == "None" ? default : v
end

pysys(python::AbstractString, var::AbstractString) = pyvar(python, "sys", var)

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

const dlprefix = is_windows() ? "" : "lib"

# return libpython name, libpython pointer
function find_libpython(python::AbstractString)
    # it is ridiculous that it is this hard to find the name of libpython
    v = pyconfigvar(python,"VERSION","")
    libs = [ dlprefix*"python"*v, dlprefix*"python" ]
    lib = pyconfigvar(python, "LIBRARY")
    lib != "None" && unshift!(libs, splitext(lib)[1])
    lib = pyconfigvar(python, "LDLIBRARY")
    lib != "None" && unshift!(unshift!(libs, basename(lib)), lib)
    libs = unique(libs)

    # it is ridiculous that it is this hard to find the path of libpython
    libpaths = [pyconfigvar(python, "LIBDIR"),
                (is_windows() ? dirname(pysys(python, "executable")) : joinpath(dirname(dirname(pysys(python, "executable"))), "lib"))]
    if is_apple()
        push!(libpaths, pyconfigvar(python, "PYTHONFRAMEWORKPREFIX"))
    end

    # `prefix` and `exec_prefix` are the path prefixes where python should look for python only and compiled libraries, respectively.
    # These are also changed when run in a virtualenv.
    exec_prefix = pysys(python, "exec_prefix")

    push!(libpaths, exec_prefix)
    push!(libpaths, joinpath(exec_prefix, "lib"))

    error_strings = String[]

    # TODO: other paths? python-config output? pyconfigvar("LDFLAGS")?

    # find libpython (we hope):
    for lib in libs
        for libpath in libpaths
            libpath_lib = joinpath(libpath, lib)
            if isfile(libpath_lib*"."*Libdl.dlext)
                try
                    return (Libdl.dlopen(libpath_lib,
                                         Libdl.RTLD_LAZY|Libdl.RTLD_DEEPBIND|Libdl.RTLD_GLOBAL),
                            libpath_lib)
                catch e
                    push!(error_strings, string("dlopen($libpath_lib) ==> ", e))
                end
            end
        end
    end

    # We do this *last* because the libpython in the system
    # library path might be the wrong one if multiple python
    # versions are installed (we prefer the one in LIBDIR):
    for lib in libs
        lib = splitext(lib)[1]
        try
            return (Libdl.dlopen(lib, Libdl.RTLD_LAZY|Libdl.RTLD_DEEPBIND|Libdl.RTLD_GLOBAL),
                    lib)
        catch e
            push!(error_strings, string("dlopen($lib) ==> ", e))
        end
    end

    if "yes" == get(ENV, "PyReCall_DEBUG_BUILD", "no") # print out extra info to help with remote debugging
        println(STDERR, "------------------------------------- exceptions -----------------------------------------")
        for s in error_strings
            print(s, "\n\n")
        end
        println(STDERR, "---------------------------------- get_config_vars ---------------------------------------")
        print(STDERR, readstring(`python -c "import distutils.sysconfig; print(distutils.sysconfig.get_config_vars())"`))
        println(STDERR, "--------------------------------- directory contents -------------------------------------")
        for libpath in libpaths
            if isdir(libpath)
                print(libpath, ":\n")
                for file in readdir(libpath)
                    if contains(file, "pyth")
                        println("    ", file)
                    end
                end
            end
        end
        println(STDERR, "------------------------------------------------------------------------------------------")
    end

    error("""
        Couldn't find libpython; check your REPYTHON environment variable.

        The python executable we tried was $python (= version $v);
        the library names we tried were $libs
        and the library paths we tried were $libpaths""")
end

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

include("depsutils.jl")

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

const python = try
    let py = get(ENV, "REPYTHON", isfile("REPYTHON") ? readchomp("REPYTHON") :
                 is_linux() || Sys.ARCH ∉ (:i686, :x86_64) ? "python" : ""),
        vers = isempty(py) ? v"0.0" : convert(VersionNumber, pyconfigvar(py,"VERSION","0.0"))
        if vers < v"2.7"
            if isempty(py)
                throw(UseCondaPython())
            else
                error("Python version $vers < 2.7 is not supported")
            end
        end

        # check word size of Python via sys.maxsize, since a common error
        # on Windows is to link a 64-bit Julia to a 32-bit Python.
        pywordsize = parse(UInt64, pysys(py, "maxsize")) > (UInt64(1)<<32) ? 64 : 32
        if pywordsize != Sys.WORD_SIZE
            error("$py is $(pywordsize)-bit, but Julia is $(Sys.WORD_SIZE)-bit")
        end

        py
    end
catch e1
    if Sys.ARCH in (:i686, :x86_64)
        if isa(e1, UseCondaPython)
            info("Using the Python distribution in the Conda package by default.\n",
                 "To use a different Python version, set ENV[\"REPYTHON\"]=\"pythoncommand\" and re-run Pkg.build(\"PyReCall\").")
        else
            info( "No system-wide Python was found; got the following error:\n",
                  "$e1\nusing the Python distribution in the Conda package")
        end
        abspath(Conda.PYTHONDIR, "python" * (is_windows() ? ".exe" : ""))
    else
        error("No system-wide Python was found; got the following error:\n",
              "$e1")
    end
end

use_conda = dirname(python) == abspath(Conda.PYTHONDIR)
if use_conda
    Conda.add("numpy")
end

const (libpython, libpy_name) = find_libpython(python)
const programname = pysys(python, "executable")

# Get PYTHONHOME, either from the environment or from Python
# itself (if it is not in the environment or if we are using Conda)
PYTHONHOME = if !haskey(ENV, "PYTHONHOME") || use_conda
    # PYTHONHOME tells python where to look for both pure python
    # and binary modules.  When it is set, it replaces both
    # `prefix` and `exec_prefix` and we thus need to set it to
    # both in case they differ. This is also what the
    # documentation recommends.  However, they are documented
    # to always be the same on Windows, where it causes
    # problems if we try to include both.
    exec_prefix = pysys(python, "exec_prefix")
    is_windows() ? exec_prefix : pysys(python, "prefix") * ":" * exec_prefix
else
    ENV["PYTHONHOME"]
end

# cache the Python version as a Julia VersionNumber
const pyversion = VersionNumber(pyvar(python, "platform", "python_version()"))

info("PyReCall is using $python (Python $pyversion) at $programname, libpython = $libpy_name")

if pyversion < v"2.7"
    error("Python 2.7 or later is required for PyReCall")
end

# A couple of key strings need to be stored as constants so that
# they persist throughout the life of the program.  In Python 3,
# they need to be wchar_t* data.
wstringconst(s) =
    VERSION < v"0.5.0-dev+4859" ?
    string("wstring(\"", escape_string(s), "\")") :
    string("Base.cconvert(Cwstring, \"", escape_string(s), "\")")

# we write configuration files only if they change, both
# to prevent unnecessary recompilation and to minimize
# problems in the unlikely event of read-only directories.
function writeifchanged(filename, str)
    if !isfile(filename) || readstring(filename) != str
        info(abspath(filename), " has been updated")
        write(filename, str)
    else
        info(abspath(filename), " has not changed")
    end
end

writeifchanged("deps.jl", """
    const python = "$(escape_string(python))"
    const libpython = "$(escape_string(libpy_name))"
    const pyprogramname = "$(escape_string(programname))"
    const wpyprogramname = $(wstringconst(programname))
    const pyversion_build = $(repr(pyversion))
    const PYTHONHOME = "$(escape_string(PYTHONHOME))"
    const wPYTHONHOME = $(wstringconst(PYTHONHOME))

    "True if we are using the Python distribution in the Conda package."
    const conda = $use_conda
    """)

# Make subsequent builds (e.g. Pkg.update) use the same Python by default:
writeifchanged("REPYTHON", isfile(programname) ? programname : python)

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

catch

# remove deps.jl (if it exists) on an error, so that PyReCall will
# not load until it is properly configured.
isfile("deps.jl") && rm("deps.jl")
rethrow()

end


text:
using Flux
using Flux.Losses: mse
using Flux.Data: DataLoader

function build_model(; hidden_layer_size::Int64 = 32)
    return Chain(Dense(2, hidden_layer_size, relu), Dense(hidden_layer_size, 2))
end

function get_data(
    x_array::Array{Array{Float64}},
    disp_array::Array{Array{Float64}};
    batchsize::Int64 = 256,
)
    Xtrain = hcat(x_array...)
    Ytrain = hcat(disp_array...)

    DataLoader((Xtrain, Ytrain), batchsize = batchsize, shuffle = true)
end

function loss_fn(data, model)
    loss = 0.0
    num = 0
    for (x, y) in data
        ŷ = model(x)
        loss += mse(ŷ, y)
        num += size(x)[end]
    end
    loss / num
end

function train(
    x_array::Array{Array{Float64}},
    disp_array::Array{Array{Float64}};
    lr::Float64 = 1e-3,
    batchsize = 256,
    epochs::Int64 = 10,
)
    data = get_data(x_array, disp_array, batchsize = batchsize)
    model = build_model() |> cpu
    ps = Flux.params(model)
    opt = ADAM(lr)

    for epoch = 1:epochs
        for (x, y) in data
            gs = gradient(() -> mse(model(x), y), ps)
            Flux.Optimise.update!(opt, ps, gs)
        end
        loss = loss_fn(data, model)
        println("Epoch = ", epoch, " Loss = ", loss)
    end
    model
end

function predict(model, x::Array{Float64})::Array{Float64}
    x + model(x)
end

function fit_ensemble(
    x_array::Array{Array{Float64}},
    disp_array::Array{Array{Float64}};
    ensemble_size::Int64 = 10,
)
    models = []
    for ensemble=1:ensemble_size
        println("Fitting network ", ensemble, " in the ensemble")
        push!(models, train(x_array, disp_array))
    end
    models
end

function predict_ensemble(
    models,
    x::Array{Float64},
)::Array{Array{Float64}}
    predictions = []
    for model in models
        push!(predictions, predict(model, x))
    end
    predictions
end

function fit_ensembles(
    x_array::Array{Array{Array{Float64}}},
    disp_array::Array{Array{Array{Float64}}},
)
    ensembles = []
    n_actions = length(x_array)
    for a=1:n_actions
        println("Fitting ensembles for action ", a)
        push!(ensembles, fit_ensemble(x_array[a], disp_array[a]))
    end
    ensembles
end

function find_max_distance(
    predictions::Array{Array{Float64}}
)
    max_distance = 0.0
    for i=1:length(predictions)
        for j=i+1:length(predictions)
            x = predictions[i]
            y = predictions[j]
            max_distance = max(max_distance, norm(x - y))
        end
    end
    max_distance
end

ext: jl


text:
using PerceptualColourMaps
using PerceptualColourMaps: ch2ab, newcolourmapdef
using Printf
using PyPlot
using ColorTypes

#=
Code for generating a set of different colormaps from pre-selected
base colors. Totally based on the code from PerceptualColourMaps
package
=#

function buildMap(
    CM;
    N::Int = 256,
    chromaK::Real = 1.0,
    shift::Real = 0.0,
    reverse::Bool = false,
    diagnostics::Bool = false,
    returnname = false,
    equalize = true,
)
    if uppercase(CM.colourspace) == "LAB"
        CM.colpts[:, 2:3] = chromaK * CM.colpts[:, 2:3]
    end

    # Colour map path is formed via a b-spline in the specified colour space
    Npts = size(CM.colpts, 1)

    if Npts < 2
        println(I)
        error("Number of input points must be 2 or more")
    elseif Npts < CM.splineorder
        error("Spline order is greater than number of data points")
    end

    # Rely on the attribute string to identify if colour map is cyclic.  We may
    # want to construct a colour map that has identical endpoints but do not
    # necessarily want continuity in the slope of the colour map path.
    if occursin("cyclic", lowercase(CM.attributeStr))
        cyclic = true
        labspline = pbspline(CM.colpts', CM.splineorder, N)
    else
        cyclic = false
        labspline = bbspline(CM.colpts', CM.splineorder, N)
    end

    # Apply contrast equalisation with required parameters. Note that sigma is
    # normalised with respect to a colour map of length 256 so that if a short
    # colour map is specified the smoothing that is applied is adjusted to suit.
    CM.sigma = CM.sigma * N / 256
    if equalize
        map = equalisecolourmap(
            CM.colourspace,
            labspline',
            CM.formula,
            CM.W,
            CM.sigma,
            cyclic,
            diagnostics,
        )
    else
        map = labspline'
    end

    # If specified apply a cyclic shift to the colour map
    if shift != 0
        if !occursin("cyclic", lowercase(CM.attributeStr))
            @warn "Colour map shifting being applied to a non-cyclic map"
        end
        map = circshift(map, round(Int, N * shift))
    end

    if reverse
        map = Base.reverse(map, dims = 1)
    end

    # Compute mean chroma of colour map for use in the name construction
    lab = srgb2lab(map)
    meanchroma = sum(sqrt.(sum(lab[:, 2:3] .^ 2, dims = 2))) / N

    # Construct lightness range description
    if uppercase(CM.colourspace) == "LAB"  # Use the control points
        L = CM.colpts[:, 1]
    else  # For RGB use the converted CIELAB values
        L = round.(lab[:, 1])
    end
    minL = minimum(L)
    maxL = maximum(L)

    if minL == maxL     # Isoluminant
        LStr = @sprintf("%d", minL)

    elseif L[1] == maxL && (
        occursin("diverging", lowercase(CM.attributeStr)) ||
        occursin("linear", lowercase(CM.attributeStr))
    )
        LStr = @sprintf("%d-%d", maxL, minL)
    else
        LStr = @sprintf("%d-%d", minL, maxL)
    end

    # Build overall colour map name
    name = @sprintf(
        "%s_%s_%s_c%d_n%d",
        CM.attributeStr,
        CM.hueStr,
        LStr,
        Int(round(meanchroma)),
        N
    )

    if shift != 0
        name = @sprintf("%s_s%d", name, round(Int, shift * 100))
    end

    if reverse
        name = @sprintf("%s_r", name)
    end

    # Build an array of ColorTypes.RGB values to return
    rgbmap = Array{ColorTypes.RGBA{Float64}}(undef, N)
    for i = 1:N
        rgbmap[i] = ColorTypes.RGBA(map[i, 1], map[i, 2], map[i, 3], 1.0)
    end

    if returnname
        return rgbmap, name, CM.desc
    else
        return rgbmap
    end
end

# define the base colors
base_colors = Dict(
    "mnfred" => [157 13 21] / 255.0,
    "blue" => [68 114 196] / 255.0,
    "darkpurple" => [86 71 148] / 255.0,
    "darkyellow" => [255 192 0] / 255.0,
    "orange" => [237 125 49] / 255.0,
    "green" => [112 173 71] / 255.0,
)

function saveCPT(cm, fname, N = 256)
    # saves the generated colormap to CPT format,
    # readable by python
    cmap = ColorMap(cm)
    b = Int.(round.(cmap(0.0) .* 255))[1:3]
    f = Int.(round.(cmap(1.0) .* 255))[1:3]
    na = [0, 0, 0]
    vals = collect(range(0.0, 1.0, length = N))
    cols =
        Int.(
            round.(
                cmap(collect(range(0.0, 1.0, length = 256)))[:, 1:3] * 255.0,
            ),
        )
    arr = [vals[1:end-1] cols[1:end-1, :] vals[2:end] cols[2:end, :]]
    open(fname, "w") do io
        println(io, "# COLOR_MODEL = RGB")
        for row in eachrow(arr)
            println(io, @sprintf("%e %3d %3d %3d %e %3d %3d %3d", row...))
        end
        println(
            io,
            @sprintf(
                "B %3d %3d %3d\nF %3d %3d %3d\nN %3d %3d %3d",
                b...,
                f...,
                na...
            )
        )
    end
end

function showMap(mymap; cyclic = false)
    clf()
    if cyclic
        sr, alpha = circlesineramp()
        imshow(sr, cmap = ColorMap(mymap))
    else
        sr = sineramp()
        imshow(sr, cmap = ColorMap(mymap), aspect = 0.2)
    end
    axis("off")
    gcf()
end
#---
# classical hot colormap, slightly adjusted to university colors
quanthot = newcolourmapdef(
    desc = "quant-hot",
    hueStr = "quanthot",
    attributeStr = "linear",
    colourspace = "RGB",
    colpts = [
        0.0 0.0 0.0
        base_colors["mnfred"]
        1.0 0.15 0.0
        base_colors["orange"]
        base_colors["darkyellow"]
        1.0 1.0 1.0
    ],
    splineorder = 2,
    formula = "CIE76",
    W = [1, 0, 0],
    sigma = 0,
)
#---
# derivative of "parula-like" colormap from PerceptualColourMaps
# adjusted to university colors
quantparrot = newcolourmapdef(
    desc = "quant-parrot",
    attributeStr = "linear",
    hueStr = "quantparrot",
    colourspace = "LAB",
    colpts = [
        20 0 0
        srgb2lab(base_colors["blue"])
        srgb2lab(base_colors["green"])
        srgb2lab(base_colors["darkyellow"])
        95 -21 92
    ],
    splineorder = 2,
    formula = "CIE76",
    W = [1, 0, 0],
    sigma = 10,
)

#---
# perceptually unifrom colormap going from
# blue to yellow in the LAB colorspace
# visually reminds plasma

function rgb2lch(rgb)
    l, a, b = srgb2lab(rgb)
    lch = [l, sqrt(a^2 + b^2), rad2deg(atan(b, a))]
    return lch
end

function norm_range(l, newmin, newmax)
    if length(l) > 1
        l .-= minimum(l)
        l ./= maximum(l)
        l .*= (newmax - newmin)
        l .+= newmin
    end
    return l
end

# Map with linear decreasing lightness, increasing chroma to blue.
# Precompute a spiral of increasing chroma down through Lab space
nsteps = 12
# starting from hue of yellow base color
ang1 = rgb2lch(base_colors["darkyellow"])[3]
# going to the hue of blue base color
ang2 = rgb2lch(base_colors["blue"])[3] + 20.0
# Linearly interpolate hue angle
ang = range(ang1, stop = ang2, length = nsteps)

# Interpolate chroma but use a 'gamma' of 0.5 to keep the colours more saturated.
sat1 = 0;
sat2 = 80;
sat = (range(sat1, stop = sat2, length = nsteps) / sat2) .^ 0.5 * sat2

l1 = 100;
l2 = 15;        # Linearly interpoate lightness
l = range(l1, stop = l2, length = nsteps) .^ 1.0
l = norm_range(l, l2, l1)
colptarray = zeros(nsteps, 3)
for n = 1:nsteps
    colptarray[n, :] = [l[n] ch2ab(sat[n], ang[n])]
end

quantplasma = newcolourmapdef(
    desc = "quant-plasma",
    attributeStr = "linear",
    hueStr = "quantplasma",
    colourspace = "LAB",
    colpts = colptarray[end:-1:1, :],
    splineorder = 2,
    formula = "CIE76",
    W = [1, 0, 0],
    sigma = 0,
)
#---
# simple diverging colormap

quantcoolwarm = newcolourmapdef(
    desc = "quant-coolwarm",
    attributeStr = "diverging",
    hueStr = "quantcoolwarm",
    colourspace = "LAB",
    colpts = [
        40 ch2ab(83, rgb2lch(base_colors["blue"])[3])
        95 0 0
        40 ch2ab(83, rgb2lch(base_colors["mnfred"])[3])
    ],
    splineorder = 2,
    formula = "CIE76",
    W = [1, 0, 0],
    sigma = 10,
)

#---
# define helper functions needed to linearly interpolate
# colors in the RGB space
function genCoefs(positions)
    K = 1
    coefs = []
    totlen = positions[end]
    push!(
        coefs,
        [
            collect(range(0.0, length = positions[1]))[end:-1:1]...
            zeros(totlen - positions[1])
        ],
    )
    for i = 1:length(positions)-1
        if i > 1
            prev = positions[i-1]
        else
            prev = 0
        end
        push!(
            coefs,
            [
                zeros(prev)...
                norm_range(
                    collect(range(0.0, length = positions[i] - prev)),
                    0.0,
                    1.0,
                )...
                norm_range(
                    collect(
                        range(0.0, length = positions[i+1] - positions[i]),
                    )[end:-1:1],
                    0.0,
                    1.0,
                )...
                zeros(totlen - positions[i+1])
            ],
        )
    end
    push!(
        coefs,
        [
            zeros(positions[end-1])...
            collect(range(0.0, length = totlen - positions[end-1]))...
        ],
    )
    coefs = Base.map(x -> norm_range(x, 0.0, 1.0), coefs)
    return coefs
end

function interpolateRGB(colors, positions)
    #=
    colors - colors to be interpolated
    positions - positions of color changes, N-1 for N colors
    0 - 100% C1
    P1 - 100% C2
    P2 - 100% C3
    ...
    P_N - 100%
    =#
    coefs = genCoefs(positions)
    res = zeros(coefs[1] |> length, 3)
    for (color, coef) in zip(colors, coefs)
        res += color .* coef
    end
    return res
end

#---
# simple black-white-red colormap, interpolated in RGB space
# not perceptually uniiform
npoints = 105
start = base_colors["blue"]
center = [0.95 0.95 0.95] .* 1.0
stop = base_colors["mnfred"]
coef1 = [collect(range(0.0, length = npoints ÷ 2))[end:-1:1]... zeros(
    npoints ÷ 2,
)...]
K = 0.7
coef1 .^= K
coef1 = norm_range(coef1, 0.0, 1.0)
coef2 = [zeros(npoints ÷ 2)... collect(range(0.0, length = npoints ÷ 2))...]
coef2 .^= K

coef2 = norm_range(coef2, 0.0, 1.0)
cpoints =
    (start .* coef1') .+ (stop .* coef2') .+
    (center .* (1.0 .- coef1 .- coef2)')

quantbwr = newcolourmapdef(
    desc = "quant-bwr",
    attributeStr = "diverging",
    hueStr = "quantbwr",
    colourspace = "RGB",
    colpts = cpoints,
    splineorder = 2,
    formula = "CIE76",
    W = [0, 0, 0], # disable "uniformization"
    sigma = 0,
)

#---
# jet-like colormap, obtained by interpolating base colors
# in RGB space, not perceptually uniform

jetcolors = [
    base_colors[c] for
    c in ["darkpurple", "blue", "green", "darkyellow", "orange", "mnfred"]
]
cpoints = interpolateRGB(jetcolors |> values, [6, 12, 18, 24, 30])
quantjet = newcolourmapdef(
    desc = "",
    attributeStr = "linear",
    hueStr = "quantjet",
    colourspace = "RGB",
    colpts = cpoints,
    splineorder = 2,
    formula = "CIE76",
    W = [0, 0, 0],
    sigma = 0.0,
)
#---
# jet-like colormap, obtained by interpolating base colors
# in RGB space, perc. uniform
quantjetuniform = newcolourmapdef(
    desc = "",
    attributeStr = "linear",
    hueStr = "quantjetuniform",
    colourspace = "RGB",
    colpts = cpoints,
    splineorder = 2,
    formula = "CIE76",
    W = [1, 0, 0],
    sigma = 0.0,
)

#---
# jet-like colormap, obtained by interpolating base colors
# in RGB space, perc. uniform
cyclic_jetcolors = [
    base_colors[c] for c in [
        "darkpurple",
        "blue",
        "green",
        "darkyellow",
        "orange",
        "mnfred",
        "darkpurple",
    ]
]
cyclic_cpoints =
    interpolateRGB(cyclic_jetcolors |> values, [8, 14, 22, 28, 34, 40])

quantcyclicjet = newcolourmapdef(
    desc = "",
    attributeStr = "cyclic",
    hueStr = "quantcyclicjet",
    colourspace = "RGB",
    colpts = cyclic_cpoints,
    splineorder = 5,
    formula = "CIE76",
    W = [0.0, 0, 0],
    sigma = 0.0,
)

#showMap(buildMap(quantcyclicjet, shift=0.0), cyclic=true)

#---
# iterate over all defined colormaps and save them to CPT
for mapdef in [
    quantbwr,
    quanthot,
    quantparrot,
    quantplasma,
    quantcoolwarm,
    quantjetuniform,
    quantjet,
    quantcyclicjet,
]
    mymap, name = buildMap(mapdef, returnname = true)
    saveCPT(mymap, "$(name).cpt")
end

#---
#generate linear colormaps, going each of the base colors to white

for (idx, (color, rgb)) in enumerate(base_colors)
    mymapdef = newcolourmapdef(
        desc = "Grey scale",
        hueStr = "$(color)_white",
        attributeStr = "linear",
        colourspace = "LAB",
        colpts = [
                  [40 srgb2lab(rgb)[2:end]...]
                  100 0 0],
        splineorder = 2,
        formula = "CIE76",
        W = [1, 0, 0],
        sigma = 0,
    )

    mymap, name = buildMap(mymapdef, returnname = true)
    saveCPT(mymap, "$(name).cpt")
end

# and from each of the base colors to black
for (idx, (color, rgb)) in enumerate(base_colors)
    mymapdef = newcolourmapdef(
        desc = "Grey scale",
        hueStr = "black_$(color)",
        attributeStr = "linear",
        colourspace = "LAB",
        colpts = [
            0 0 0
            [70 srgb2lab(rgb)[2:end]...]
        ],
        splineorder = 2,
        formula = "CIE76",
        W = [1, 0, 0],
        sigma = 0,
    )

    mymap, name = buildMap(mymapdef, returnname = true)
    saveCPT(mymap, "$(name).cpt")

end


text:
@testset "Measurements" begin
    @testset "Wetted liquid solid" begin
        area_size = [1.0]
        maxheight = [1.0]
        height = zeros(5,5)
        height[2:3, 2:3] .= 1.0
        Swalbe.wetted!(area_size, maxheight, height, 1)
        @test area_size[1] == 4.0
        @test maxheight[1] == 1.0

        area_size = ones(5)
        height = zeros(5,5)
        maxheight = zeros(5)
        for t in 1:5
            height[t, t] = t
            Swalbe.wetted!(area_size, maxheight, height, t)
        end
        sol1 = [1.0, 2.0, 3.0, 4.0, 5.0]

        for t in 1:5
            @test area_size[t] == sol1[t]
            @test maxheight[t] == sol1[t]
        end
    end

    @testset "Position liquid" begin
        height = zeros(5,5)
        fluid = falses(1,25)
        height[2:3, 2:3] .= 1.0
        dummy = falses(5,5)
        Swalbe.fluid_dry!(fluid, dummy, height, 1)

        @test all(fluid[1, 1:6] .== false)
        @test all(fluid[1, 7:8] .== true)

        fluid = falses(5,25)
        height = zeros(5,5)
        maxheight = zeros(5)
        # Desired solution
        sol2 = [1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
                1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
                1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0;
                1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0;
                1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1]
        for t in 1:5
            height[t, t] = t
            Swalbe.fluid_dry!(fluid, dummy, height, t)
        end
        # println("This is fluid $fluid\nAnd sol $(reshape(sol2, 5, 25))")
        for t in 1:5
            @test all(fluid[t,:] .== sol2[t,:])
        end
    end

    @testset "Surface area liquid vapor and reduced energy" begin
        height = reshape(collect(1.0:25.0),5,5)
        dummy = ones(5,5,9)
        surface = ones(5,5)
        dx = ones(5,5)
        dy = ones(5,5)
        area_lv = [0.0]
        red_energy = [0.0]
        θ = fill(0.5, 5, 5)
        Swalbe.surfacearea!(area_lv, red_energy, height, θ, dx, dy, dummy, surface, 1)
        # Solutions gradient
        Xgrad = [-1.5 -1.5 -1.5 -1.5 -1.5;
                  1.0 1.0 1.0 1.0 1.0;
                  1.0 1.0 1.0 1.0 1.0;
                  1.0 1.0 1.0 1.0 1.0;
                 -1.5 -1.5 -1.5 -1.5 -1.5]

        Ygrad = [-7.5 5.0 5.0 5.0 -7.5;
                 -7.5 5.0 5.0 5.0 -7.5;
                 -7.5 5.0 5.0 5.0 -7.5;
                 -7.5 5.0 5.0 5.0 -7.5;
                 -7.5 5.0 5.0 5.0 -7.5]

        solsurf = 0.0
        solsurf = sum(sqrt.(Xgrad.^2 .+ Ygrad.^2 .+ 1))
        solener = 0.0
        solener = solsurf - sum(cospi.(θ))
        println("area is $(area_lv)")
        @test solsurf ≈ area_lv[1] atol=1e-10
        @test solener ≈ red_energy[1] atol=1e-10
        # With a non vanishing cosine
        θ = fill(1/9, 5, 5)
        Swalbe.surfacearea!(area_lv, red_energy, height, θ, dx, dy, dummy, surface, 1)
        solener = solsurf - sum(cospi.(θ))
        @test solsurf ≈ area_lv[1] atol=1e-10
        @test solener ≈ red_energy[1] atol=1e-10
        # More time steps
        d1 = zeros(5,5)
        d2 = zeros(5,5)
        area_lv = zeros(5)
        red_energy = zeros(5)
        surf_cal = zeros(5)
        ener_cal = zeros(5)
        for t in 1:5
            height[t,t] = t
            θ[t,t] = 0.0
            Swalbe.∇f_simple!(d1, d2, height, dummy)
            Swalbe.surfacearea!(area_lv, red_energy, height, θ, dx, dy, dummy, surface, t)
            surf_cal[t] = sum(sqrt.(d1.^2 .+ d2.^2 .+ 1))
            ener_cal[t] = surf_cal[t] - sum(cospi.(θ))
        end
        @test all(area_lv .== surf_cal)
        @test all(red_energy .== ener_cal)
    end

    @testset "Characteristic time scale" begin
        t_0 = Swalbe.t0()
        @test t_0 == 177428.32340802532
    end

    @testset "Field snapshots" begin
        h1 = reshape(collect(1:25),5,5)
        h2 = reshape(collect(5:5:125),5,5)
        snapshot = zeros(2, 25);
        Swalbe.snapshot!(snapshot,h1,10,dumping=10)
        Swalbe.snapshot!(snapshot,h2,20,dumping=10)
        @test all(h1 .== reshape(snapshot[1,:],5,5))
        @test all(h2 .== reshape(snapshot[2,:],5,5))
    end
end


text:
push!(LOAD_PATH, "@stdlib")  # NOTE: this is needed to enable this test to run from the Pkg manager
push!(LOAD_PATH, joinpath(@__DIR__, ".."))
using Test
using Pkg
import ParallelStencil: SUPPORTED_PACKAGES, PKG_CUDA
TEST_PACKAGES = SUPPORTED_PACKAGES
@static if PKG_CUDA in TEST_PACKAGES
    import CUDA
    if !CUDA.functional() TEST_PACKAGES = filter!(x->x≠PKG_CUDA, TEST_PACKAGES) end
end

@static for package in TEST_PACKAGES  eval(:(
    @testset "$(basename(@__FILE__)) (package: $(nameof($package)))" begin
        @testset "incremental compilation" begin
            Pkg.activate(joinpath(@__DIR__, "test_projects", "Diffusion3D_$(nameof($package))"))
            using $(Symbol("Diffusion3D_$package"))
            @test $(Symbol("Diffusion3D_$package")).diffusion3D()
        end;
    end;
)) end == nothing || true;


text:
using VortexDistributions, Distributions, Statistics, BenchmarkTools

include("../util.jl")

##
N = 2^4
Lx = 200; Ly = Lx
x = LinRange(-Lx/2,Lx/2, N+1)[1:end-1]; y = LinRange(-Ly/2,Ly/2, N+1)[1:end-1];
psi0 = one.(x*y') |> complex
psi4 = Torus(psi0,x,y)

##
N=2^5
x = LinRange(-Lx/2,Lx/2, N+1)[1:end-1]; y = LinRange(-Ly/2,Ly/2, N+1)[1:end-1];
psi0 = one.(x*y') |> complex
psi5 = Torus(psi0,x,y)

##
N=2^6
x = LinRange(-Lx/2,Lx/2, N+1)[1:end-1]; y = LinRange(-Ly/2,Ly/2, N+1)[1:end-1];
psi0 = one.(x*y') |> complex
psi6 = Torus(psi0,x,y)

##
N=2^7
x = LinRange(-Lx/2,Lx/2, N+1)[1:end-1]; y = LinRange(-Ly/2,Ly/2, N+1)[1:end-1];
psi0 = one.(x*y') |> complex
psi7 = Torus(psi0,x,y)

##
N=2^8
x = LinRange(-Lx/2,Lx/2, N+1)[1:end-1]; y = LinRange(-Ly/2,Ly/2, N+1)[1:end-1];
psi0 = one.(x*y') |> complex
psi8 = Torus(psi0,x,y)

##
N=2^9
x = LinRange(-Lx/2,Lx/2, N+1)[1:end-1]; y = LinRange(-Ly/2,Ly/2, N+1)[1:end-1];
psi0 = one.(x*y') |> complex
psi9 = Torus(psi0,x,y)

##
N=2^10
x = LinRange(-Lx/2,Lx/2, N+1)[1:end-1]; y = LinRange(-Ly/2,Ly/2, N+1)[1:end-1];
psi0 = one.(x*y') |> complex
psi10 = Torus(psi0,x,y)

##
N=2^11
x = LinRange(-Lx/2,Lx/2, N+1)[1:end-1]; y = LinRange(-Ly/2,Ly/2, N+1)[1:end-1];
psi0 = one.(x*y') |> complex
psi11 = Torus(psi0,x,y)

##
N=2^12
x = LinRange(-Lx/2,Lx/2, N+1)[1:end-1]; y = LinRange(-Ly/2,Ly/2, N+1)[1:end-1];
psi0 = one.(x*y') |> complex
psi12 = Torus(psi0,x,y)

##
rand_vorts = rand_pointvortex(2, psi)

vortex!(psi4, rand_vorts)
vortex!(psi5, rand_vorts)
vortex!(psi6, rand_vorts)
vortex!(psi7, rand_vorts)
vortex!(psi8, rand_vorts)
vortex!(psi9, rand_vorts)
vortex!(psi10, rand_vorts)
vortex!(psi11, rand_vorts)
vortex!(psi12, rand_vorts)

bench4 = @benchmark naivePlaquette(psi4.ψ, [psi4.x, psi4.y], 0)
bench5 = @benchmark naivePlaquette(psi5.ψ, [psi5.x, psi5.y], 0)
bench6 = @benchmark naivePlaquette(psi6.ψ, [psi6.x, psi6.y], 0)
bench7 = @benchmark naivePlaquette(psi7.ψ, [psi7.x, psi7.y], 0)
bench8 = @benchmark naivePlaquette(psi8.ψ, [psi8.x, psi8.y], 0)
bench9 = @benchmark naivePlaquette(psi9.ψ, [psi9.x, psi9.y], 0)
bench10 = @benchmark naivePlaquette(psi10.ψ, [psi10.x, psi10.y], 0)
bench11 = @benchmark naivePlaquette(psi11.ψ, [psi11.x, psi11.y], 0)
bench12 = @benchmark naivePlaquette(psi12.ψ, [psi12.x, psi12.y], 0)

bench_all = [bench4, bench5, bench6, bench7, bench8, bench9, bench10, bench11, bench12]
N_all = [2^x for x=4:12]

using Plots, LaTeXStrings
bench_all_mean = [mean(bench_all[i].times) for i in 1:length(bench_all)]
plot(log.(2, N_all), log.(bench_all_mean), marker='*', xlab = L"n = 2^x", ylab = L"time = ln(10^9 s)", title = "log-log plaquette n x n grid", label = "naive")

fench4 = @benchmark remove_vortices_edge(findvortices_jumps(psi4), psi4)
fench5 = @benchmark remove_vortices_edge(findvortices_jumps(psi5), psi5)
fench6 = @benchmark remove_vortices_edge(findvortices_jumps(psi6), psi6)
fench7 = @benchmark remove_vortices_edge(findvortices_jumps(psi7), psi7)
fench8 = @benchmark remove_vortices_edge(findvortices_jumps(psi8), psi8)
fench9 = @benchmark remove_vortices_edge(findvortices_jumps(psi9), psi9)
fench10 = @benchmark remove_vortices_edge(findvortices_jumps(psi10), psi10)
fench11 = @benchmark remove_vortices_edge(findvortices_jumps(psi11), psi11)
fench12 = @benchmark remove_vortices_edge(findvortices_jumps(psi12), psi12)

fench_all = [fench4, fench5, fench6, fench7, fench8, fench9, fench10, fench11, fench12]
fench_all_mean = [mean(fench_all[i].times) for i in 1:length(fench_all)]

plot!(log.(2, N_all), log.(fench_all_mean), marker='*', label = "optimised")


text:
function CIR_smoothing_logscale_internals(δ, γ, σ, λ, data; silence = false)
    β = γ/σ^2
    α = δ/2#Alternative parametrisation

    if !silence
        println("Filtering")
    end
    Λ_of_t, logwms_of_t, θ_of_t = filter_CIR_logweights(δ, γ, σ, λ, data; silence = silence)
    if !silence
        println("Cost to go")
    end
    Λ_tilde_prime_of_t, logwms_tilde_of_t, θ_tilde_prime_of_t = compute_all_log_cost_to_go_functions_CIR(δ, γ, σ, λ, data; silence = silence)

    return merge_filtering_and_cost_to_go_logweights_CIR(Λ_of_t, logwms_of_t, θ_of_t, Λ_tilde_prime_of_t, logwms_tilde_of_t, θ_tilde_prime_of_t, data, α, β)

end

function filter_CIR_logweights(δ, γ, σ, λ, data; silence = true)
    Λ_of_t, wms_of_t, θ_of_t = filter_CIR(δ, γ, σ, λ, data; silence = silence)
    logwms_of_t = wms_of_t |> convert_weights_to_logweights
    return Λ_of_t, logwms_of_t, θ_of_t
end

compute_all_log_cost_to_go_functions_CIR(δ, γ, σ, λ, data; silence = false) = compute_all_log_cost_to_go_functions_CIR_pruning(δ, γ, σ, λ, data, (x,y) -> (x,y); silence = silence)

function compute_all_log_cost_to_go_functions_CIR_pruning(δ, γ, σ, λ, data, pruning_function::Function; silence = false)
    times = data |> keys |> collect |> sort
    reversed_times = reverse(times)
    α = δ/2#Alternative parametrisation
    θ0_CIR = γ/σ^2
    Λ_tilde_prime_of_t = Dict()
    logwms_tilde_of_t = Dict()
    θ_tilde_prime_of_t = Dict()

    yT = data[times[end]]

    Λ_tilde_prime_kp2 = [0]
    Λ_tilde_kp1 = t_CIR(yT, Λ_tilde_prime_kp2)
    θ_tilde_prime_kp2 = θ0_CIR
    logwms_tilde_kp2 = [0.]

    for k in 2:length(reversed_times)
        if !silence
            println("(Cost to go) Step index: $k")
            println("Number of components: $(length(Λ_tilde_kp1))")
        end
        # Change of notation for clarity
        prev_t = reversed_times[k-1]
        t = reversed_times[k]
        Δk = prev_t - t
        ykp1 = data[prev_t]

        # New weight computation
        θ_tilde_kp1 = θ_tilde_k_from_θ_tilde_prime_kp1(ykp1, θ_tilde_prime_kp2)

        logwms_tilde_kp1, Λ_tilde_prime_kp1 = logwms_tilde_kp1_from_logwms_tilde_kp2_pruning(logwms_tilde_kp2, Λ_tilde_prime_kp2, θ_tilde_kp1, θ_tilde_prime_kp2, ykp1, Δk, α, γ, σ, λ; pruning_function = pruning_function, return_indices = true)

        #Storage of the results
        θ_tilde_prime_kp1 = θ_tilde_prime_k_from_θ_tilde_k_CIR(θ_tilde_kp1, Δk, γ, σ)

        yk = data[t]
        Λ_tilde_kp = Λ_tilde_k_from_Λ_tilde_prime_kp1_CIR(yk, Λ_tilde_prime_kp1) #Not stored, but better for consistency of notations.

        θ_tilde_prime_of_t[prev_t] = θ_tilde_prime_kp1
        Λ_tilde_prime_of_t[prev_t] = Λ_tilde_prime_kp1
        logwms_tilde_of_t[prev_t] = logwms_tilde_kp1

        #Preparation of next iteration
        θ_tilde_prime_kp2 = θ_tilde_prime_kp1
        logwms_tilde_kp2 = logwms_tilde_kp1
        Λ_tilde_prime_kp2 = Λ_tilde_prime_kp1
        Λ_tilde_kp1 = Λ_tilde_kp
    end

    return Λ_tilde_prime_of_t, logwms_tilde_of_t, θ_tilde_prime_of_t
end

function θ_tilde_k_from_θ_tilde_prime_kp1(yk, θ_tilde_prime_kp1; λ = 1)
    return T_CIR(yk, θ_tilde_prime_kp1; λ = λ)
end

function logwms_tilde_kp1_from_logwms_tilde_kp2_pruning(logwms_tilde_kp2, Λ_tilde_prime_kp2, θ_tilde_kp1, θ_tilde_prime_kp2, ykp1, Δk, α, γ, σ, λ; pruning_function::Function, return_indices = false)

    updated_logwms_tilde_kp2 = update_logweights_cost_to_go_CIR(logwms_tilde_kp2, Λ_tilde_prime_kp2, θ_tilde_prime_kp2, α, ykp1, λ)

    Λ_tilde_prime_kp2_pruned, updated_logwms_tilde_kp2_pruned = pruning_function(Λ_tilde_prime_kp2, updated_logwms_tilde_kp2)

    logwms_tilde_kp1 = predict_logweights_cost_to_go_CIR(updated_logwms_tilde_kp2_pruned, Λ_tilde_prime_kp2_pruned, θ_tilde_kp1, θ_tilde_prime_kp2, ykp1, Δk, α, γ, σ, λ)

    if return_indices
        return logwms_tilde_kp1, 0:maximum(t_CIR(ykp1, maximum(Λ_tilde_prime_kp2_pruned)))
    else
        return logwms_tilde_kp1
    end
end

function θ_tilde_prime_k_from_θ_tilde_k_CIR(θ_tilde_k, Δt, γ, σ)
    return θ_prime_from_θ_CIR(θ_tilde_k, Δt, γ, σ)
end

function Λ_tilde_k_from_Λ_tilde_prime_kp1_CIR(yk, Λ_tilde_prime_kp1)
    return t_CIR(yk, Λ_tilde_prime_kp1)
end

function update_logweights_cost_to_go_CIR(logwms_tilde_kp2, Λ_tilde_prime_kp2, θ_tilde_prime_kp2, α, ykp1, λ)
    return update_logweights_cost_to_go(logwms_tilde_kp2, Λ_tilde_prime_kp2, θ_tilde_prime_kp2, ykp1, (α = α, λ = λ), (m, θ, y, params) -> logμπh(m, θ, params[:α], y; λ = params[:λ]))
end

function predict_logweights_cost_to_go_CIR(updated_logwms_tilde_kp2, Λ_tilde_prime_kp2, θ_tilde_kp1, θ_tilde_prime_kp2, ykp1, Δk, α, γ, σ, λ)
    logwms_tilde_kp1 = fill(-Inf, maximum(t_CIR(ykp1, maximum(Λ_tilde_prime_kp2)))+1)

    p = γ/σ^2*1/(θ_tilde_kp1*exp(2*γ*Δk) + γ/σ^2 - θ_tilde_kp1)

    for k in eachindex(Λ_tilde_prime_kp2)
        n = Λ_tilde_prime_kp2[k]
        t_ykp1_n = t_CIR(ykp1, n)

        for m in 0:t_ykp1_n
            idx = m+1
            # Careful, logwms_tilde_kp2[k] has index k because we assume that the weights and the indices are in the same order.
            logwms_tilde_kp1[idx] = logaddexp(logwms_tilde_kp1[idx], updated_logwms_tilde_kp2[k] + logpmn_CIR(t_ykp1_n, m, p))
        end
    end
    return logwms_tilde_kp1
end

function Λ_tilde_prime_k_from_Λ_tilde_k_CIR(Λ_tilde_k)
    return 0:maximum(Λ_tilde_k)
end

function logpmn_CIR(m, n, p)
    return logpdf(Binomial(m, p), n)
end

function logpmn_CIR(m, n, Δt, θ, σ, γ)
    p = γ/σ^2*1/(θ*exp(2*γ*Δt) + γ/σ^2 - θ)
    return logpmn_CIR(m, n, p)
end

function logCmn_CIR(m, n, θback_k, θ_k, α, β)
    return lgamma_local(α) - lgamma_local(α+m) - lgamma_local(α+n) - α*log(β) + (α+m)*log(θback_k) + (α+n)*log(θ_k) + lgamma_local(α+m+n) - (α+m+n) * log(θback_k + θ_k-β)
end

function merge_filtering_and_cost_to_go_logweights_CIR(Λ_of_t, logwms_of_t, θ_of_t, Λ_tilde_prime_of_t, logwms_tilde_of_t, θ_tilde_prime_of_t, data, α, β)
    times = Λ_of_t |> keys |> collect |> sort

    Λ_of_t_smooth = Dict()
    wms_of_t_smooth = Dict()
    θ_of_t_smooth = Dict()

    log_Λ_weights = Array{Float64,1}(undef, 2*sum(sum(values(data))))

    for k in 1:(length(times)-1)
        fill!(log_Λ_weights, -Inf)

        Λk = Λ_of_t[times[k]]
        logwk = logwms_of_t[times[k]]
        Λ_tilde_prime_kp1 = Λ_tilde_prime_of_t[times[k+1]]
        logw_tilde_kp1 = logwms_tilde_of_t[times[k+1]]
        for i in eachindex(Λk)
            n = Λk[i]
            for j in eachindex(Λ_tilde_prime_kp1)
                m = Λ_tilde_prime_kp1[j]
                log_Λ_weights[d_CIR(m, n)] = logaddexp(log_Λ_weights[d_CIR(m, n)], logw_tilde_kp1[j] + logwk[i] + logCmn_CIR(m, n, θ_tilde_prime_of_t[times[k+1]], θ_of_t[times[k]], α, β))
            end
        end
        Λ_weights = exp.(log_Λ_weights .- logsumexp(log_Λ_weights))

        Λ_of_t_smooth[times[k]] = [n for n in eachindex(Λ_weights) if Λ_weights[n] > 0.]
        wms_of_t_smooth[times[k]] = Λ_weights[Λ_weights .> 0.]
        θ_of_t_smooth[times[k]] = e_CIR(θ_tilde_prime_of_t[times[k+1]], θ_of_t[times[k]], β)
    end

        #The last smoothing distribution is a filtering distribution
        Λ_of_t_smooth[times[end]] = Λ_of_t[times[end]]
        wms_of_t_smooth[times[end]] = exp.(logwms_of_t[times[end]])
        θ_of_t_smooth[times[end]] = θ_of_t[times[end]]

    return Λ_of_t_smooth, wms_of_t_smooth, θ_of_t_smooth
end

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering


text:
"""
    generate(pkg::AbstractString, t::Template) -> Nothing
    generate(t::Template, pkg::AbstractString) -> Nothing

Generate a package named `pkg` from `t`. If `git` is `false`, no Git repository is created.
"""
function generate(
    pkg::AbstractString,
    t::Template;
    git::Bool=true,
    gitconfig::Union{GitConfig, Nothing}=nothing,
)
    pkg = splitjl(pkg)
    pkg_dir = joinpath(t.dir, pkg)
    ispath(pkg_dir) && throw(ArgumentError("$pkg_dir already exists"))

    try
        # Create the directory with some boilerplate inside.
        Pkg.generate(pkg_dir)

        # Add a [compat] section for Julia.
        open(joinpath(pkg_dir, "Project.toml"), "a") do io
            println(io, "\n[compat]\njulia = $(repr_version(t.julia_version))")
        end

        # Replace the authors field with the template's authors.
        if !isempty(t.authors)
            path = joinpath(pkg_dir, "Project.toml")
            project = read(path, String)
            authors = string("[", join(map(repr ∘ strip, split(t.authors, ",")), ", "), "]")
            write(path, replace(project, r"authors = .*" => "authors = $authors"))
        end

        if git
            # Initialize the repo.
            repo = LibGit2.init(pkg_dir)
            @info "Initialized Git repo at $pkg_dir"

            if gitconfig !== nothing
                # Configure the repo.
                repoconfig = GitConfig(repo)
                for c in LibGit2.GitConfigIter(gitconfig)
                    LibGit2.set!(repoconfig, unsafe_string(c.name), unsafe_string(c.value))
                end
            end

            # Commit and set the remote.
            LibGit2.commit(repo, "Initial commit")
            rmt = if t.ssh
                "git@$(t.host):$(t.user)/$pkg.jl.git"
            else
                "https://$(t.host)/$(t.user)/$pkg.jl"
            end
            # We need to set the remote in a strange way, see #8.
            close(LibGit2.GitRemote(repo, "origin", rmt))
            @info "Set remote origin to $rmt"

            # Create the gh-pages branch if necessary.
            if haskey(t.plugins, GitHubPages)
                LibGit2.branch!(repo, "gh-pages")
                LibGit2.commit(repo, "Initial commit")
                @info "Created empty gh-pages branch"
                LibGit2.branch!(repo, "master")
            end
        end

        # Generate the files.
        files = vcat(
            "src/", "Project.toml",  # Created by Pkg.generate.
            gen_tests(pkg_dir, t),
            gen_readme(pkg_dir, t),
            gen_license(pkg_dir, t),
            vcat(map(p -> gen_plugin(p, t, pkg), values(t.plugins))...),
        )

        if git
            append!(files, gen_gitignore(pkg_dir, t))
            LibGit2.add!(repo, files...)
            LibGit2.commit(repo, "Files generated by PkgTemplates")
            @info "Committed $(length(files)) files/directories: $(join(files, ", "))"

            if length(collect(LibGit2.GitBranchIter(repo))) > 1
                @info "Remember to push all created branches to your remote: git push --all"
            end
        end

        if t.dev
            # Add the new package to the current environment.
            Pkg.develop(PackageSpec(path=pkg_dir))
        end

        @info "New package is at $pkg_dir"
    catch e
        rm(pkg_dir; recursive=true)
        rethrow(e)
    end
end

function generate(
    t::Template,
    pkg::AbstractString;
    git::Bool=true,
    gitconfig::Union{GitConfig, Nothing}=nothing,
)
    generate(pkg, t; git=git, gitconfig=gitconfig)
end

"""
    generate_interactive(pkg::AbstractString; fast::Bool=false, git::Bool=true) -> Template

Interactively create a template, and then generate a package with it. Arguments and
keywords are used in the same way as in [`generate`](@ref) and
[`interactive_template`](@ref).
"""
function generate_interactive(
    pkg::AbstractString;
    fast::Bool=false,
    git::Bool=true,
    gitconfig::Union{GitConfig, Nothing}=nothing,
)
    t = interactive_template(; git=git, fast=fast)
    generate(pkg, t; git=git, gitconfig=gitconfig)
    return t
end

"""
    gen_tests(pkg_dir::AbstractString, t::Template) -> Vector{String}

Create the test entrypoint in `pkg_dir`.

# Arguments
* `pkg_dir::AbstractString`: The package directory in which the files will be generated
* `t::Template`: The template whose tests we are generating.

Returns an array of generated file/directory names.
"""
function gen_tests(pkg_dir::AbstractString, t::Template)
    # TODO: Silence Pkg for this section? Adding and removing Test creates a lot of noise.
    proj = Base.current_project()
    try
        Pkg.activate(pkg_dir)
        Pkg.add("Test")

        # Move the Test dependency into the [extras] section.
        toml = read(joinpath(pkg_dir, "Project.toml"), String)
        lines = split(toml, "\n")
        idx = findfirst(l -> startswith(l, "Test = "), lines)
        testdep = lines[idx]
        deleteat!(lines, idx)
        toml = join(lines, "\n") * """
        [extras]
        $testdep

        [targets]
        test = ["Test"]
        """
        gen_file(joinpath(pkg_dir, "Project.toml"), toml)
        Pkg.update()  # Regenerate Manifest.toml (this cleans up Project.toml too).
    finally
        proj === nothing ? Pkg.activate() : Pkg.activate(proj)
    end

    pkg = basename(pkg_dir)
    text = """
        using $pkg
        using Test

        @testset "$pkg.jl" begin
            # Write your own tests here.
        end
        """

    gen_file(joinpath(pkg_dir, "test", "runtests.jl"), text)
    return ["test/"]
end

"""
    gen_readme(pkg_dir::AbstractString, t::Template) -> Vector{String}

Create a README in `pkg_dir` with badges for each enabled plugin.

# Arguments
* `pkg_dir::AbstractString`: The directory in which the files will be generated.
* `t::Template`: The template whose README we are generating.

Returns an array of generated file/directory names.
"""
function gen_readme(pkg_dir::AbstractString, t::Template)
    pkg = basename(pkg_dir)
    text = "# $pkg\n"
    done = []
    # Generate the ordered badges first, then add any remaining ones to the right.
    for plugin_type in BADGE_ORDER
        if haskey(t.plugins, plugin_type)
            text *= "\n"
            text *= join(
                badges(t.plugins[plugin_type], t.user, pkg),
                "\n",
            )
            push!(done, plugin_type)
        end
    end
    for plugin_type in setdiff(keys(t.plugins), done)
        text *= "\n"
        text *= join(
            badges(t.plugins[plugin_type], t.user, pkg),
            "\n",
        )
    end
    if haskey(t.plugins, Citation) && t.plugins[Citation].readme_section
        text *= "\n## Citing\n\nSee `CITATION.bib` for the relevant reference(s).\n"
    end

    gen_file(joinpath(pkg_dir, "README.md"), text)
    return ["README.md"]
end

"""
    gen_gitignore(pkg_dir::AbstractString, t::Template) -> Vector{String}

Create a `.gitignore` in `pkg_dir`.

# Arguments
* `pkg_dir::AbstractString`: The directory in which the files will be generated.
* `t::Template`: The template whose .gitignore we are generating.

Returns an array of generated file/directory names.
"""
function gen_gitignore(pkg_dir::AbstractString, t::Template)
    pkg = basename(pkg_dir)
    init = [".DS_Store", "/dev/"]
    entries = mapfoldl(p -> p.gitignore, append!, values(t.plugins); init=init)
    if !t.manifest && !in("Manifest.toml", entries)
        push!(entries, "/Manifest.toml")  # Only ignore manifests at the repo root.
    end
    unique!(sort!(entries))
    text = join(entries, "\n")

    gen_file(joinpath(pkg_dir, ".gitignore"), text)
    files = [".gitignore"]
    t.manifest && push!(files, "Manifest.toml")
    return files
end

"""
    gen_license(pkg_dir::AbstractString, t::Template) -> Vector{String}

Create a license in `pkg_dir`.

# Arguments
* `pkg_dir::AbstractString`: The directory in which the files will be generated.
* `t::Template`: The template whose LICENSE we are generating.

Returns an array of generated file/directory names.
"""
function gen_license(pkg_dir::AbstractString, t::Template)
    if isempty(t.license)
        return String[]
    end

    text = "Copyright (c) $(year(today())) $(t.authors)\n"
    text *= read_license(t.license)

    gen_file(joinpath(pkg_dir, "LICENSE"), text)
    return ["LICENSE"]
end

"""
    gen_file(file::AbstractString, text::AbstractString) -> Int

Create a new file containing some given text. Always ends the file with a newline.

# Arguments
* `file::AbstractString`: Path to the file to be created.
* `text::AbstractString`: Text to write to the file.

Returns the number of bytes written to the file.
"""
function gen_file(file::AbstractString, text::AbstractString)
    mkpath(dirname(file))
    if !endswith(text , "\n")
        text *= "\n"
    end
    return write(file, text)
end

"""
    version_floor(v::VersionNumber=VERSION) -> String

Format the given Julia version.

# Keyword arguments
* `v::VersionNumber=VERSION`: Version to floor.

Returns "major.minor" for the most recent release version relative to v. For prereleases
with v.minor == v.patch == 0, returns "major.minor-".
"""
function version_floor(v::VersionNumber=VERSION)
    return if isempty(v.prerelease) || v.patch > 0
        "$(v.major).$(v.minor)"
    else
        "$(v.major).$(v.minor)-"
    end
end

"""
    substitute(template::AbstractString, view::Dict{String, Any}) -> String
    substitute(
        template::AbstractString,
        pkg_template::Template;
        view::Dict{String, Any}=Dict{String, Any}(),
    ) -> String

Replace placeholders in `template` with values in `view` via
[`Mustache`](https://github.com/jverzani/Mustache.jl). `template` is not modified.
If `pkg_template` is supplied, some default replacements are also performed.

For information on how to structure `template`, see "Defining Template Files" section in
[Custom Plugins](@ref).

**Note**: Conditionals in `template` without a corresponding key in `view` won't error,
but will simply be evaluated as false.
"""
substitute(template::AbstractString, view::Dict{String, Any}) = render(template, view)

function substitute(
    template::AbstractString,
    pkg_template::Template;
    view::Dict{String, Any}=Dict{String, Any}(),
)
    # Don't use version_floor here because we don't want the trailing '-' on prereleases.
    v = pkg_template.julia_version
    d = Dict{String, Any}(
        "USER" => pkg_template.user,
        "VERSION" => "$(v.major).$(v.minor)",
        "DOCUMENTER" => any(map(p -> isa(p, Documenter), values(pkg_template.plugins))),
        "CODECOV" => haskey(pkg_template.plugins, Codecov),
        "COVERALLS" => haskey(pkg_template.plugins, Coveralls),
    )
    # d["AFTER"] is true whenever something needs to occur in a CI "after_script".
    d["AFTER"] = d["DOCUMENTER"] || d["CODECOV"] || d["COVERALLS"]
    # d["COVERAGE"] is true whenever a coverage plugin is enabled.
    # TODO: This doesn't handle user-defined coverage plugins.
    # Maybe we need an abstract CoveragePlugin <: GenericPlugin?
    d["COVERAGE"] = d["CODECOV"] || d["COVERALLS"]
    return substitute(template, merge(d, view))
end

splitjl(pkg::AbstractString) = endswith(pkg, ".jl") ? pkg[1:end-3] : pkg

# Format a version in a way suitable for a Project.toml file.
function repr_version(v::VersionNumber)
    s = string(v.major)
    v.minor == 0 || (s *= ".$(v.minor)")
    v.patch == 0 || (s *= ".$(v.patch)")
    return repr(s)
end


text:
# Partial Pose3 constraints

"""
    $(TYPEDEF)

Partial prior belief on Z, Roll, and Pitch of a `Pose3`.
"""
mutable struct PriorPose3ZRP{T1,T2} <: IncrementalInference.FunctorSingleton where {T1 <: SamplableBelief, T2 <: SamplableBelief}
  z::T1
  rp::T2
  partial::Tuple{Int,Int,Int}
  PriorPose3ZRP{T1,T2}() where {T1, T2} = new{T1,T2}()
  PriorPose3ZRP{T1,T2}(z::T1,rp::T2) where {T1 <: IIF.SamplableBelief, T2 <: IIF.SamplableBelief} = new{T1,T2}(z, rp, (3,4,5))
end
PriorPose3ZRP(z::T1,rp::T2) where {T1 <: SamplableBelief, T2 <: SamplableBelief} = PriorPose3ZRP{T1,T2}(z, rp)
function getSample(pprz::PriorPose3ZRP, N::Int=1)
  return ([rand(pprz.z,N)[:]';rand(pprz.rp,N)], )
end

"""
    $SIGNATURES

Deprecated in favor of PriorPose3ZRP
"""
function convert(::Type{RoME.PriorPose3ZRP}, fact::Dict{String, Any})
    rp = fact["measurement"][1]
    z = fact["measurement"][2]
    rp = convert(_evalType(rp["distType"]), rp)
    z = convert(_evalType(z["distType"]), z)
    return PriorPose3ZRP(rp, z)
end

"""
Converter: PriorPose3ZRP::Dict{String, Any} -> PriorPose3ZRP

Deprecated in favor of PriorPose3ZRP
"""
function convert(::Type{Dict{String, Any}}, fact::RoME.PriorPose3ZRP)
    pf = Dict{String, Any}(
        "measurement" => [
            convert(Dict{String, Any}, fact.rp),
            convert(Dict{String, Any}, fact.z)
        ],
        "factorType" => "PriorPose3ZRP"
    )
    return pf
end

"""
    $(TYPEDEF)

Partial prior belief on Roll Pitch and Z of a `Pose3` variable.
"""
mutable struct PartialPriorRollPitchZ{T1,T2} <: IncrementalInference.FunctorSingleton where {T1 <: SamplableBelief, T2 <: SamplableBelief}
  rp::T1
  z::T2
  partial::Tuple{Int,Int,Int}
  PartialPriorRollPitchZ{T1,T2}() where {T1, T2} = new{T1,T2}()
  PartialPriorRollPitchZ{T1,T2}(rp::T1,z::T2) where {T1 <: IIF.SamplableBelief, T2 <: IIF.SamplableBelief} = new{T1,T2}(rp, z, (3,4,5))
end
PartialPriorRollPitchZ(rp::T1,z::T2) where {T1 <: SamplableBelief, T2 <: SamplableBelief} = PartialPriorRollPitchZ{T1,T2}(rp, z)
function getSample(pprz::PartialPriorRollPitchZ, N::Int=1)
  return ([rand(pprz.z,N)[:]';rand(pprz.rp,N)], )
end

"""
    $TYPEDEF

Serialization type of `PartialPriorRollPitchZ`.
"""
mutable struct PackedPartialPriorRollPitchZ <: IncrementalInference.PackedInferenceType
  rpdata::String
  zdata::String
  PackedPartialPriorRollPitchZ() = new()
  PackedPartialPriorRollPitchZ(x1::AS,x2::AS) where {AS <:AbstractString} = new(x1,x2)
end
function convert(::Type{PartialPriorRollPitchZ}, d::PackedPartialPriorRollPitchZ)
  # TODO: Change out for extractdistributionJson
  PartialPriorRollPitchZ( extractdistribution(d.rpdata), extractdistribution(d.zdata)  )
end
function convert(::Type{PackedPartialPriorRollPitchZ}, d::PartialPriorRollPitchZ)
  PackedPartialPriorRollPitchZ( string(d.rp), string(d.z) )
end

"""
Converter: PartialPriorRollPitchZ::Dict{String, Any} -> PartialPriorRollPitchZ
"""
function convert(::Type{RoME.PartialPriorRollPitchZ}, fact::Dict{String, Any})
    rp = fact["measurement"][1]
    z = fact["measurement"][2]
    rp = convert(_evalType(rp["distType"]), rp)
    z = convert(_evalType(z["distType"]), z)
    return PartialPriorRollPitchZ(rp, z)
end

"""
Converter: PartialPriorRollPitchZ::Dict{String, Any} -> PartialPriorRollPitchZ
"""
function convert(::Type{Dict{String, Any}}, fact::RoME.PartialPriorRollPitchZ)
    pf = Dict{String, Any}(
        "measurement" => [
            convert(Dict{String, Any}, fact.rp),
            convert(Dict{String, Any}, fact.z)
        ],
        "factorType" => "PartialPriorRollPitchZ"
    )
    return pf
end

function compare(a::Normal, b::Normal; tol::Float64=1e-10)
  TP = true
  TP = TP && abs(a.μ - b.μ) < tol
  TP = TP && abs(a.σ - b.σ) < tol
  TP
end

function compare(a::PartialPriorRollPitchZ, b::PartialPriorRollPitchZ; tol::Float64=1e-10)
  TP = true
  TP = TP && compare(a.rp, b.rp)
  TP = TP && compare(a.z, b.z)
  TP = TP && norm(collect(a.partial)-collect(b.partial)) < tol
  return TP
end

# Partial pairwise constraint between poses X,Y,Yaw
# ------------------------------------------------------------------------------

# TODO duplicate name Pose3Pose3XYYaw

"""
    $TYPEDEF

Partial factor between XY and Yaw of two Pose3 variables.

To be deprecated: use Pose3Pose3XYYaw instead.
"""
mutable struct PartialPose3XYYaw{T1,T2} <: FunctorPairwise where {T1 <: SamplableBelief, T2 <: SamplableBelief}
  xy::T1
  yaw::T2
  partial::Tuple{Int,Int,Int}
  PartialPose3XYYaw{T1,T2}() where {T1, T2} = new()
  PartialPose3XYYaw{T1,T2}(xy::T1, yaw::T2) where {T1 <: IIF.SamplableBelief, T2 <: IIF.SamplableBelief} = new(xy, yaw, (1,2,6))
end
PartialPose3XYYaw(xy::T1, yaw::T2) where {T1 <: IIF.SamplableBelief, T2 <: IIF.SamplableBelief} =  PartialPose3XYYaw{T1,T2}(xy, yaw)

function getSample(pxyy::PartialPose3XYYaw, N::Int=1)
  return ([rand(pxyy.xy,N);rand(pxyy.yaw,N)[:]'], )
end
function (pxyy::PartialPose3XYYaw)(res::Array{Float64},
            userdata,
            idx::Int,
            meas::Tuple{Array{Float64,2}},
            wXi::Array{Float64,2},
            wXj::Array{Float64,2}  )
  #
  wXjhat = SE2(wXi[[1;2;6],idx]) * SE2(meas[1][1:3,idx])
  jXjhat = SE2(wXj[[1;2;6],idx]) \ wXjhat
  se2vee!(res, jXjhat)
  nothing
end

"""
    $TYPEDEF

Partial factor between XY and Yaw of two Pose3 variables.
"""
mutable struct Pose3Pose3XYYaw{T1,T2} <: FunctorPairwise where {T1 <: SamplableBelief, T2 <: SamplableBelief}
  xy::T1
  yaw::T2
  partial::Tuple{Int,Int,Int}
  Pose3Pose3XYYaw{T1,T2}() where {T1, T2} = new()
  Pose3Pose3XYYaw{T1,T2}(xy::T1, yaw::T2) where {T1 <: IIF.SamplableBelief, T2 <: IIF.SamplableBelief} = new(xy, yaw, (1,2,6))
end
Pose3Pose3XYYaw(xy::T1, yaw::T2) where {T1 <: IIF.SamplableBelief, T2 <: IIF.SamplableBelief} =  Pose3Pose3XYYaw{T1,T2}(xy, yaw)

function getSample(pxyy::Pose3Pose3XYYaw, N::Int=1)
  return ([rand(pxyy.xy,N);rand(pxyy.yaw,N)[:]'], )
end
function (pxyy::Pose3Pose3XYYaw)(res::Array{Float64},
            userdata,
            idx::Int,
            meas::Tuple{Array{Float64,2}},
            wXi::Array{Float64,2},
            wXj::Array{Float64,2}  )
  #
  wXjhat = SE2(wXi[[1;2;6],idx]) * SE2(meas[1][1:3,idx])
  jXjhat = SE2(wXj[[1;2;6],idx]) \ wXjhat
  se2vee!(res, jXjhat)
  nothing
end

"""
    $TYPEDEF

Serialization type of PartialPose3XYYaw.
"""
mutable struct PackedPartialPose3XYYaw <: IncrementalInference.PackedInferenceType
  xydata::String
  yawdata::String
  PackedPartialPose3XYYaw() = new()
  PackedPartialPose3XYYaw(xy::String, yaw::String) = new(xy, yaw)
end
function convert(::Type{PartialPose3XYYaw}, d::PackedPartialPose3XYYaw)
  return PartialPose3XYYaw( extractdistribution(d.xydata), extractdistribution(d.yawdata) )
end
function convert(::Type{PackedPartialPose3XYYaw}, d::PartialPose3XYYaw)
  return PackedPartialPose3XYYaw( string(d.xy), string(d.yaw) )
end

"""
    $SIGNATURES

Converter: Dict{String, Any} -> PartialPose3XYYaw
"""
function convert(::Type{Dict{String, Any}}, fact::RoME.PartialPose3XYYaw)
    pf = Dict{String, Any}(
        "measurement" => [
            convert(Dict{String, Any}, fact.xy),
            convert(Dict{String, Any}, fact.yaw)
        ],
        "factorType" => "PartialPose3XYYaw"
    )
    return pf
end

"""
    $SIGNATURES

Converter: PartialPose3XYYaw -> Dict{String, Any}
"""
function convert(::Type{RoME.PartialPose3XYYaw}, fact::Dict{String, Any})
    xy = fact["measurement"][1]
    yaw = fact["measurement"][2]
    xy = convert(_evalType(xy["distType"]), xy)
    yaw = convert(_evalType(yaw["distType"]), yaw)
    return PartialPriorRollPitchZ(xy, yaw)
end

function compare(a::PartialPose3XYYaw, b::PartialPose3XYYaw; tol::Float64=1e-10)
  TP = true
  TP = TP && compare(a.xy, b.xy)
  TP = TP && compare(a.yaw, b.yaw)
  TP = TP && norm(collect(a.partial)-collect(b.partial)) < tol
  return TP
end

"""
OBSOLETE: see https://github.com/JuliaRobotics/IncrementalInference.jl/issues/237
"""
mutable struct PartialPose3XYYawNH <: IncrementalInference.FunctorPairwiseNH
  xyy::Distributions.MvNormal
  partial::Tuple{Int, Int, Int}
  nullhypothesis::Distributions.Categorical
  PartialPose3XYYawNH() = new()
  PartialPose3XYYawNH(xyy::MvNormal, vh::Vector{Float64}) = new(xyy, (1,2,6),  Distributions.Categorical(vh))
end
function getSample(pxyy::PartialPose3XYYawNH, N::Int=1)
  return (rand(pxyy.xyy,N), )
end
function (pxyy::PartialPose3XYYawNH)(res::Array{Float64},
            userdata,
            idx::Int,
            meas::Tuple{Array{Float64,2}},
            wXi::Array{Float64,2},
            wXj::Array{Float64,2}  )
  #
  wXjhat = SE2(wXi[[1;2;6],idx])*SE2(meas[1][:,idx]) #*SE2(pp2.Zij[:,1])*SE2(meas[1][:,idx])
  jXjhat = SE2(wXj[[1;2;6],idx]) \ wXjhat
  se2vee!(res, jXjhat)
  nothing
end

mutable struct PackedPartialPose3XYYawNH <: IncrementalInference.PackedInferenceType
  vecZij::Array{Float64,1} # 3translations, 3rotation
  vecCov::Array{Float64,1}
  nullhypothesis::Vector{Float64}
  PackedPartialPose3XYYawNH() = new()
  PackedPartialPose3XYYawNH(x1::Vector{Float64}, x2::Array{Float64}, x3::Vector{Float64}) = new(x1, x2[:], x3)
end
function convert(::Type{PartialPose3XYYawNH}, d::PackedPartialPose3XYYawNH)
  return PartialPose3XYYawNH( Distributions.MvNormal(d.vecZij,
               reshapeVec2Mat(d.vecCov, 3)), d.nullhypothesis  )
end
function convert(::Type{PackedPartialPose3XYYawNH}, d::PartialPose3XYYawNH)
  return PackedPartialPose3XYYawNH(d.xyy.μ, d.xyy.Σ.mat, d.nullhypothesis.p )
end

function compare(a::PartialPose3XYYawNH, b::PartialPose3XYYawNH; tol::Float64=1e-10)
  TP = true
  TP = TP && norm(a.xyy.μ-b.xyy.μ) < tol
  TP = TP && norm(a.xyy.Σ.mat[:]-b.xyy.Σ.mat[:]) < tol
  TP = TP && norm(collect(a.partial)-collect(b.partial)) < tol
  TP = TP && norm(a.nullhypothesis.p-b.nullhypothesis.p) < tol
  return TP
end

#


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


text:
"""
myLS(y, X)

A simple least squares implementation.
"""
struct myLS <: myEstimator
    β::Array{Float64} # coefficient
    y::Array{Float64} # response
    X::Array{Float64} # features

    # Define constructor function
    function myLS(y::Array{Float64}, X::Array{Float64})

        # Calculate LS
        β = (X' * X) \ (X' * y)

        # Organize and return output
        new(β, y, X)
    end #MYLS
end #MYLS

# Methods ======================================================================
"""
inference(fit::myLS; heteroskedastic, print_df)

A method to calculate standard errors of a myLS object.

To do: clustered se
"""
function inference(fit::myLS; heteroskedastic::Bool=false, print_df::Bool=true)
    # Obtain data parameters
    N = length(fit.y)
    K = size(fit.X, 2)

    # Calculate the covariance matrix
    u = fit.y - predict(fit) # residuals
    XX_inv = inv(fit.X' * fit.X)
    if !heteroskedastic
        # homoskedastic se
        covar = sum(u.^2) * XX_inv
	    covar = covar .* (1 / (N - K)) # dof adjustment
    else
        # heteroskedastic se
        covar = XX_inv * ((fit.X .* (u.^2))' * fit.X) * XX_inv
		covar = covar .* (N / (N - K)) # dof adjustment
    end

    # Get standard errors, t-statistics and p-values
    se = sqrt.(covar[diagind(covar)])
    t_stat = fit.β ./ se
    p_val = 2 * cdf.(Normal(), -abs.(t_stat))

    # Print estimates
    if print_df
        out_df = DataFrame(hcat(fit.β, se, t_stat, p_val), :auto)
        rename!(out_df, ["coef", "se", "t-stat", "p-val"])
        display(out_df)
    end

    # Organize and return output
    output = (β = fit.β, se = se, t = t_stat, p = p_val)
    return output
end #INFERENCE.MYLS


text:
module Kwon

using Parameters
using TightBinding: SKHamiltonian, hop, TBModel, NullPotential, ORTHOGONAL
using JuLIP.Potentials: fcut, fcut_d, SplineCutoff,
   PairPotential, ZeroPairPotential, EAM, @analytic

import JuLIP: cutoff
import TightBinding: hop, onsite!, onsite_grad!

"""
`KwonHamiltonian <: SKHamiltonian{ORTHOGONAL, 4}`

Hamiltonian for an orthogonal sp TB model of Si developed by Kwon et al [1].

This implementation deviates  from [1] in how the cut-off is applied:
instead of "patching" a cubic spline between r1 and rcut, we simply multiply
with a quintic spline on the interval [0.5 (rcut + r0), rcut].

[1] I. Kwon, R. Biswas, C. Z. Wang, K. M. Ho and C. M. Soukoulis.
Transferable tight-binding models for silicon.
Phys Rev B 49 (11), 1994.
"""
@with_kw mutable struct KwonHamiltonian <: SKHamiltonian{ORTHOGONAL, 4}
   r0::Float64 = 2.360352   # Å
   Es::Float64 = -5.25       # eV
   Ep::Float64 = 1.2        # eV
   E0::Float64 = 8.7393204  # eV   8.7393204  is the original value; but it is just a constant
   # ------------------------------- Electronic Parameters
   # α   1     2    3     4
   #    ssσ   spσ  ppσ   ppπ
   hr0::NTuple{4, Float64} = (-2.038, 1.745, 2.75, -1.075)
   nc::NTuple{4, Float64} = (9.5, 8.5, 7.5, 7.5)
   rc::NTuple{4, Float64} = (3.4, 3.55, 3.7, 3.7)
   # ------------------------------- 2-body parameters
   m::Float64 = 6.8755
   mc::Float64 = 13.017
   dc::Float64 = 3.66995     # Å
   C::NTuple{4, Float64} = (2.1604385, -0.1384393, 5.8398423e-3, -8.0263577e-5)
   # --------------------------------
   r1::Float64 = 3.260176     # Å  (start of cut-off) >>> different r1 from Kwon paper
   rcut::Float64 = 4.16       # Å  (end of cut-off)
end

cutoff(H::KwonHamiltonian) = H.rcut

kwon_hop(H::KwonHamiltonian, r::Real, α) = ( H.hr0[α] * (H.r0 / r)^2 *
               exp( - 2 * (r/H.rc[α])^H.nc[α] + 2 * (H.r0/H.rc[α])^H.nc[α] ) )

hop(H::KwonHamiltonian, r::Real, α) = kwon_hop(H, r, α) * fcut(r, H.r1, H.rcut)

function onsite!(H::KwonHamiltonian, _r, _R, H_nn)
   fill!(H_nn, 0.0)
   H_nn[1,1] = H.Es
   H_nn[2,2] = H_nn[3,3] = H_nn[4,4] = H.Ep
   return H_nn
end

onsite_grad!(H::KwonHamiltonian, _r, _R, dH_nn) = fill!(dH_nn, 0.0)

# ============ Repulsive potential

# embedding function C1 x + C2 x^2 + C3 x^3 + C3 x^4
KwonEmbedding(H::KwonHamiltonian) =
   let E0=H.E0, C1=H.C[1], C2=H.C[2], C3=H.C[3], C4=H.C[4]
      @analytic r -> E0 + r * (C1 + r * (C2 + r * (C3 + r * C4)))
   end

# electron density: (r0/r)^m * exp( - m * (r/dc)^mc + m * (r0/dc)^mc )
KwonElDensity(H::KwonHamiltonian) =
   let r0 = H.r0, m = H.m, dc = H.dc, mc = H.mc
      @analytic r -> (r0/r)^m * exp( - m * (r/dc)^mc + m * (r0/dc)^mc )
   end * SplineCutoff(H.r1, H.rcut)

KwonEAM(H::KwonHamiltonian) = EAM(ZeroPairPotential(), KwonElDensity(H), KwonEmbedding(H))

function KwonTBModel(; potential = NullPotential(), bzquad = NullBZQ())
   H = KwonHamiltonian()
   return TBModel(H, KwonEAM(H), potential, bzquad, 1e-6)
end

end


text:
Coord = Tuple{Float64,Float64,Float64}

struct Facet
    normal::Coord
    loop::Vector{Coord}
end
function Facet(normal::Coord)
    Facet(normal, Vector{Coord}())
end
function Facet()
    Facet(Coord([0., 0., 0.]))
end
function Facet(coord::Vector{Float64})
    if length(coord) == 3
        Facet(Coord(coord))
    else
        Facet()
    end
end

function add_vertex!(facet::Facet, vertex::Coord)
    push!(facet.loop, vertex)
end
function add_vertex!(facet::Facet, vertex::Vector{Float64})
    if length(vertex) == 3
        add_vertex!(facet, Coord(vertex))
    end
end

struct Design
    name::String
    facets::Vector{Facet}
end
function Design(name::String = "OBJECT")
    Design(name, Vector{Facet}())
end

function add_facet!(design::Design, facet::Facet)
    push!(design.facets, facet)
end


text:
mutable struct Atoms
    Natoms::Int64
    Nspecies::Int64
    positions::Array{Float64,2}
    atm2species::Array{Int64,1}
    atsymbs::Array{String,1}
    SpeciesSymbols::Array{String,1}  # unique symbols
end

# Overload println for Atoms
import Base.println
function println( a::Atoms )
    @printf("\n")
    @printf("                                     -----\n")
    @printf("                                     Atoms\n")
    @printf("                                     -----\n")
    @printf("\n")
    @printf("Natoms   = %5d\n", a.Natoms)
    @printf("Nspecies = %5d\n", a.Nspecies)
    @printf("\n")
    for isp = 1:a.Nspecies
        @printf("Species %3d: %s\n", isp, a.SpeciesSymbols[isp])
    end
    @printf("\n")
    @printf("Coordinates in bohr:\n")
    @printf("\n")
    for ia = 1:a.Natoms
        @printf("%5s %18.10f %18.10f %18.10f\n", a.atsymbs[ia],
                a.positions[1,ia], a.positions[2,ia], a.positions[3,ia])
    end
end

# dummy atoms, contains only one atom
function Atoms()
    Natoms = 1
    Nspecies = 1
    positions = zeros(3,Natoms)
    atm2species = [1]
    atsymbs = ["X"]
    SpeciesSymbols = ["X"]  # unique symbols
    return Atoms( Natoms, Nspecies, positions, atm2species, atsymbs, SpeciesSymbols )
end

"""
Initialize from xyz file
"""
function init_atoms_xyz(filexyz; in_bohr=false, verbose=false)
    f = open(filexyz, "r")
    l = readline(f)
    Natoms = parse(Int64, l)
    positions = zeros(3,Natoms)
    atsymbs = Array{String}(undef,Natoms)
    l = readline(f)
    for ia = 1:Natoms
        ll = split( readline(f) )
        atsymbs[ia] = ll[1]
        positions[1,ia] = parse( Float64, ll[2] )
        positions[2,ia] = parse( Float64, ll[3] )
        positions[3,ia] = parse( Float64, ll[4] )
    end
    close(f)

    # convert from angstrom to bohr
    if !in_bohr
        if verbose
            println("Coordinate in xyz file is assumed to be given in angstrom")
            println("It will be converted to bohr")
        end
        positions[:,:] = positions[:,:] * ANG2BOHR
    else
        if verbose
            println("Coordinate in xyz file is assumed to be given in bohr")
        end
    end

    # Determine number of species
    Nspecies = 0
    for ia = 1:Natoms
        k2 = 0
        for k1 = 1:ia-1
            if atsymbs[k1] == atsymbs[ia]
                k2 = 1
            end
        end
        # find different
        if k2 == 0
            Nspecies = Nspecies + 1
        end
    end

    SpeciesSymbols = Array{String}(undef,Nspecies)

    idx1 = 0
    for ia = 1:Natoms
        k2 = 0
        for k1 = 1:ia-1
            if atsymbs[k1] == atsymbs[ia]
                k2 = 1
            end
        end
        # Found different species
        if k2==0
            idx1 = idx1 + 1
            SpeciesSymbols[idx1] = atsymbs[ia]
        end
    end

    # Mapping of atoms to species index
    atm2species = Array{Int64}(undef,Natoms)
    for ia = 1:Natoms
        for isp = 1:Nspecies
            if atsymbs[ia] == SpeciesSymbols[isp]
                atm2species[ia] = isp
            end
        end
    end

    return Atoms(Natoms, Nspecies, positions, atm2species, atsymbs, SpeciesSymbols )

    #ALLOCATE( AtomicMasses(Nspecies) )
    #AtomicMasses(:) = 0.d0 ! FIXME: Use internal database to set this

end

const ZATOMS = Dict(
    "H"  => 1,
    "He" => 2,
    "Li" => 3,
    "Be" => 4,
    "B"  => 5,
    "C"  => 6,
    "N"  => 7,
    "O"  => 8,
    "F"  => 9,
    "Ne" => 10,
    "Na" => 11,
    "Mg" => 12,
    "Al" => 13,
    "Si" => 14,
    "P"  => 15,
    "S"  => 16,
    "Cl" => 17,
    "Ar" => 18,
    "K"  => 19,
    "Ca" => 20,
    "Sc" => 21,
    "Ti" => 22,
    "V"  => 23,
    "Cr" => 24,
    "Mn" => 25,
    "Fe" => 26,
    "Co" => 27,
    "Ni" => 28,
    "Cu" => 29,
    "Zn" => 30,
    "Ga" => 31,
    "Ge" => 32,
    "As" => 33,
    "Se" => 34,
    "Br" => 34,
    "Kr" => 36,
    "Rb" => 37,
    "Sr" => 38,
    "Y"  => 39,
    "Zr" => 40,
    "Nb" => 41,
    "Mo" => 42,
    "Tc" => 43,
    "Ru" => 44,
    "Rh" => 45,
    "Pd" => 46,
    "Ag" => 47,
    "Cd" => 48,
    "In" => 49,
    "Sn" => 50,
    "Sb" => 51,
    "Te" => 52,
    "I"  => 53,
    "Xe" => 54, )

function get_Zatoms( atoms::Atoms )
    Nspecies = atoms.Nspecies
    SpeciesSymbols = atoms.SpeciesSymbols
    Zatoms = zeros(Float64,Nspecies)
    for isp = 1:Nspecies
        Zatoms[isp] = ZATOMS[SpeciesSymbols[isp]]
    end
    return Zatoms
end

function get_Nelectrons( atoms:: Atoms )
    Zatoms = get_Zatoms(atoms)
    Natoms = atoms.Natoms
    atm2species = atoms.atm2species
    Nelectrons = 0.0
    for ia = 1:Natoms
        isp = atm2species[ia]
        Nelectrons = Nelectrons + Zatoms[isp]
    end
    return Nelectrons
end


text:
########################
### Filter functions ###----------------------------------------------------------------------------
########################

#' Basic Eddy Covariance Data Filtering
#'
#' Filters time series of EC data for high-quality values and specified
#'              meteorological conditions.
#'
#' - data            DataFrame or matrix containing all required input variables in
#'                        half-hourly or hourly resolution. Including year, month, day information
#' - quality_control Should quality control be applied? Defaults to `true`.
#' - filter_growseas Should data be filtered for growing season? Defaults to `false`.
#' - filter_precip   Should precipitation filtering be applied? Defaults to `false`.
#' - filter_vars     Additional variables to be filtered. Vector of type character.
#' - filter_vals_min Minimum values of the variables to be filtered. Numeric vector of
#'                        the same length than `filter_vars`. Set to `missing` to be ignored.
#' - filter_vals_max Maximum values of the variables to be filtered. Numeric vector of
#'                        the same length than `filter_vars`. Set to `missing` to be ignored.
#' - NA_as_invalid   If `true` (the default) missing data are filtered out (applies to all variables).
#' - vars_qc         Character vector indicating the variables for which quality filter should
#'                        be applied. Ignored if `quality_control = false`.
#' - quality_ext     The extension to the variables' names that marks them as
#'                        quality control variables. Ignored if `quality_control = false`.
#' - good_quality    Which values indicate good quality (i.e. not to be filtered)
#'                        in the quality control (qc) variables? Ignored if `quality_control = false`.
#' - missing_qc_as_bad If quality control variable is `missing`, should the corresponding data point be
#'                          treated as bad quality? Defaults to `true`. Ignored if `quality_control = false`.
#' - precip          Precipitation (mm time-1)
#' - GPP             Gross primary productivity (umol m-2 s-1); Ignored if `filter_growseas = false`.
#' - doy             Day of year; Ignored if `filter_growseas = false`.
#' - year            Year; Ignored if `filter_growseas = false`.
#' - tGPP            GPP threshold (fraction of 95th percentile of the GPP time series).
#'                        Must be between 0 and 1. Ignored if `filter_growseas` is `false`.
#' - ws              Window size used for GPP time series smoothing.
#'                        Ignored if `filter_growseas = false`.
#' - min_int         Minimum time interval in days for a given state of growing season.
#'                        Ignored if `filter_growseas = false`.
#' - tprecip         Precipitation threshold used to identify a precipitation event (mm).
#'                        Ignored if `filter_precip = false`.
#' - precip_hours    Number of hours removed following a precipitation event (h).
#'                        Ignored if `filter_precip = false`.
#' - records_per_hour Number of observations per hour. I_e. 2 for half-hourly data.
#' - filtered_data_to_NA Logical. If `true` (the default), all variables in the input
#'                              DataFrame/matrix are set to `missing` for the time step where ANY of the
#'                              `filter_vars` were beyond their acceptable range (as
#'                              determined by `filter_vals_min` and `filter_vals_max`).
#'                              If `false`, values are not filtered, and an additional column 'valid'
#'                              is added to the DataFrame/matrix, indicating if any value of a row
#'                              did (1) or did not fulfill the filter criteria (0).
#' - constants frac2percent - conversion between fraction and percent
#'
#' # Details
 This routine consists of two parts:
#'
#'          1) Quality control: All variables included in `vars_qc` are filtered for
#'             good quality data. For these variables a corresponding quality variable with
#'             the same name as the variable plus the extension as specified in `quality_ext`
#'             must be provided. For time steps where the value of the quality indicator is not included
#'             in the argument `good_quality`, i.e. the quality is not considered as 'good',
#'             its value is set to `missing`.
#'
#'          2) Meteorological filtering. Under certain conditions (e.g. low ustar), the assumptions
#'             of the EC method are not fulfilled. Further, some data analysis require certain meteorological
#'             conditions, such as periods without rainfall, or active vegetation (growing season, daytime).
#'             The filter applied in this second step serves to exclude time periods that do not fulfill the criteria
#'             specified in the arguments. More specifically, time periods where one of the variables is higher
#'             or lower than the specified thresholds (`filter_vals_min` and `filter_vals_max`)
#'             are set to `missing` for all variables. If a threshold is set to `missing`, it will be ignored.
#'
#' # Value
 If `filtered_data_to_NA = true` (default), the input DataFrame/matrix with
#'         observations which did not fulfill the filter criteria set to `missing`.
#'         If `filtered_data_to_NA = false`, the input DataFrame/matrix with an additional
#'         column "valid", which indicates whether all the data of a time step fulfill the
#'         filtering criteria (1) or not (0).
#'
#' # Note
#' The thresholds set with `filter_vals_min` and `filter_vals_max` filter all data
#'       that are smaller than ("<"), or greater than (">") the specified thresholds. That means
#'       if a variable has exactly the same value as the threshold, it will not be filtered. Likewise,
#'       `tprecip` filters all data that are greater than `tprecip`.
#'
#'       Variables considered of bad quality (as specified by the corresponding quality control variables)
#'       will be set to `missing` by this routine. Data that do not fulfill the filtering criteria are set to
#'       `missing` if `filtered_data_to_NA = true`. Note that with this option *all* variables of the same
#'       time step are set to `missing`. Alternatively, if `filtered_data_to_NA = false` data are not set to `missing`,
#'       and a new column "valid" is added to the DataFrame/matrix, indicating if any value of a row
#'       did (1) or did not fulfill the filter criteria (0).
#'
#'
#' ```@example; output = false
#' ```
#' # Example of data filtering; data are for a month within the growing season,
#' # hence growing season is not filtered.
#' # If filtered_data_to_NA=true, all values of a row are set to missing if one filter
#' # variable is beyond its bounds.
#' DE_Tha_Jun_2014_2 = filter_data(DE_Tha_Jun_2014,quality_control=false,
#'                                  vars_qc=c("Tair","precip","H","LE"),
#'                                  filter_growseas=false,filter_precip=true,
#'                                  filter_vars=c("Tair","PPFD","ustar"),
#'                                  filter_vals_min=c(5,200,0.2),
#'                                  filter_vals_max=c(missing,missing,missing),NA_as_invalid=true,
#'                                  quality_ext="_qc",good_quality=c(0,1),
#'                                  missing_qc_as_bad=true,GPP="GPP",doy="doy",
#'                                  year="year",tGPP=0.5,ws=15,min_int=5,precip="precip",
#'                                  tprecip=0.1,precip_hours=24,records_per_hour=2,
#'                                  filtered_data_to_NA=true)
#'
#'  ## same, but with filtered_data_to_NA=false
#'  DE_Tha_Jun_2014_3 = filter_data(DE_Tha_Jun_2014,quality_control=false,
#'                                  vars_qc=c("Tair","precip","H","LE"),
#'                                  filter_growseas=false,filter_precip=true,
#'                                  filter_vars=c("Tair","PPFD","ustar"),
#'                                  filter_vals_min=c(5,200,0.2),
#'                                  filter_vals_max=c(missing,missing,missing),NA_as_invalid=true,
#'                                  quality_ext="_qc",good_quality=c(0,1),
#'                                  missing_qc_as_bad=true,GPP="GPP",doy="doy",
#'                                  year="year",tGPP=0.5,ws=15,min_int=5,precip="precip",
#'                                  tprecip=0.1,precip_hours=24,records_per_hour=2,
#'                                  filtered_data_to_NA=false)
#'
#'  # note the additional column 'valid' in DE_Tha_Jun_2014_3.
#'  # To remove time steps marked as filtered out (i.e. 0 values in column 'valid'):
#'  DE_Tha_Jun_2014_3[DE_Tha_Jun_2014_3["valid"] == 0,] = missing
#'
#'
#' @importFrom stats aggregate
#' @export
function filter_data(data,quality_control=true,filter_growseas=false,
                        filter_precip=false,filter_vars=nothing,
                        filter_vals_min,filter_vals_max,NA_as_invalid=true,
                        vars_qc=nothing,quality_ext="_qc",good_quality=c(0,1),
                        missing_qc_as_bad=true,GPP="GPP",doy="doy",
                        year="year",tGPP=0.5,ws=15,min_int=5,precip="precip",
                        tprecip=0.01,precip_hours=24,records_per_hour=2,
                        filtered_data_to_NA=true,constants=BigleafConstants())

  ### I) Quality control filter
  if (quality_control)

    if (isnothing(vars_qc))
      stop("quality_control (qc) is true, but no qc variables are provided!")
end

    if (any(!vars_qc in colnames(data)))

      missing_vars = vars_qc[which(!vars_qc in colnames(data))]
      stop(paste("Variable ",missing_vars," is included in 'vars_qc', but does not exist in the input data!"))

end

    vars_qc_qc = paste0(vars_qc,quality_ext)
    if (any(!vars_qc_qc in colnames(data)))

      missing_vars_qc = vars_qc_qc[which(!vars_qc_qc in colnames(data))]
      missing_vars2   = substr(missing_vars_qc,1,nchar(missing_vars_qc) - nchar(quality_ext))
      stop(paste("Quality control for variable ",missing_vars2,"(",missing_vars_qc,") does not exist in the input data!"))
end

    ## data quality
    cat("Quality control:",fill=true)
    for (var in vars_qc)
      var_qc = paste0(var,quality_ext)
      check_input(data,var)
      check_input(data,var_qc)

      if (missing_qc_as_bad)
        data[get(paste0(var,quality_ext)) > max(good_quality) | ismissing(get(paste0(var,quality_ext))),var] = missing   # exclude bad quality data or those where qc flag is not available
        qc_invalid      = sum(get(paste0(var,quality_ext)) > max(good_quality) | ismissing(get(paste0(var,quality_ext)))) # count & report
else { # same, but consider missing quality flag variables as good
        data[get(paste0(var,quality_ext)) > max(good_quality) & !ismissing(get(paste0(var,quality_ext))),var] = missing
        qc_invalid      = sum(get(paste0(var,quality_ext)) > max(good_quality) & !ismissing(get(paste0(var,quality_ext))))
end

      qc_invalid_perc = round((qc_invalid/nrow(data))*constants.frac2percent,2)

      cat(var,": ",qc_invalid," data points (",qc_invalid_perc,"%) set to missing",fill=true,sep="")
end
end

  #### II) Data filter
  valid = rep(1L,nrow(data))

  # 1) GPP
  growseas_invalid = numeric()
  if(filter_growseas)
    check_input(data,doy,year,GPP)
    date             = strptime(paste0(year,"-",doy),format="%Y-%j")
    GPP_daily        = aggregate(GPP,by=list(strftime(date)),mean,na_rm=true)[,2]
    growing_season   = filter_growingseason(GPP_daily,tGPP=tGPP,ws=ws,min_int=min_int)
    growseas_invalid = which(sapply(growing_season,rep,48) == 0)
end

  # 2) precipitation
  precip_invalid = numeric()
  if (filter_precip)
    check_input(data,precip)
    if (NA_as_invalid)
      precip_events = which(precip > tprecip | ismissing(precip))
else
      precip_events = which(precip > tprecip)
end
    precip_invalid = unique(as_numeric(unlist(sapply(precip_events, function(x) x:(min(x+precip_hours*records_per_hour,nrow(data),na_rm=true))))))
end

  # 3) all other filter variables (as defined in filter_vars)
  invalids = list(growseas_invalid,precip_invalid)

  if (!isnothing(filter_vars))
    for (var in filter_vars)
      v  = which(filter_vars == var)
      vf = v + 2
      check_input(data,var)
      if (NA_as_invalid)
        invalids[[vf]] = which(get(var) < filter_vals_min[v] | get(var) > filter_vals_max[v] | ismissing(get(var)))
else
        invalids[[vf]] = which(get(var) < filter_vals_min[v] | get(var) > filter_vals_max[v] & !ismissing(get(var)))
end
end
end

  # 4) calculate number and percentage of filtered values
  invalids_perc = sapply(invalids, function(x) round((length(x)/nrow(data))*constants.frac2percent,2))

  additional_invalids = sapply(2:length(invalids), function(x)
    length(setdiff(invalids[[x]],unique(unlist(invalids[1:(x-1)])))))

  additional_invalids_perc = round(additional_invalids/nrow(data)*constants.frac2percent,2)

  # 5) write to output
  if (filter_growseas | filter_precip | length(filter_vars) > 0)

    var_names = c("growing season","precipitation",filter_vars)

    if (quality_control)
      cat("-------------------------------------------------------------------",fill=true)
end

    cat("Data filtering:",fill=true)

    cat(length(growseas_invalid)," data points (",invalids_perc[1],"%) excluded by growing season filter",fill=true,sep="")

    invisible(sapply(c(1:(length(invalids)-1)), function(x) cat(additional_invalids[x]," additional data points (",
                                                                additional_invalids_perc[x],"%) excluded by ",var_names[x+1],
                                                                " filter (",length(unlist(invalids[x+1]))," data points = ",
                                                                invalids_perc[x+1]," % in total)",fill=true,sep="")))

    invalid        = unique(unlist(invalids))
    valid[invalid] = 0

    excl_perc = round((length(invalid)/nrow(data))*constants.frac2percent,2)

    cat(length(invalid)," data points (",excl_perc,"%) excluded in total",fill=true,sep="")
    cat(nrow(data) - length(invalid)," valid data points (",constants.frac2percent-excl_perc,"%) remaining.",fill=true,sep="")

    # 6) return input data frame with filtered time steps set to missing or an additional 'valid' column
    if (filtered_data_to_NA)
      data_filtered = data
      data_filtered[valid < 1,] = missing
else
      data_filtered = DataFrame(data,valid)
end

else

    data_filtered = data

end

  return(data_filtered)
end

#' GPP-based Growing Season Filter
#'
#' Filters annual time series for growing season based on smoothed daily GPP data.
#'
#' - GPPd    daily GPP (any unit)
#' - tGPP    GPP threshold (fraction of 95th percentile of the GPP time series).
#'                Takes values between 0 and 1.
#' - ws      window size used for GPP time series smoothing
#' - min_int minimum time interval in days for a given state of growing season
#'
#' # Details
 The basic idea behind the growing season filter is that vegetation is
#'          considered to be active when its carbon uptake (GPP) is above a specified
#'          threshold, which is defined relative to the peak GPP (95th percentile)
#'          observed in the year.
#'          The GPP-threshold is calculated as:
#'
#'          ``GPP_threshold = quantile(GPPd,0.95)*tGPP``
#'
#'          GPPd time series are smoothed with a moving average to avoid fluctuations
#'          in the delineation of the growing season. The window size defaults to 15
#'          days, but depending on the ecosystem, other values can be appropriate.
#'
#'          The argument `min_int` serves to avoid short fluctuations in the
#'          status growing season vs. no growing season by defining a minimum length
#'          of the status. If a time interval shorter than `min_int` is labeled
#'          as growing season or non-growing season, it is changed to the status of
#'          the neighboring values.
#'
#' # Value
 a vector of type integer of the same length as the input GPPd in which 0 indicate
#'         no growing season (dormant season) and 1 indicate growing season.
#'
#' @importFrom stats quantile filter
#' @export
function filter_growingseason(GPPd,tGPP,ws=15,min_int=5)

  if(sum(ismissing(GPPd)) < 0.5*length(GPPd))

    growseas      = rep(1,length(GPPd))
    GPP_threshold = quantile(GPPd,probs=0.95,na_rm=true)*tGPP

    ## smooth GPP
    GPPd_smoothed = filter(GPPd,method="convolution",filter=rep(1/ws,ws))

    ## set values at the beginning and end of the time series to the mean of the original values
    wsd = floor(ws/2)
    GPPd_smoothed[1:wsd] = mean(GPPd[1:(2*wsd)],na_rm=true)
    GPPd_smoothed[(length(GPPd)-(wsd-1)):length(GPPd)] = mean(GPPd[(length(GPPd)-(2*wsd-1)):length(GPPd)],na_rm=true)

    # check for occurence of missing values and set them to mean of the values surrounding them
    missing = which(ismissing(GPPd_smoothed))
    if (length(missing) > 0)
      if (length(missing) > 10){@warn"Attention, there is a gap in 'GPPd' of length n = ",length(missing))}
      replace_val = mean(GPPd_smoothed[max(1,missing[1] - 4):min((missing[length(missing)] + 4),length(GPPd_smoothed))],na_rm=true)
      GPPd_smoothed[missing] = replace_val
end

    # filter daily GPP
    growseas[GPPd_smoothed < GPP_threshold] = 0

    ## change short intervals to the surrounding values to avoid 'wrong' fluctuations
    intervals = rle(growseas)
    short_int = which(intervals$lengths <= min_int)

    if (length(short_int) > 0)
      start = numeric()
      end   = numeric()

      for (i in 1:length(short_int))

        start[i] = sum(intervals$lengths[1:short_int[i]-1]) + 1
        end[i]   = start[i]+intervals$lengths[short_int[i]] - 1

        val = unique(growseas[start[i]:end[i]])

        if (val == 0 & growseas[start[i]-1] == 1)
          growseas[start[i]:end[i]] = 1
elseif (val == 1 & growseas[start[i]-1] == 0)
          growseas[start[i]:end[i]] = 0
end
end
end

    growseas = as_integer(growseas)

else

    @warn"number of available GPPd data is less than half the total number of days per year. Filter is not applied!")
    growseas = as_integer(rep(1,length(GPPd)))

end

  return(growseas)
end


text:
using Plots
using Polynomials
using DataFrames, DataFramesMeta
using Latexify
texpyplot()

using Random; Random.seed!(1234)
x = rand(-10:20,10)
y = 2 * x.^2 .+ rand(-80:80,length(x))
#x = [-3, 4, 5, -1, 2, -1, 0, -3, 3, -2]
#y = [-6, 6, 10, -1, 5, -2, -2, -7, 4, -2]

x[end], y[end]

df = @linq DataFrame(x = x[1:end-1], y = y[1:end-1]) |>
    by(:x, y = mean(:y))

g = polyfit(df.x, df.y)

xs = range(minimum(x), maximum(x), length = 200)

plot(xs, g.(xs), legend = false, linewidth = 2, grid = true)
plot!(x[1:end-1],y[1:end-1], seriestype = :scatter, markercolor = :lightgreen, markersize = 8)

savefig(joinpath(@__DIR__,"overfit.pdf"))

linreg(x, y) = hcat(fill!(similar(x), 1), x) \ y
a,b = linreg(x[1:end-1], y[1:end-1])

plot(xs, b*xs .+ a, legend = false, linewidth = 2, grid = true)
plot!(x[1:end-1],y[1:end-1], seriestype = :scatter, markercolor = :lightgreen, markersize = 8)
savefig(joinpath(@__DIR__,"underfit.pdf"))

plot(xs, 2*xs.^2, legend = false, linewidth = 2, grid = true)
plot!(x[1:end-1],y[1:end-1], seriestype = :scatter, markercolor = :lightgreen, markersize = 8)
savefig(joinpath(@__DIR__,"justright.pdf"))

clipboard(mdtable(hcat(x,y)))
clipboard(string(x))
clipboard(string(y))

max_stars_repo_name: seanrattana/courses

max_issues_repo_name: seanrattana/courses

max_forks_repo_name: seanrattana/courses


text:
@testset "Constraint Derivatives" begin

    # Test gradient and Hessian of constraints
    T = Float64
    N = 10
    dt = 0.1
    model = UnicycleGame(p=3)
    probsize = ProblemSize(N,model)
    pdtraj = PrimalDualTraj(probsize, dt, f=ones, amplitude=0.1)
    game_con = GameConstraintValues(probsize)

    u_max =  ones(SVector{model.m,T})
    u_min = -ones(SVector{model.m,T})
    add_control_bound!(game_con, u_max, u_min)

    cval = game_con.control_conval[1]
    # Change the duals
    for k = 1:N-2
        cval.λ[k] = k*ones(2model.m)
    end
    # Evaluate the constraint and its Jacobian
    TrajectoryOptimization.evaluate!(cval, pdtraj.pr)
    @test cval.vals[1] == [-0.9*ones(model.m); -1.1*ones(model.m)]
    @test cval.vals[end] == [-0.9*ones(model.m); -1.1*ones(model.m)]
    TrajectoryOptimization.jacobian!(cval, pdtraj.pr)
    @test cval.jac[1] == [I(model.m) ;-I(model.m)]
    @test cval.jac[end] == [I(model.m) ;-I(model.m)]
    TrajectoryOptimization.cost_expansion!(cval)
    Iρ = Diagonal((cval.vals[1] .>= 0) .| (cval.λ[1] .> 0))*cval.μ[1][1]
    @test cval.grad[1] == (cval.jac[1]'*cval.λ[1] + cval.jac[1]'*Iρ*cval.vals[1])
    @test cval.hess[1] == cval.jac[1]'*Iρ*cval.jac[1]
    Iρ = Diagonal((cval.vals[end] .>= 0) .| (cval.λ[end] .> 0))*cval.μ[end][1]
    @test cval.grad[end] == (cval.jac[end]'*cval.λ[end] + cval.jac[end]'*Iρ*cval.vals[end])
    @test cval.hess[end] == cval.jac[end]'*Iρ*cval.jac[end]

end


text:
Base.iterate(a::AbstractIndex) = isempty(a) ? nothing : (first(a), 1)
function Base.iterate(a::AbstractIndex, state::Int)
    if state < length(a)
        newstate = state + 1
        return @inbounds(values(a)[newstate]), newstate
    else
        return nothing
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ e9dcc3e9-7c0f-40aa-b633-1242aacca66c
using Pkg; Pkg.activate("tmp", shared=true);

# ╔═╡ 545b2750-1bc5-11ec-1d6a-275133eec22d
using PythonCall, RDatasets;

# ╔═╡ dbf0b8fc-d9d1-492a-892e-79f1ea86c345
iris = dataset("datasets", "iris");

# ╔═╡ 41f6b405-3436-4058-979e-bcf7be661e9a
sns = pyimport("seaborn"); sns.set_theme();

# ╔═╡ 9f8ddcfe-778e-4e33-9308-c4acb409770f
sns.pairplot(pytable(iris), hue="Species")

# ╔═╡ Cell order:
# ╠═e9dcc3e9-7c0f-40aa-b633-1242aacca66c
# ╠═545b2750-1bc5-11ec-1d6a-275133eec22d
# ╠═dbf0b8fc-d9d1-492a-892e-79f1ea86c345
# ╠═41f6b405-3436-4058-979e-bcf7be661e9a
# ╠═9f8ddcfe-778e-4e33-9308-c4acb409770f


text:
# Configuration for running spatial simulation
@everywhere simtype = 2
@everywhere const N_list = [25,50,100]        # Meta-population size list
const num_trials = 1
const mutation_stddev_list = [0.005,0.01,0.02,0.04]
#const N_mut_list = [0.1,0.2,0.3 ]
const num_subpops=1
const num_attributes_list = [5]        # number attributes for quantitative representation
const ngens = 50       # Generations after burn-in
const burn_in= 2.0    # generations of burn_in as a multiple of N
const ideal = 1.0
const fit_slope = 1.0
const neutral=false


text:
#Get current time since t0
time_since(t0) = time() - t0

#String to time
function time_string_to_seconds( time::String)

    #Split time string by ':'
    v = split(time , ":" , limit = 3)

    #Calculate time in seconds
    t_seconds = 0
    for i in 1:3
        t_seconds += parse(Float64 , v[i])*60^(3-i)
    end

    return t_seconds
end

#String to time
function time_seconds_to_string( seconds::Float64)

    #
    hour = round( Int, seconds / (60*60))
    seconds %= (60*60)
    minutes = round(Int , seconds / 60, )
    seconds %= 60

    return string(  lpad(hour,2,"0") , ":" ,lpad(minutes,2,"0") , ":" , lpad(round(Int , seconds),2,"0") )

end


text:
config = MOI.Test.Config()
optimize!(mock) = MOIU.mock_optimize!(
    mock, MOI.INFEASIBLE, MOI.NO_SOLUTION, MOI.INFEASIBILITY_CERTIFICATE)
for mock in mocks(optimize!)
    Tests.choi_test(mock, config)
end

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
@test !Base.is_unix(:Windows)
@test !Base.is_linux(:Windows)
@test Base.is_linux(:Linux)
@test Base.is_windows(:Windows)
@test Base.is_windows(:NT)
@test !Base.is_windows(:Darwin)
@test Base.is_apple(:Darwin)
@test Base.is_apple(:Apple)
@test !Base.is_apple(:Windows)
@test Base.is_unix(:Darwin)
@test Base.is_unix(:FreeBSD)
@test_throws ArgumentError Base.is_unix(:BeOS)
if !is_windows()
    @test Sys.windows_version() === (0, 0)
else
    @test (Sys.windows_version()::Tuple{Int,Int})[1] > 0
end

@test (@static true ? 1 : 2) === 1
@test (@static false ? 1 : 2) === 2
@test (@static if true 1 end) === 1
@test (@static if false 1 end) === nothing
@test (@static true && 1) === 1
@test (@static false && 1) === false
@test (@static true || 1) === true
@test (@static false || 1) === 1