text:
using Mechanisms
using Test

# tests here
function mytest()
    p =  [25.806923155616076,7.206840543271378,25.433307805225706,
          25.581872903631492,43.76124192997665,17.31237766430696,
          4.328061708805246,-18.454401816641795,56.391854888638704,
          1.7725777135059568,2.406537344482028,2.905924030160018,
          3.3833860876551376,3.8836626455580485,4.520954645032934]

    D, f, g = getOptimizationProblem(PTS_VERTICAL_LINE)

    return f( p ) + sum(g( p ))
end

@test mytest() < 1e-10


text:
"""
    Poisson(λ)

A *Poisson distribution* descibes the number of independent events occurring within a unit time interval, given the average rate of occurrence `λ`.

```math
P(X = k) = \\frac{\\lambda^k}{k!} e^{-\\lambda}, \\quad \\text{ for } k = 0,1,2,\\ldots.
```

```julia
Poisson()        # Poisson distribution with rate parameter 1
Poisson(lambda)       # Poisson distribution with rate parameter lambda

params(d)        # Get the parameters, i.e. (λ,)
mean(d)          # Get the mean arrival rate, i.e. λ
```

External links:

* [Poisson distribution on Wikipedia](http://en.wikipedia.org/wiki/Poisson_distribution)

"""
struct Poisson{T<:Real} <: DiscreteUnivariateDistribution
    λ::T

    Poisson{T}(λ::Real) where {T} = (@check_args(Poisson, λ >= zero(λ)); new{T}(λ))
end

Poisson(λ::T) where {T<:Real} = Poisson{T}(λ)
Poisson(λ::Integer) = Poisson(Float64(λ))
Poisson() = Poisson(1.0)

@distr_support Poisson 0 (d.λ == zero(typeof(d.λ)) ? 0 : Inf)

#### Conversions
convert(::Type{Poisson{T}}, λ::S) where {T <: Real, S <: Real} = Poisson(T(λ))
convert(::Type{Poisson{T}}, d::Poisson{S}) where {T <: Real, S <: Real} = Poisson(T(d.λ))

### Parameters

params(d::Poisson) = (d.λ,)
@inline partype(d::Poisson{T}) where {T<:Real} = T

rate(d::Poisson) = d.λ

### Statistics

mean(d::Poisson) = d.λ

mode(d::Poisson) = floor(Int,d.λ)

function modes(d::Poisson)
    λ = d.λ
    isinteger(λ) ? [round(Int, λ) - 1, round(Int, λ)] : [floor(Int, λ)]
end

var(d::Poisson) = d.λ

skewness(d::Poisson) = one(typeof(d.λ)) / sqrt(d.λ)

kurtosis(d::Poisson) = one(typeof(d.λ)) / d.λ

function entropy(d::Poisson{T}) where T<:Real
    λ = rate(d)
    if λ == zero(T)
        return zero(T)
    elseif λ < 50
        s = zero(T)
        λk = one(T)
        for k = 1:100
            λk *= λ
            s += λk * lgamma(k + 1) / gamma(k + 1)
        end
        return λ * (1 - log(λ)) + exp(-λ) * s
    else
        return log(2 * pi * ℯ * λ)/2 -
               (1 / (12 * λ)) -
               (1 / (24 * λ * λ)) -
               (19 / (360 * λ * λ * λ))
    end
end

### Evaluation

@_delegate_statsfuns Poisson pois λ

struct RecursivePoissonProbEvaluator <: RecursiveProbabilityEvaluator
    λ::Float64
end

RecursivePoissonProbEvaluator(d::Poisson) = RecursivePoissonProbEvaluator(rate(d))
nextpdf(s::RecursivePoissonProbEvaluator, p::Float64, x::Integer) = p * s.λ / x

Base.broadcast!(::typeof(pdf), r::AbstractArray, d::Poisson, rgn::UnitRange) =
    _pdf!(r, d, rgn, RecursivePoissonProbEvaluator(d))
function Base.broadcast(::typeof(pdf), d::Poisson, X::UnitRange)
    r = similar(Array{promote_type(partype(d), eltype(X))}, axes(X))
    r .= pdf.(Ref(d),X)
end

function mgf(d::Poisson, t::Real)
    λ = rate(d)
    return exp(λ * (exp(t) - 1))
end

function cf(d::Poisson, t::Real)
    λ = rate(d)
    return exp(λ * (cis(t) - 1))
end

### Fitting

struct PoissonStats <: SufficientStats
    sx::Float64   # (weighted) sum of x
    tw::Float64   # total sample weight
end

suffstats(::Type{<:Poisson}, x::AbstractArray{T}) where {T<:Integer} = PoissonStats(sum(x), length(x))

function suffstats(::Type{<:Poisson}, x::AbstractArray{T}, w::AbstractArray{Float64}) where T<:Integer
    n = length(x)
    n == length(w) || throw(DimensionMismatch("Inconsistent array lengths."))
    sx = 0.
    tw = 0.
    for i = 1 : n
        @inbounds wi = w[i]
        @inbounds sx += x[i] * wi
        tw += wi
    end
    PoissonStats(sx, tw)
end

fit_mle(::Type{<:Poisson}, ss::PoissonStats) = Poisson(ss.sx / ss.tw)

## samplers

# TODO: remove RFunctions dependency once Poisson has been fully implemented
# Currently depends on a quantile function for one option
@rand_rdist(Poisson)
rand(d::Poisson) = convert(Int, StatsFuns.RFunctions.poisrand(d.λ))

# algorithm from:
#   J.H. Ahrens, U. Dieter (1982)
#   "Computer Generation of Poisson Deviates from Modified Normal Distributions"
#   ACM Transactions on Mathematical Software, 8(2):163-179
# TODO: implement poisson sampler
function rand(rng::AbstractRNG, d::Poisson)
    μ = d.λ
    if μ >= 10.0  # Case A

        s = sqrt(μ)
        d = 6.0*μ^2
        L = floor(Int64, μ-1.1484)

        # Step N
        T = randn(rng)
        G = μ + s*T

        if G >= 0.0
            K = floor(Int64, G)
            # Step I
            if K >= L
                return K
            end

            # Step S
            U = rand(rng)
            if d*U >= (μ-K)^3
                return K
            end

            # Step P
            px,py,fx,fy = procf(μ,K,s)

            # Step Q
            if fy*(1-U) <= py*exp(px-fx)
                return K
            end
        end

        while true
            # Step E
            E = randexp(rng)
            U = rand(rng)
            U = 2.0*U-1.0
            T = 1.8+copysign(E,U)
            if T <= -0.6744
                continue
            end

            K = floor(Int64, μ + s*T)
            px,py,fx,fy = procf(μ,K,s)
            c = 0.1069/μ

            # Step H
            if c*abs(U) <= py*exp(px+E)-fy*exp(fx+E)
                return K
            end
        end
    else # Case B
        # Ahrens & Dieter use a sequential method for tabulating and looking up quantiles.
        # TODO: check which is more efficient.
        return quantile(d,rand(rng))
    end
end

# Procedure F
function procf(μ,K,s)
    ω = 0.3989422804014327/s
    b1 = 0.041666666666666664/μ
    b2 = 0.3*b1^2
    c3 = 0.14285714285714285*b1*b2
    c2 = b2 - 15.0*c3
    c1 = b1 - 6.0*b2 + 45.0*c3
    c0 = 1.0 - b1 + 3.0*b2 - 15.0*c3

    if K < 10
        px = -μ
        py = μ^K/factorial(K) # replace with loopup?
    else
        δ = 0.08333333333333333/K
        δ -= 4.8*δ^3
        V = (μ-K)/K
        px = K*log1pmx(V) - δ # avoids need for table
        py = 0.3989422804014327/sqrt(K)

    end
    X = (K-μ+0.5)/s
    X2 = X^2
    fx = -0.5*X2 # missing negation in paper
    fy = ω*(((c3*X2+c2)*X2+c1)*X2+c0)
    return px,py,fx,fy
end

ext: jl


text:
using Tools
using Test

@testset "Tools.jl" begin
    # Write your tests here.
end


text:
using Logging
using Printf
using Random
using Statistics

using ArgParse
using LoggingExtras
using DataDeps
using Flux
using JLD2
using OrdinaryDiffEq
using Zygote

using Oceananigans
using OceanParameterizations
using FreeConvection

ENV["DATADEPS_ALWAYS_ACCEPT"] = "true"
ENV["GKSwstype"] = "100"

function parse_command_line_arguments()
    settings = ArgParseSettings()

    @add_arg_table! settings begin
        "--grid-points"
            help = "Number of vertical grid points in the trained neural differential equation (LES data will be coarse-grained to this resolution)."
            default = 32
            arg_type = Int

        "--base-parameterization"
            help = "Base parameterization to use for the NDE. Options: nothing, convective_adjustment"
            default = "nothing"
            arg_type = String

        "--time-stepper"
            help = "DifferentialEquations.jl time stepping algorithm to use."
            default = "Tsit5"
            arg_type = String

        "--name"
            help = "Experiment name (also determines name of output directory)."
            default = "testing_free_conection_nde"
            arg_type = String

        "--epochs"
            help = "Number of epochs per optimizer to train on the full time series."
            default = 10
            arg_type = Int

        "--conv"
            help = "Toggles filter dim/if a convolutional layer is included in the NN architecture. conv > 1 --> layer is added"
            default = 0
            arg_type = Int

        "--spatial_causality"
            help = "Toggles how/if spatial causality is enforced in dense layer models. Empty string -> not enforced."
            default = ""
            arg_type = String

        "--training-simulations"
            help = "Simulation IDs (list of integers separated by spaces) to train the neural differential equation on." *
                   "All other simulations will be used for testing/validation."
            action = :append_arg
            nargs = '+'
            arg_type = Int
            range_tester = (id -> id in FreeConvection.SIMULATION_IDS)
    end

    return parse_args(settings)
end

@info "Parsing command line arguments..."

args = parse_command_line_arguments()

nde_type = Dict(
    "nothing" => FreeConvectionNDE,
    "convective_adjustment" => ConvectiveAdjustmentNDE
)

Nz = args["grid-points"]
NDEType = nde_type[args["base-parameterization"]]
algorithm = Meta.parse(args["time-stepper"] * "()") |> eval
experiment_name = args["name"]
full_epochs = args["epochs"]

conv = args["conv"]
spatial_causality = args["spatial_causality"]

ids_train = args["training-simulations"][1]
ids_test = setdiff(FreeConvection.SIMULATION_IDS, ids_train)
validate_simulation_ids(ids_train, ids_test)

output_dir = joinpath(@__DIR__, experiment_name)
mkpath(output_dir)

@info "Planting loggers..."

log_filepath = joinpath(output_dir, "$(experiment_name)_training.log")
TeeLogger(
    OceananigansLogger(),
    MinLevelLogger(FileLogger(log_filepath), Logging.Info)
) |> global_logger

@info "Architecting neural network..."

if conv > 1
    NN = Chain(
           x -> reshape(x, Nz, 1, 1, 1),
           Conv((conv, 1), 1 => 1, relu),
           x -> reshape(x, Nz-conv+1),
           Dense(Nz-conv+1, 4Nz, relu),
           Dense(4Nz, 4Nz, relu),
           Dense(4Nz, Nz-1))
else
    NN = Chain(Dense(Nz, 4Nz, relu),
               Dense(4Nz, 4Nz, relu),
               Dense(4Nz, Nz-1))
end

function free_convection_neural_network(input)
    wT_interior = NN(input.temperature)
    wT = cat(input.bottom_flux, wT_interior, input.top_flux, dims=1)
    return wT
end

@info "Loading training data..."

data = load_data(ids_train, ids_test, Nz)

training_datasets = data.training_datasets
coarse_training_datasets = data.coarse_training_datasets

@info "Animating ⟨T⟩(z,t) and ⟨w'T⟩(z,t) training data..."

for id in keys(training_datasets)
    filepath = joinpath(output_dir, "free_convection_training_data_$id")
    if !isfile(filepath * ".mp4") || !isfile(filepath * ".gif")
        animate_training_data(training_datasets[id], coarse_training_datasets[id]; filepath, frameskip=5)
    end
end

@info "Wrangling (T, wT) training data..."

input_training_data = wrangle_input_training_data(coarse_training_datasets)
output_training_data = wrangle_output_training_data(coarse_training_datasets)

@info "Scaling features..."

T_training_data = reduce(hcat, input.temperature for input in input_training_data)
wT_training_data = output_training_data

@assert size(wT_training_data, 1) == size(T_training_data, 1) + 1
@assert size(wT_training_data, 2) == size(T_training_data, 2)

T_scaling = ZeroMeanUnitVarianceScaling(T_training_data)
wT_scaling = ZeroMeanUnitVarianceScaling(wT_training_data)

input_training_data = [rescale(i, T_scaling, wT_scaling) for i in input_training_data]
output_training_data = wT_scaling.(output_training_data)

@info "Batching training data..."

n_training_data = length(input_training_data)
training_data = [(input_training_data[n], output_training_data[:, n]) for n in 1:n_training_data] |> shuffle
data_loader = Flux.Data.DataLoader(training_data, batchsize=n_training_data, shuffle=true)

n_obs = data_loader.nobs
batch_size = data_loader.batchsize
n_batches = ceil(Int, n_obs / batch_size)
@info "Training data loader contains $n_obs pairs of observations (batch size = $batch_size)."

@info "Training neural network on fluxes: ⟨T⟩(z) -> ⟨w'T⟩(z) mapping..."

causal_penalty = nothing

if spatial_causality == "soft"
    ps = Flux.params(NN)

    dense_layer_idx = 1 + Int(conv > 1) * 3
    dense_layer_params_idx = 1 + Int(conv > 1) * 2

    nrows, ncols = size(ps[dense_layer_params_idx])
    mask = [x < y ? true : false for x in 1:nrows, y in 1:ncols]

    causal_penalty() = sum(abs2, NN[dense_layer_idx].W[mask])

    nn_loss(input, output) = Flux.mse(free_convection_neural_network(input), output) + causal_penalty()
else
    nn_loss(input, output) = Flux.mse(free_convection_neural_network(input), output)
end

nn_training_set_loss(training_data) = mean(nn_loss(input, output) for (input, output) in training_data)

function nn_callback()
    μ_loss = nn_training_set_loss(training_data)
    @info @sprintf("Training free convection neural network... training set MSE loss: μ_loss::%s = %.10e", typeof(μ_loss), μ_loss)
    return μ_loss
end

epochs = 10
optimizers = [ADAM(), Descent(1e-2)]

for opt in optimizers, e in 1:epochs, (i, mini_batch) in enumerate(data_loader)
    @info "Training heat flux neural network with $(typeof(opt))(η=$(opt.eta))... (epoch $e/$epochs, mini-batch $i/$n_batches)"
    Flux.train!(nn_loss, Flux.params(NN), mini_batch, opt, cb=Flux.throttle(nn_callback, 5))
end

@info "Animating what the neural network has learned so far (T -> wT mapping)..."

for (id, ds) in coarse_training_datasets
    filepath = joinpath(output_dir, "learned_free_convection_initial_guess_$id")
    if !isfile(filepath * ".mp4") || !isfile(filepath * ".gif")
        animate_learned_free_convection(ds, NN, free_convection_neural_network, NDEType, algorithm, T_scaling, wT_scaling,
                                        filepath=filepath, frameskip=5)
    end
end

@info "Saving initial neural network weights to disk..."

initial_nn_filepath = joinpath(output_dir, "free_convection_initial_neural_network.jld2")

jldopen(initial_nn_filepath, "w") do file
    file["grid_points"] = Nz
    file["neural_network"] = NN
    file["T_scaling"] = T_scaling
    file["wT_scaling"] = wT_scaling
end

@info "Training neural differential equation on incrementally increasing time spans..."

nn_history_filepath = joinpath(output_dir, "neural_network_history.jld2")

training_iterations = (1:20, 1:5:101, 1:10:201, 1:20:401, 1:40:801)
training_epochs     = (10,   10,      10,       10,       10)
opt = ADAM()

for (iterations, epochs) in zip(training_iterations, training_epochs)
    @info "Training free convection NDE with iterations=$iterations for $epochs epochs  with $(typeof(opt))(η=$(opt.eta))..."
    train_neural_differential_equation!(NN, NDEType, algorithm, coarse_training_datasets, T_scaling, wT_scaling,
                                        iterations, opt, epochs, history_filepath=nn_history_filepath, causal_penalty=causal_penalty)
end

@info "Training the neural differential equation on the entire solution while decreasing the learning rate..."

burn_in_iterations = 1:9:1153
optimizers = [ADAM(1e-3), ADAM(5e-4)]

for opt in optimizers
    @info "Training free convection NDE with iterations=$burn_in_iterations for $full_epochs epochs with $(typeof(opt))(η=$(opt.eta))..."
    train_neural_differential_equation!(NN, NDEType, algorithm, coarse_training_datasets, T_scaling, wT_scaling,
                                        burn_in_iterations, opt, full_epochs, history_filepath=nn_history_filepath, causal_penalty=causal_penalty)
end

@info "Saving final neural network weights to disk..."

final_nn_filepath = joinpath(output_dir, "free_convection_final_neural_network.jld2")

jldopen(final_nn_filepath, "w") do file
    file["grid_points"] = Nz
    file["neural_network"] = NN
    file["T_scaling"] = T_scaling
    file["wT_scaling"] = wT_scaling
end


text:
__precompile__(true)
module KFEstimate

# Write your package code here.
using LinearAlgebra
using ForwardDiff
using Zygote
using Random
using Revise
using Plots
using ProgressBars
using Flux, Flux.Optimise
using Statistics
using Distributions
import Random: rand

export
    AbstractFilter,
    KalmanFilter,
    ExtendedKalmanFilter,
    AbstractParamFilter,
    ParamKalmanFilter,
    ExtendedParamKalmanFilter,
    State
include("filters.jl")

export
    dynamic,
    param_dynamic,
    prediction,
    param_prediction,
    observation,
    param_observation,
    correction,
    param_correction,
    state_mse
include("ekf_functions.jl")
include("kf_functions.jl")

export
    kf_likelihood,
    ekf_likelihood
include("likelihood.jl")

export
    run_simulation,
    run_filter,
    run_param_filter,
    run_kf_gradient,
    run_ekf_gradient
include("run.jl")

export
    unpack
include("utils.jl")
end


text:
"""
    rstar([rng ,]classif::Supervised, chains::Chains; kwargs...)
    rstar([rng ,]classif::Supervised, x::AbstractMatrix, y::AbstractVector; kwargs...)

Compute the ``R^*`` convergence diagnostic of MCMC.

This implementation is an adaption of Algorithm 1 & 2, described by [^LambertVehtari2020].
Note that the correctness of the statistic depends on the convergence of the classifier used
internally in the statistic. You can inspect the training of the classifier by adjusting the
verbosity level.

[^LambertVehtari2020]: Lambert & Vehtari (2020). ``R^*``: A robust MCMC convergence diagnostic with uncertainty using gradient-boosted machines. arXiv preprint <https://arxiv.org/abs/2003.07900>.

# Keyword Arguments
* `subset = 0.8` ... Subset used to train the classifier, i.e. 0.8 implies 80% of the samples are used.
* `iterations = 10` ... Number of iterations used to estimate the statistic. If the classifier is not probabilistic, i.e. does not return class probabilities, it is advisable to use a value of one.
* `verbosity = 0` ... Verbosity level used during fitting of the classifier.

# Example
```jldoctest rstar; output = false, filter = r".*"s
using MLJModels

XGBoost = @load XGBoostClassifier verbosity=0
chn = Chains(fill(4, 100, 2, 3))

Rs = rstar(XGBoost(), chn; iterations=20)
R = round(mean(Rs); digits=0)

# output

1.0
```
"""
function rstar(rng::Random.AbstractRNG, classif::MLJModelInterface.Supervised, x::AbstractMatrix, y::AbstractVector{Int}; iterations = 10, subset = 0.8, verbosity = 0)

    size(x,1) != length(y) && throw(DimensionMismatch())
    iterations >= 1 && ArgumentError("Number of iterations has to be positive!")

    if iterations > 1 && classif isa MLJModelInterface.Deterministic
        @warn("Classifier is not a probabilistic classifier but number of iterations is > 1.")
    elseif iterations == 1 && classif isa MLJModelInterface.Probabilistic
        @warn("Classifier is probabilistic but number of iterations is equal to one.")
    end

    N = length(y)
    K = length(unique(y))

    # randomly sub-select training and testing set
    Ntrain = round(Int, N*subset)
    Ntest = N - Ntrain

    ids = Random.randperm(rng, N)
    train_ids = view(ids, 1:Ntrain)
    test_ids = view(ids, (Ntrain+1):N)

    # train classifier using XGBoost
    fitresult, _ = MLJModelInterface.fit(classif, verbosity, Tables.table(x[train_ids,:]), MLJModelInterface.categorical(y[train_ids]))

    xtest = Tables.table(x[test_ids,:])
    ytest = view(y, test_ids)

    Rstats = map(i -> K*rstar_score(rng, classif, fitresult, xtest, ytest), 1:iterations)
    return Rstats
end

function rstar(classif::MLJModelInterface.Supervised, x::AbstractMatrix, y::AbstractVector{Int}; kwargs...)
    rstar(Random.GLOBAL_RNG, classif, x, y; kwargs...)
end

function rstar(classif::MLJModelInterface.Supervised, chn::Chains; kwargs...)
    return rstar(Random.GLOBAL_RNG, classif, chn; kwargs...)
end

function rstar(rng::Random.AbstractRNG, classif::MLJModelInterface.Supervised, chn::Chains; kwargs...)
    nchains = size(chn, 3)
    nchains <= 1 && throw(DimensionMismatch())

    # collect data
    x = Array(chn)
    y = repeat(chains(chn); inner = size(chn,1))

    return rstar(rng, classif, x, y; kwargs...)
end

function rstar_score(rng::Random.AbstractRNG, classif::MLJModelInterface.Probabilistic, fitresult, xtest, ytest)
    pred = get.(rand.(Ref(rng), MLJModelInterface.predict(classif, fitresult, xtest)))
    return mean(((p,y),) -> p == y, zip(pred, ytest))
end

function rstar_score(rng::Random.AbstractRNG, classif::MLJModelInterface.Deterministic, fitresult, xtest, ytest)
    pred = MLJModelInterface.predict(classif, fitresult, xtest)
    return mean(((p,y),) -> p == y, zip(pred, ytest))
end


text:
################################################################################

"""

    getSignals(shParams::Dict)

# Description
Read EDF file from shell arguments. Return dataframe of signals.

"""
function getSignals(shParams::Dict)
  if haskey(shParams, "indir") && haskey(shParams, "file")
    return getSignals( string(shParams["indir"], shParams["file"]) )
  else
    @error "Variables are not defined in dictionary"
  end
end

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

"""

    getSignals(edfFile::S) where S <: String

# Description
Read EDF file. Return dataframe of signals.

"""
function getSignals(edfFile::S) where S <: String
  # read edf file
  @info "Reading EDF file..."
  edfRecord = EDF.read(edfFile)

  # load signals
  edfDf = DataFrames.DataFrame()
  signalLength = length(edfRecord.signals)

  for ι ∈ 1:signalLength
    try
      if edfRecord.signals[ι].header.label ∉ ["A", "MK", "ECG", "EKG"]
        edfDf[!, edfRecord.signals[ι].header.label] = EDF.decode(edfRecord.signals[ι])
      end
    catch
      # no catch
    end
  end

  #  collect record frecuency
  recordFreq = getedfRecordFreq(edfRecord)

  # collect timestamp
  startTime = getedfStart(edfRecord)

  return edfDf, startTime, recordFreq
end

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

"obtain start time"
function getedfStart(edfRecord)
  # collect timestamp
  startTime = edfRecord.header.start |> Dates.Time
end

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

"obtain record frequency"
function getedfRecordFreq(edfRecord)
  # collect recording frecuency
  signalLength = length(edfRecord.signals)

  recordFreq = [edfRecord.signals[1].header.samples_per_record]
  for ι ∈ 2:signalLength
    try
      if edfRecord.signals[ι].header.label ∉ ["A", "MK", "ECG", "EKG"]
        recordFreq = [recordFreq; edfRecord.signals[ι].header.samples_per_record]
      end
    catch
      # no catch
    end
  end
  # explicit type conversion
  return convert.(Int64, recordFreq)
end

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

max_stars_repo_name: DanielRivasMD/MindReader

max_issues_repo_name: DanielRivasMD/MindReader

max_forks_repo_name: DanielRivasMD/MindReader


text:
    using Dart, Optim
    using Test
    function Hosaki(x)
        T = eltype(x)
        pol = Base.Math.@horner(x[1], T(1), -T(8), T(7), -T(7)/3, T(1)/4)

        fx = exp(-x[2])*pol*x[2]^2
    end

    L = (f, lb, ub, x) -> optimize(f, lb, ub, x, Fminbox(BFGS()))

    hosaki_res = mlsl(Hosaki, L, [0.0, 0.0], [10.0, 10.0]; N=5, sigma=4.0, try_ratio=0.1, maxiter=100, getmin=Optim.minimizer, getminf=Optim.minimum)

    function HolderTable(x)
        fact1 = sin(x[1])*cos(x[2])
        fact2 = exp(abs(1-sqrt(x[1]^2+x[2]^2)/π))
        -abs(fact1*fact2)
    end
    holder_res=mlsl(HolderTable, L, [-10.0, -10.0], [10.0, 10.0]; N=5, sigma=4.0, try_ratio=0.1, maxiter=100, getmin=Optim.minimizer, getminf=Optim.minimum)

    function Rastrigin(x)
        T = eltype(x)
        d = length(x)
        10*d + sum(xᵢ^2-10*cos(2π*xᵢ) for xᵢ in x)::T
    end
    n_rastrigin = 3
    rastrigin_x, rastrigin_fx, from_x, from_fx=mlsl(Rastrigin, L, fill(-5.12, n_rastrigin), fill(5.12, n_rastrigin); N=5, sigma=4.0, try_ratio=0.1, maxiter=100, getmin=Optim.minimizer, getminf=Optim.minimum)

    using LinearAlgebra
    function Dropwave(x)
        normx = norm(x)
        num = 1 + cos(12*norm(x))
        den = normx^2/2 + 2
        -num/den
    end
    dropwave_x, dropwave_fx, from_x, from_fx=mlsl(Dropwave, L, fill(-5.12, 2), fill(5.12, 2); N=5, sigma=4.0, try_ratio=0.1, maxiter=100, getmin=Optim.minimizer, getminf=Optim.minimum, try_center=false)

    function Damavandi(x)
        (1-abs(sin(π*(x[1]-2))*sin(π*(x[2]-2))/(π^2*(x[1]-2)*(x[2]-2)))^5)*(2+(x[1]-7)^2+2*(x[2]-7)^2)
    end
    Damavandi_x, Damavandi_fx, from_x, from_fx=mlsl(Damavandi, L, fill(0.0, 2), fill(14.0, 2); N=5, sigma=2.0, try_ratio=0.1, maxiter=400, getmin=Optim.minimizer, getminf=Optim.minimum, try_center=false)

    @testset "Dart.jl" begin

    end


text:
#=
Copyright 2022 INSIGNEO Institute for in silico Medicine

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

   http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

"""
	saveTempData(t :: Float64, vessels :: Array{Vessel, 1}, counter :: Int)

Save current solution to `temporary` arrays.
"""
function saveTempData(t :: Float64, vessels :: Array{Vessel, 1}, counter :: Int)
	for v in vessels
		saveTempData(t, v, counter)
	end
end

function saveTempData(t :: Float64, v :: Vessel, counter :: Int)
	v.P_t[counter,:] = [t, v.P[1], v.P[v.node2], v.P[v.node3], v.P[v.node4], v.P[end]]
	v.A_t[counter,:] = [t, v.A[1], v.A[v.node2], v.A[v.node3], v.A[v.node4], v.A[end]]
	v.Q_t[counter,:] = [t, v.Q[1], v.Q[v.node2], v.Q[v.node3], v.Q[v.node4], v.Q[end]]
	v.u_t[counter,:] = [t, v.u[1], v.u[v.node2], v.u[v.node3], v.u[v.node4], v.u[end]]
	v.c_t[counter,:] = [t, v.c[1], v.c[v.node2], v.c[v.node3], v.c[v.node4], v.c[end]]
end

"""
	transferTempToLast(vessels :: Array{Vessel, 1})

Copy (deepcopy) `temporary` arrays content to `last` arrays.
"""
function transferTempToLast(vessels :: Array{Vessel, 1})
	for v in vessels
		transferTempToLast(v)
	end
end

function transferTempToLast(v :: Vessel)
	v.A_l = deepcopy(v.A_t)
	v.P_l = deepcopy(v.P_t)
	v.Q_l = deepcopy(v.Q_t)
	v.u_l = deepcopy(v.u_t)
	v.c_l = deepcopy(v.c_t)
end

"""
	transferLastToOut(vessels :: Array{Vessel, 1})

Write `last` arrays content to `.out` files.
"""
function transferLastToOut(vessels :: Array{Vessel, 1})
	for v in vessels
		transferLastToOut(v)
	end
end

function transferLastToOut(v :: Vessel)
	lastP = v.P_l
	lastQ = v.Q_l
	lastA = v.A_l
	lastc = v.c_l
	lastu = v.u_l
	lasts = [lastP, lastQ, lastA, lastc, lastu]

	outP  = v.out_P_name
	outQ  = v.out_Q_name
	outA  = v.out_A_name
	outc  = v.out_c_name
	outu  = v.out_u_name
	outs = [outP, outQ, outA, outc, outu]

	for (a, b) in zip(outs, lasts)
		out_file = open(a, "a")
		writedlm(out_file, b, " ")
		close(out_file)
	end
end

"""
	writeResults(vessels :: Array{Vessel, 1})

Write `last` arrays content to `.last` files.
"""
function writeResults(vessels :: Array{Vessel, 1})
	for v in vessels
		writeResults(v)
	end
end

function writeResults(v :: Vessel)
	lastP = v.P_l
	lastQ = v.Q_l
	lastA = v.A_l
	lastc = v.c_l
	lastu = v.u_l
	lasts = [lastP, lastQ, lastA, lastc, lastu]

	resP  = v.last_P_name
	resQ  = v.last_Q_name
	resA  = v.last_A_name
	resc  = v.last_c_name
	resu  = v.last_u_name
	ress = [resP, resQ, resA, resc, resu]

	for (a, b) in zip(ress, lasts)
		writedlm(a, b, " ")
	end
end

"""
	writeConv(data :: Dict{Any,Any}, passed_cycles :: Int)

Write `.conv` file
"""
function writeConv(data :: Dict{Any,Any}, passed_cycles :: Int)
	f_name = join([data["project name"],".conv"])
	conv_file = open(f_name, "w")
	write(conv_file, join(["Converged after ", string(passed_cycles)," iterations\n"]))
	close(conv_file)
end


text:
struct Weight <: EdgeModifier
    w
end

struct Observed <: NodeModifier end

struct Label <: NodeModifier l end
import Base.==
==(x::Label, y::Label) = x.l == y.l

@testset "Swap Node <-> Node" begin
    n1 = Node(:a)
    n2 = Node(:b)
    @test swap_node(n1, n1, n2) == swap_node(n1, n1 => n2) == n2
    @test swap_node(n1 ^ Observed(), n1, n2) == swap_node(n1 ^ Observed(), n1 => n2) == n2 ^ Observed()
    @test swap_node(n1 * Weight(1), n1, n2) == n2 * Weight(1)
    @test swap_node(n1 ^ Label("hi"), n1 ^ Label("hi"), n2 * Observed()) == n2 * Observed()
    @test swap_node(n1, n1, n1 ^ Label("hi")) == n1 ^ Label("hi")
end

@testset "Swap Nodes in Edges" begin
    for e in (:DirectedEdge, :UndirectedEdge)
        @eval begin
            n1 = Node(:a)
            n2 = Node(:b)
            n3 = Node(:c)
            e1 = $e(n1, n2)
            e3 = $e(n3, n2)
            @test swap_node(e1, n1, n3) == e3
            @test swap_node(e1, n2, n3) == $e(n1, n3)
            @test swap_node(e1 * Weight(1), n1, n3) == e3 * Weight(1)
            @test swap_node($e(n1 ^ Label("hi"), n2), n1 ^ Label("hi"), n3 ^ Label("ho")) == $e(n3 ^ Label("ho"), n2)
        end
    end
end

@testset "Swap Edges" begin
    for e in (:DirectedEdge, :UndirectedEdge)
        @eval begin
            n1 = Node(:a)
            n2 = Node(:b)
            n3 = Node(:c)
            e1 = $e(n1, n2)
            e3 = $e(n3, n2)
            @test swap_edge(e1, e1, e3) == swap_edge(e1, e1 => e3) == e3
            @test swap_edge(e1 * Weight(1), e1, e3) == e3 * Weight(1)
            @test swap_edge(e1, e1, e1 * Weight(1)) == e1 * Weight(1)
        end
    end
end

@testset "Swap nonexistent nodes" begin
    n1 = Node(:a)
    @test_throws KeyError swap_node(n1, Node(:c), Node(:d))

    e1 = Edge(Node(:a), Node(:b))
    @test_throws KeyError swap_node(e1, Node(:c), Node(:d))
end


text:
function train_light!(x, y, model, ps, loss, opt)
    ps = Flux.Params(ps)
    gs = Flux.gradient(ps) do
        loss(x, y, model)
    end
    Flux.update!(opt, ps, gs)
end

scalar_yield(iter::Int; every = 100) = (iter % every == 0) && yield()

function log2_yield(iter::Int; every = 100)
    l2 = log2(iter)
    if (iter >= 0) && (trunc(l2) ≈ l2)
        yield()
    end
end

function smart_yield(iter::Int; every = 100)
    if (iter % every == 0)
        yield()
    else
        l2 = log2(iter)
        if (iter > 0) && (trunc(l2) ≈ l2)
            yield()
        end
    end
end

function train_loop!(   X, Y, model, loss, opt;
                        max_iters = 1000, yeild_fn = identity, record_every = 10)
    loss_record = []
    ps = Flux.params(model)
    for (iter,x,y) in zip(1:max_iters,X,Y)
        ps = Flux.Params(ps)
        gs = Flux.gradient(ps) do
            loss(x, y, model)
        end
        Flux.update!(opt, ps, gs)
        #for reporting out
        if iter % record_every == 0
            push!(loss_record)
        end
        #allow Ctrl+C interrupts after N iterations
        yeild_fn(iter)
    end
    return loss_record
end

#Do early stopping...


text:
@testset "Model matrix" begin

    using StreamModels
    using Test
    using DataStreams

    using StreamModels: @formula

    source = Data.Table((a = collect(1:10),
                         b = rand(10),
                         c = repeat(["a", "b"], inner=5)))

    @testset "single term RHS" begin

        @test all(modelmatrix(source, @formula(~ 1)) .== 1)
        @test all(modelmatrix(source, @formula(~ a)) .== float(source[:a]))
        @test all(modelmatrix(source, @formula(~ b)) .== source[:b])
        @test all(modelmatrix(source, @formula(~ c)) .== [ones(5)  zeros(5)
                                                          zeros(5) ones(5)])

    end

    @testset "star expansion" begin
        f = @formula( ~ a*b*c)

        mm = modelmatrix(source, f)

        @test size(mm) == (10, 8)
        @test mm[:,1] == source[:a]
        @test mm[:,2] == source[:b]
        @test mm[:,3] == [ones(5); zeros(5)]
        @test mm[:,4] == [zeros(5); ones(5)]
        @test mm[:,5] == source[:a] .* source[:b]
        @test mm[:,6] == source[:a] .* [zeros(5); ones(5)]
        @test mm[:,7] == source[:b] .* [zeros(5); ones(5)]
        @test mm[:,8] == source[:a] .* source[:b] .* [zeros(5); ones(5)]
    end

    @testset "arbitrary functions" begin

        @test modelmatrix(source, @formula(~ log(a)))[:,1] == log.(float(source[:a]))
        @test modelmatrix(source, @formula(~ log(1+a)))[:,1] == log.(1 .+ float(source[:a]))
        @test modelmatrix(source, @formula(~ 1 + log(a))) == [ones(10) log.(float(source[:a]))]
        @test modelmatrix(source, @formula(~ log(a+b)))[:,1] == log.(source[:a].+source[:b])
        ## @test modelmatrix(source, @formula(~ log(a*b)))[:,1] == log.(source[:a].*source[:b])

    end

end


text:
"Sets up the master problem"
function setupmaster(subproblems,patients,resources,timeslots,subMastercalendar,sets)
    env = Gurobi.Env()
    Gurobi.setparams!(env, OutputFlag=0)
    master = Model(with_optimizer(Gurobi.Optimizer,env))

    K = length(patients)
    Pg = sets.Pg
    Gp = [i for i in 1:length(Pg)]
    D = JuliaDB.select(resources,:intID)
    J = keys(subMastercalendar)
    Jd = sets.Jd
    I = sets.I

    @variable(master,lambda[1:K] >= 0)
    #@variable(master,closingtime[d in D, j in Jd[d].j] >= 0)
    # for d in D, j in Jd[d].j
    #     set_lower_bound(closingtime[d,j],getTlowerbound(timeslots,d,j))
    # end

    #@objective(master, Min, sum(closingtime[d,j] for d in D, j in Jd[d].j) + sum(1000000*lambda[m] for m in 1:K))
    @objective(master, Min, sum(1000000*lambda[m] for m in 1:K))
    #@constraint(master,consref_offtime[d in D, j in Jd[d].j], sum(lambda[m]*1000 for m in 1:K) <= closingtime[d,j]) #TODO Only for consultations

    @constraint(master,convexitycons[g in Gp],
       sum(lambda[m] for m in 1:K if m in Pg[g].patients) == length(Pg[g].patients) )

    @constraint(master,consref_onepatient[ d in D, j in Jd[d].j, i in I[d,j].i],
        0 <=1 )

    return Masterproblem(master,[],consref_onepatient,convexitycons,lambda,[],I,env)
    #return Masterproblem(master,consref_offtime,consref_onepatient,convexitycons,lambda,closingtime,I,env)

end
"Return the sum of the maximal path of visits that visit x must preceed"
function sumpath(Tdelta,x)
    sum = []
    for y in (Tdelta |> @filter(_.visit1 == x) |> @map(_.visit2) |> collect)
        push!(sum,Tdelta[x,y].delta+sumpath(Tdelta,y))
    end
    if length(sum)> 0
        return maximum(sum)
    else
        return 0
    end
end

"Sorts visits in in the scheduling order for the heuristic"
function sortvisit(V_input,Tdelta_input)
    V = copy(V_input)
    result = []
    Tdelta = Tdelta_input |> @filter(_.visit1 in V_input && _.visit2 in V_input) |> NDSparse
    curr = "Consultation"
    V = filter!(x-> x != curr,V)
    push!(result,curr)
    while length(V)>0
        temp = Tdelta |> @filter(_.visit1 in V && _.visit2 in result) |> @orderby_descending(sumpath(Tdelta,_.visit1))|> @map(_.visit1)|> collect
        if length(temp)> 0
            curr = first(temp)
        else
            curr = rand(V)
        end
        while true
            temp2 = Tdelta[curr,:] |> @filter(!(_.visit2 in result))|> @orderby_descending(sumpath(Tdelta,_.visit2))|> @map(_.visit2)|> collect
            if length(temp2)> 0
                curr = first(temp2)
            else
                break
            end
        end
        push!(result,curr)
        V = filter!(x-> x != curr,V)
    end
    result
end
"used to check if a timeslot on dayID, with startTime and endTime is feasible with the current already plannedslots for the patient, such that there are no overlap and the visit order given in timedelta is kept."
function isfree(plannedslots,timeDelta,type,dayID,startTime,endTime)
    for p in plannedslots
        if p.dayID == dayID && ((startTime< p.endTime && endTime > p.startTime) || (haskey(timeDelta,(type,p.type)) && endTime > p.startTime))
            return false
        end
    end
    true
end

"Checks if the day with dayID, is already used in the plan, such that the patient is already at hospital"
function daystoopen(plannedslots,dayID)
    if dayID in (x->x.dayID).(plannedslots)
        return 0
    else
        return 1
    end
end
"Calculates the latest day that a visit of the type 'type' can be scheduled according to the data of timeDelta"
function lastestday(plannedslots,timeDelta,type)
    if length(plannedslots) > 0
        x = (x-> x.dayID -get(timeDelta,(type,x.type),(delta = -1,)).delta).(plannedslots)
        return minimum(x)
    end
    return 100000
end

"Builds schedule using a greedy heuristic. For each type of visit (in sorted order), the algorithm places the last one and then plans the most constrained visit (according to the already planned visits) next, until all visits are planned."
function buildschedule(timeslots,appointments,types, timeDelta)
    plannedslots = []
    for type in types
        latest = lastestday(plannedslots,timeDelta,type)
        x = @from i in timeslots begin
            @where i.type == type && i.dayID <= latest && isfree(plannedslots,timeDelta,i.type,i.dayID,i.startTime,i.endTime)#TODO her mangler vi at tjekke for starttid
            @left_outer_join j in appointments on i.timeslotID equals j.timeslotID
            @where i.timeslotID != j.timeslotID
            @orderby daystoopen(plannedslots,i.dayID), descending(i.dayID)
            @select a=i.timeslotID
            @collect
        end
        length(x) == 0 && return []
        y = first(x)
        push!(plannedslots,timeslots[y])
    end
    return plannedslots
end

"Generates initial columns for the masterproblem using the buildschedule greedy heuristic that ensures that no columns overlap."
function generateInitialColumns!(masterproblem::Masterproblem,sets,timeslots,timeDelta)

    appointments = table((timeslots |> @filter(_.booked) |> @map((timeslotID = _.timeslotID ,resourceID=_.resourceID ,visitID= 0)) |> collect);pkey=[:timeslotID])
    while true
        done = true
        for group in 1:length(sets.Pg)
            types = sets.Pg[group].types
            for patient in sets.Pg[group].patients
                plannedslots = buildschedule(timeslots,appointments,types,timeDelta)
                if length(plannedslots)> 0
                    done = false
                    appointments = append!(rows(appointments),(timeslotID = (x->x.timeslotID).(plannedslots),resourceID = (x->x.resourceID).(plannedslots),visitID  =-1 ))
                    addcolumntomaster!(masterproblem,plannedslots,group,0)
                end
            end
        end
        done && break

    end
end

"adds columns from suproblems that have objective value beloew -EPSVALUE"
function addcolumntomaster!(masterproblem::Masterproblem,subproblems::Subproblems,iteration::Int64,EPSVALUE)
    done = true
    for sub in values(subproblems.pricingproblems)
        if objective_value(sub.model) < -EPSVALUE
            println("Subproblem objective value = $(JuMP.objective_value(sub.model))")
            addcolumntomaster!(masterproblem,sub,iteration)
            done = false
        end
    end
    return done
end

"adds column from a single pricingproblems if the objective is below -EPSVALUE"
function addcolumntomaster!(masterproblem::Masterproblem,pricingproblem::PricingProblem,iteration::Int64,EPSVALUE)
        if objective_value(pricingproblem.model) < -EPSVALUE
            println("Subproblem objective value = $(JuMP.objective_value(pricingproblem.model))")
            addcolumntomaster!(masterproblem,pricingproblem,iteration)
            return false
        end
    return true
end

"adds column from a single pricingproblems"
function addcolumntomaster!(masterproblem::Masterproblem,pricingproblem::PricingProblem,iteration::Int64)

    touchedconstraints = ConstraintRef[]
    constraint_coefficients = Float64[]

    for t in getIndexofPositiveVariables(pricingproblem.tvars)

            push!(touchedconstraints,masterproblem.consref_offtime[t])
            push!(constraint_coefficients,value(pricingproblem.tvars[t]))
    end
    for x in getIndexofPositiveVariables(pricingproblem.xvars)

        push!(touchedconstraints,masterproblem.consref_onepatient[x[2:end]])
        push!(constraint_coefficients,1)
    end

    push!(touchedconstraints,masterproblem.convexitycons[pricingproblem.intID])
    push!(constraint_coefficients,1)

    push!(masterproblem.lambda,@variable(
        masterproblem.model,
        lower_bound = 0,
        base_name = "lambda_new[$(pricingproblem.intID),$(iteration)]_$(length(masterproblem.lambda))" ,

    ))
    JuMP.set_objective_coefficient(masterproblem.model,masterproblem.lambda[end],sum(value.(pricingproblem.yvars)))
    JuMP.set_normalized_coefficient.(touchedconstraints,masterproblem.lambda[end],constraint_coefficients)
end

"Adds column to master problem that uses the timeslots given in 'timeslots' for a patientgroup. Used to add columns from Heuristic"
function addcolumntomaster!(masterproblem::Masterproblem,timeslots,patientgroup::Int64,iteration)
    touchedconstraints = ConstraintRef[]
    constraint_coefficients = Float64[]
    for timeslot in timeslots
        d = timeslot.resourceID
        j = timeslot.dayID
        i = timeslot.timeslotID
        push!(touchedconstraints,masterproblem.consref_onepatient[d,j,i])
        push!(constraint_coefficients,1)
    end
    push!(touchedconstraints,masterproblem.convexitycons[patientgroup])
    push!(constraint_coefficients,1)

    push!(masterproblem.lambda,@variable(
        masterproblem.model,
        lower_bound = 0,
        base_name = "lambda_[$(patientgroup),$(iteration)]_$(length(masterproblem.lambda))" ,

    ))
    columnprice = float(length(timeslots|> @groupby(_.dayID) |> collect))

    JuMP.set_objective_coefficient(masterproblem.model,masterproblem.lambda[end],columnprice)
    JuMP.set_normalized_coefficient.(touchedconstraints,masterproblem.lambda[end],constraint_coefficients)
end


text:
export Tensor

struct Tensor <: AbstractTensor
dim::Int          # "dimension"
T::SparseVector{Float64,Int}
get::Function
op::AbstractOrthoPoly

    # inner constructors
    function Tensor(dim::Int,mop::MultiOrthoPoly)
        tensorEntries = computeTensorizedSP(dim, mop)
        getfun(ind) = getentry(ind, tensorEntries, mop.ind, dim)
        new(dim, tensorEntries, getfun, mop)
    end

    function Tensor(dim::Int,opq::AbstractOrthoPoly)
        tensorEntries = computeTensorizedSP(dim, opq)
        getfun(ind) = getentry(ind, tensorEntries, calculateMultiIndices(1, opq.deg), dim)
        new(dim, tensorEntries, getfun, opq)
    end
end


text:
using Distributions

include("assigndist_normal.jl")
include("assigndist_uniform.jl")
include("assigndist_beta.jl")
include("assigndist_betaprime.jl")
include("assigndist_betabinomial.jl")
include("assigndist_binomial.jl")
include("assigndist_gamma.jl")
include("assigndist_frechet.jl")


text:
module Caesar

# import RoME: initfg # collision on RoME.initfg() since no parameters are given in both RoME and Caesar
import Distributions: Normal
import RoME: getRangeKDEMax2D, getLastPose, initfg

using
  RoME,
  IncrementalInference,
  Graphs,
  KernelDensityEstimate,
  Distributions,
  TransformUtils,
  CoordinateTransformations,
  Rotations,
  JLD,
  HDF5,
  JSON,
  FileIO,
  DataStructures,
  ProgressMeter,
  ImageMagick,
  ImageCore,
  DocStringExtensions,
  CloudGraphs,
  Neo4j,
  Mongo,
  LibBSON

export
  # pass through from KDE
  kde!,
  getPoints,
  getBW,
  Ndim,
  Npts,

  # pass through from IIF and RoME
  ls,
  FactorGraph,
  writeGraphPdf,
  getVert,
  getVal,
  saveplot,
  wipeBuildNewTree!,
  inferOverTree!,
  inferOverTreeR!,
  # callbacks for datalayer changes
  localapi,
  dlapi,
  # Victoria Park example -- batch
  loadVicPrkDataset,

  # passthrough variable and factor types
  Pose2,
  Point2,

  # passthrough RoME factor types
  PriorPose2,
  Pose2Pose2,
  Pose2DPoint2DBearingRange,

  # insitu component
  GenericInSituSystem,
  makeGenericInSituSys,
  InSituSystem,
  makeInSituSys,
  triggerPose,
  poseTrigAndAdd!,
  processTreeTrackersUpdates!,
  advOdoByRules,
  progressExamplePlot,
  plotPoseDict,
  plotTrckStep,
  SLAMWrapper,

  # servers
  tcpStringSLAMServer,
  tcpStringBRTrackingServer,

  # save and load data
  saveSlam,
  loadSlam,
  haselement,

  # user functions
  identitypose6fg,
  projectrbe,
  hasval,

  # repeats from RoME and IIF
  initfg,
  addNode!,
  addFactor!,

  # CloudGraphs helper functions
  insertnodefromcv!,
  checkandinsertedges!,
  getbinarraymongo,
  gettopoint,
  getdotwothree,
  bin2arr,
  fetchsubgraph!,
  getVertNeoIDs!,
  insertrobotdatafirstpose!,
  tryunpackalltypes!,
  fetchrobotdatafirstpose,
  getExVertexNeoIDs,
  db2jld,

  # Robot Utils
  getRangeKDEMax2D,

  # would be CloudGraphs calls
  hasBigDataElement,
  getBigDataElement,
  removeNeo4jID,

  # solver service SLAMinDB
  getcredentials,
  startSlamInDb,
  runSlamInDbOnSession,
  slamindb,
  convertdb,
  resetconvertdb,
  getmaxfactorid,

  # webserver
  SolverStatus,
  CaesarConfig,
  VisualizationConfig,

  # multisession utils
  multisessionquery,
  parsemultisessionqueryresult!,
  getLandmOtherSessNeoIDs,
  getAllLandmarkNeoIDs,
  getLocalSubGraphMultisession,
  findExistingMSConstraints,
  getprpt2kde,
  rmInstMultisessionPriors!,
  removeMultisessions!

VoidUnion{T} = Union{Void, T}

include("BearingRangeTrackingServer.jl")

include("SlamServer.jl")
include("DataUtils.jl")
include("UserFunctions.jl")

# Configuration
include("config/CaesarConfig.jl")

# using CloudGraphs
include("cloudgraphs/SolverStatus.jl")
include("cloudgraphs/CloudGraphIntegration.jl") # Work in progress code
include("cloudgraphs/ConvertGeneralSlaminDB.jl")
include("cloudgraphs/slamindb.jl")
include("cloudgraphs/MultisessionUtils.jl")
include("cloudgraphs/FoveationUtils.jl")

end


text:
"""
EXAMPLE: Generate structures using the `generate_structs()` function from the
InfrastructureSystems module.

Specifically, generate the structures required to edit input files into the required formats.
These all share the supertype `Edit`, defined in src/SLiDE.jl.
The new structures are defined in the descriptor files.
"""

# See ../src/utils/generate_structs.jl.
# The entirity of this example migrated here and is called by src/SLiDE.jl.
# to ensure that all of the necessary structs are generated
# upon compilatiion of the SLiDE module.

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
diff_backends() = [:Symbolic, :ForwardDiff, :ReverseDiff, :Zygote, :FiniteDiff]

abstract type DFunction <: Function end

## Differentiation
"""
    GetDeriv(ADmode::Symbol; kwargs...) -> Function
Returns a function which generates gradients via the method specified by `ADmode`.
`GetDeriv()` currently only available via ForwardDiff.
This outputted function has argument structure `(F::Function, x::Number) -> Number`.
"""
GetDeriv(ADmode::Symbol, args...; kwargs...) = GetDeriv(Val(ADmode), args...; kwargs...)
"""
    GetGrad(ADmode::Symbol; kwargs...) -> Function
Returns a function which generates gradients via the method specified by `ADmode`.
For available backends, see `InformationGeometry.diff_backends()`.
This outputted function has argument structure `(F::Function, x::AbstractVector) -> AbstractVector`.
"""
GetGrad(ADmode::Symbol, args...; kwargs...) = GetGrad(Val(ADmode), args...; kwargs...)
"""
    GetJac(ADmode::Symbol; kwargs...) -> Function
Returns a function which generates jacobians via the method specified by `ADmode`.
For available backends, see `InformationGeometry.diff_backends()`.
This outputted function has argument structure `(F::Function, x::AbstractVector) -> AbstractMatrix`.
"""
GetJac(ADmode::Symbol, args...; kwargs...) = GetJac(Val(ADmode), args...; kwargs...)
"""
    GetHess(ADmode::Symbol; kwargs...) -> Function
Returns a function which generates hessians via the method specified by `ADmode`.
For available backends, see `InformationGeometry.diff_backends()`.
This outputted function has argument structure `(F::Function, x::Number) -> AbstractMatrix`.
"""
GetHess(ADmode::Symbol, args...; kwargs...) = GetHess(Val(ADmode), args...; kwargs...)
"""
    GetDoubleJac(ADmode::Symbol; kwargs...) -> Function
Returns second derivatives of a vector-valued function via the method specified by `ADmode`.
For available backends, see `InformationGeometry.diff_backends()`.
This outputted function has argument structure `(F::Function, x::AbstractVector) -> AbstractArray{3}`.

THIS FEATURE IS STILL EXPERIMENTAL.
"""
GetDoubleJac(ADmode::Symbol, args...; kwargs...) = GetDoubleJac(Val(ADmode), args...; kwargs...)
"""
    GetMatrixJac(ADmode::Symbol; kwargs...) -> Function
Returns second derivatives of an array-valued function via the method specified by `ADmode`.
For available backends, see `InformationGeometry.diff_backends()`.
This outputted function has argument structure `(F::Function, x::AbstractVector) -> AbstractArray{n+1}` if `F` outputs `AbstractArray{n}`.

THIS FEATURE IS STILL EXPERIMENTAL.
"""
GetMatrixJac(ADmode::Symbol, args...; kwargs...) = GetMatrixJac(Val(ADmode), args...; kwargs...)

# Overload GetDeriv, GetGrad and so on with methods for symbolic inference!
# Currently: No symbolic passthrough for GetDeriv(ADmode), only for GetDeriv(ADmode, F). Does this result in better performance for GetDeriv(ADmode)?

GetDeriv(ADmode::Val; Kwargs...) = EvaluateDerivative(F::Function, X; kwargs...) = _GetDeriv(ADmode; Kwargs...)(F, X; kwargs...)
GetDeriv(ADmode::Val, F::DFunction; Kwargs...) = EvaldF(F)
function GetDeriv(ADmode::Val, F::Function; Kwargs...)
    EvaluateDeriv(X::Number) = _GetDeriv(ADmode; Kwargs...)(F, X)
    EvaluateDeriv(X::Num) = _GetDerivPass(F, X)
end
function GetDeriv(ADmode::Val{:Symbolic}, F::Function; kwargs...)
    M = try GetSymbolicDerivative(F, 1, :derivative; kwargs...)   catch;  nothing  end
    if isnothing(M)
        @warn "Unable to compute symbolic derivative of $F, falling back to ForwardDiff."
        GetDeriv(Val(:ForwardDiff), F)
    else M end
end

GetGrad(ADmode::Val; Kwargs...) = EvaluateGradient(F::Function, X; kwargs...) = _GetGrad(ADmode; Kwargs...)(F, X; kwargs...)
GetGrad(ADmode::Val, F::DFunction; Kwargs...) = EvaldF(F)
function GetGrad(ADmode::Val, F::Function; Kwargs...)
    EvaluateGradient(X::AbstractVector{<:Number}) = _GetGrad(ADmode; Kwargs...)(F, X)
    EvaluateGradient(X::AbstractVector{<:Num}) = _GetGradPass(F, X)
end
function GetGrad(ADmode::Val{:Symbolic}, F::Function; kwargs...)
    M = try GetSymbolicDerivative(F, :gradient; kwargs...)   catch;  nothing  end
    if isnothing(M)
        @warn "Unable to compute symbolic derivative of $F, falling back to ForwardDiff."
        GetGrad(Val(:ForwardDiff), F)
    else M end
end

GetJac(ADmode::Val; Kwargs...) = EvaluateJacobian(F::Function, X; kwargs...) = _GetJac(ADmode; Kwargs...)(F, X; kwargs...)
GetJac(ADmode::Val, F::DFunction; Kwargs...) = EvaldF(F)
function GetJac(ADmode::Val, F::Function; Kwargs...)
    EvaluateJacobian(X::AbstractVector{<:Number}) = _GetJac(ADmode; Kwargs...)(F, X)
    EvaluateJacobian(X::AbstractVector{<:Num}) = _GetJacPass(F, X)
end
function GetJac(ADmode::Val{:Symbolic}, F::Function; kwargs...)
    M = try GetSymbolicDerivative(F, :jacobian; kwargs...)   catch;  nothing  end
    if isnothing(M)
        @warn "Unable to compute symbolic derivative of $F, falling back to ForwardDiff."
        GetJac(Val(:ForwardDiff), F)
    else M end
end

GetHess(ADmode::Val; Kwargs...) = EvaluateHessian(F::Function, X; kwargs...) = _GetHess(ADmode; Kwargs...)(F, X; kwargs...)
GetHess(ADmode::Val, F::DFunction; Kwargs...) = EvalddF(F)
function GetHess(ADmode::Val, F::Function; Kwargs...)
    EvaluateHess(X::AbstractVector{<:Number}) = _GetHess(ADmode; Kwargs...)(F, X)
    EvaluateHess(X::AbstractVector{<:Num}) = _GetHessPass(F, X)
end
function GetHess(ADmode::Val{:Symbolic}, F::Function; kwargs...)
    M = try GetSymbolicDerivative(F, :hessian; kwargs...)   catch;  nothing  end
    if isnothing(M)
        @warn "Unable to compute symbolic derivative of $F, falling back to ForwardDiff."
        GetHess(Val(:ForwardDiff), F)
    else M end
end

GetMatrixJac(ADmode::Val; Kwargs...) = EvaluateMatrixJacobian(F::Function, X; kwargs...) = _GetMatrixJac(ADmode; Kwargs...)(F, X; kwargs...)
GetMatrixJac(ADmode::Val, F::DFunction; Kwargs...) = EvaldF(F)
function GetMatrixJac(ADmode::Val, F::Function; Kwargs...)
    m = GetArgLength(F);    f = size(F(ones(m)))
    EvaluateMatrixJacobian(X::AbstractVector{<:Number}) = reshape(_GetJac(ADmode)(vec∘F, X), f..., m)
    EvaluateMatrixJacobian(X::AbstractVector{<:Num}) = _GetMatrixJacPass(F, X)
end
function GetMatrixJac(ADmode::Val{:Symbolic}, F::Function; kwargs...)
    M = try GetSymbolicDerivative(F, :matrixjacobian; kwargs...)   catch;  nothing  end
    if isnothing(M)
        @warn "Unable to compute symbolic derivative of $F, falling back to ForwardDiff."
        GetMatrixJac(Val(:ForwardDiff), F)
    else M end
end
# For emergencies: needs an extra evaluation of function to determine length(Func(p))
function _GetMatrixJac(ADmode::Val; kwargs...)
    Functor(Func::Function, X::AbstractVector{<:Number}) = reshape(_GetJac(ADmode)(vec∘Func, X), size(Func(X))..., length(p))
end

GetDoubleJac(ADmode::Val; Kwargs...) = EvaluateDoubleJacobian(F::Function, X; kwargs...) = _GetDoubleJac(ADmode; Kwargs...)(F, X; kwargs...)
GetDoubleJac(ADmode::Val, F::DFunction; Kwargs...) = EvalddF(F)
function GetDoubleJac(ADmode::Val, F::Function; Kwargs...)
    m = GetArgLength(F);    f = length(F(ones(m)))
    if f == 1
        EvaluateDoubleJac(X::AbstractVector{<:Number}) = reshape(_GetJac(ADmode)(vec∘(z->_GetJac(ADmode)(F,z)), X), m, m)
        EvaluateDoubleJac(X::AbstractVector{<:Num}) = _GetDoubleJacPass(F, X)
    else
        EvaluateDoubleJacobian(X::AbstractVector{<:Number}) = reshape(_GetJac(ADmode)(vec∘(z->_GetJac(ADmode)(F,z)), X), f, m, m)
        EvaluateDoubleJacobian(X::AbstractVector{<:Num}) = _GetDoubleJacPass(F, X)
    end
end
function GetDoubleJac(ADmode::Val{:Symbolic}, F::Function; kwargs...)
    M = try GetSymbolicDerivative(F, :doublejacobian; kwargs...)   catch;  nothing  end
    if isnothing(M)
        @warn "Unable to compute symbolic derivative of $F, falling back to ForwardDiff."
        GetDoubleJac(Val(:ForwardDiff), F)
    else M end
end
# For emergencies: needs an extra evaluation of function to determine length(Func(p))
function _GetDoubleJac(ADmode::Val; kwargs...)
    Functor(Func::Function, X) = reshape(_GetJac(ADmode)(vec∘(z->_GetJac(ADmode)(Func,z)), X), length(Func(X)), length(X), length(X))
end

# Fall back to ForwardDiff as standard
_GetDeriv(ADmode::Val{true}; kwargs...) = _GetDeriv(Val(:ForwardDiff); kwargs...)
_GetGrad(ADmode::Val{true}; kwargs...) = _GetGrad(Val(:ForwardDiff); kwargs...)
_GetJac(ADmode::Val{true}; kwargs...) = _GetJac(Val(:ForwardDiff); kwargs...)
_GetHess(ADmode::Val{true}; kwargs...) = _GetHess(Val(:ForwardDiff); kwargs...)
_GetDoubleJac(ADmode::Val{true}; kwargs...) = _GetDoubleJac(Val(:ForwardDiff); kwargs...)

_GetDeriv(ADmode::Val{:ForwardDiff}; kwargs...) = ForwardDiff.derivative
_GetGrad(ADmode::Val{:ForwardDiff}; kwargs...) = ForwardDiff.gradient
_GetJac(ADmode::Val{:ForwardDiff}; kwargs...) = ForwardDiff.jacobian
_GetHess(ADmode::Val{:ForwardDiff}; kwargs...) = ForwardDiff.hessian

# Deriv not available for ReverseDiff
_GetGrad(ADmode::Val{:ReverseDiff}; kwargs...) = ReverseDiff.gradient
_GetJac(ADmode::Val{:ReverseDiff}; kwargs...) = ReverseDiff.jacobian
_GetHess(ADmode::Val{:ReverseDiff}; kwargs...) = ReverseDiff.hessian
# Deriv not available for Zygote
_GetGrad(ADmode::Val{:Zygote}; order::Int=-1, kwargs...) = (Func::Function,p;Kwargs...) -> Zygote.gradient(Func, p; kwargs...)[1]
_GetJac(ADmode::Val{:Zygote}; order::Int=-1, kwargs...) = (Func::Function,p;Kwargs...) -> Zygote.jacobian(Func, p; kwargs...)[1]
_GetHess(ADmode::Val{:Zygote}; order::Int=-1, kwargs...) = (Func::Function,p;Kwargs...) -> Zygote.hessian(Func, p; kwargs...)
# Deriv not available for FiniteDifferences
_GetGrad(ADmode::Val{:FiniteDiff}; order::Int=2, kwargs...) = (Func::Function,p;Kwargs...) -> FiniteDifferences.grad(central_fdm(order,1), Func, p; kwargs...)[1]
_GetJac(ADmode::Val{:FiniteDiff}; order::Int=2, kwargs...) = (Func::Function,p;Kwargs...) -> FiniteDifferences.jacobian(central_fdm(order,1), Func, p; kwargs...)[1]
_GetHess(ADmode::Val{:FiniteDiff}; order::Int=5, kwargs...) = (Func::Function,p;Kwargs...) -> FiniteDifferences.jacobian(central_fdm(order,1), z->FiniteDifferences.grad(central_fdm(order,1), Func, z)[1], p)[1]

## User has passed either Num or Vector{Num} to function, try to perfom symbolic passthrough
# Still need to extend this to functions F which are in-place
_GetDerivPass(F::Function, X) = SymbolicPassthrough(F(X), X, :derivative)
_GetGradPass(F::Function, X) = SymbolicPassthrough(F(X), X, :gradient)
_GetJacPass(F::Function, X) = SymbolicPassthrough(F(X), X, :jacobian)
_GetHessPass(F::Function, X) = SymbolicPassthrough(F(X), X, :hessian)

_GetDoubleJacPass(F::Function, X) = SymbolicPassthrough(F(X), X, :doublejacobian)
_GetMatrixJacPass(F::Function, X) = SymbolicPassthrough(F(X), X, :matrixjacobian)

GetGrad!(ADmode::Symbol, args...; kwargs...) = GetGrad!(Val(ADmode), args...; kwargs...)
GetJac!(ADmode::Symbol, args...; kwargs...) = GetJac!(Val(ADmode), args...; kwargs...)
GetHess!(ADmode::Symbol, args...; kwargs...) = GetHess!(Val(ADmode), args...; kwargs...)
GetMatrixJac!(ADmode::Symbol, args...; kwargs...) = GetMatrixJac!(Val(ADmode), args...; kwargs...)

# No Passthrough for these
GetGrad!(ADmode::Val, args...; kwargs...) = _GetGrad!(ADmode, args...; kwargs...)
GetJac!(ADmode::Val, args...; kwargs...) = _GetJac!(ADmode, args...; kwargs...)
GetHess!(ADmode::Val, args...; kwargs...) = _GetHess!(ADmode, args...; kwargs...)
GetMatrixJac!(ADmode::Val, args...; kwargs...) = _GetMatrixJac!(ADmode, args...; kwargs...)

function GetGrad!(ADmode::Val, F::Function; Kwargs...)
    EvaluateGradient!(Y::AbstractVector{<:Number}, X::AbstractVector{<:Number}) = _GetGrad!(ADmode; Kwargs...)(Y, F, X)
    EvaluateGradient!(Y::AbstractVector{<:Num}, X::AbstractVector{<:Num}) = _GetGradPass!(Y, F, X)
end
function GetJac!(ADmode::Val, F::Function; Kwargs...)
    EvaluateJacobian!(Y::AbstractMatrix{<:Number}, X::AbstractVector{<:Number}) = _GetJac!(ADmode; Kwargs...)(Y, F, X)
    EvaluateJacobian!(Y::AbstractMatrix{<:Num}, X::AbstractVector{<:Num}) = _GetJacPass!(Y, F, X)
end
function GetHess!(ADmode::Val, F::Function; Kwargs...)
    EvaluateHess!(Y::AbstractMatrix{<:Number}, X::AbstractVector{<:Number}) = _GetHess!(ADmode; Kwargs...)(Y, F, X)
    EvaluateHess!(Y::AbstractMatrix{<:Num}, X::AbstractVector{<:Num}) = _GetHessPass!(Y, F, X)
end

# Looks like for in-place functor jacobian! there is no difference for any array shape
function GetMatrixJac!(ADmode::Val, F::Function; Kwargs...)
    # m = GetArgLength(F);    f = size(F(ones(m)))
    EvaluateMatrixJacobian(Y::AbstractArray{<:Number}, X::AbstractVector{<:Number}) = _GetJac!(ADmode)(Y, F, X) # DELIBERATE!!!! GetJac!() recognizes output format from given Array
    EvaluateMatrixJacobian(Y::AbstractArray{<:Num}, X::AbstractVector{<:Num}) = _GetMatrixJacPass!(Y, F, X)
end

# Fall back to ForwarDiff as standard
_GetGrad!(ADmode::Val{true}; kwargs...) = _GetGrad!(Val(:ForwardDiff); kwargs...)
_GetJac!(ADmode::Val{true}; kwargs...) = _GetJac!(Val(:ForwardDiff); kwargs...)
_GetHess!(ADmode::Val{true}; kwargs...) = _GetHess!(Val(:ForwardDiff); kwargs...)

_GetGrad!(ADmode::Val{:ForwardDiff}; kwargs...) = ForwardDiff.gradient!
_GetJac!(ADmode::Val{:ForwardDiff}; kwargs...) = ForwardDiff.jacobian!
_GetHess!(ADmode::Val{:ForwardDiff}; kwargs...) = ForwardDiff.hessian!
_GetGrad!(ADmode::Val{:ReverseDiff}; kwargs...) = ReverseDiff.gradient!
_GetJac!(ADmode::Val{:ReverseDiff}; kwargs...) = ReverseDiff.jacobian!
_GetHess!(ADmode::Val{:ReverseDiff}; kwargs...) = ReverseDiff.hessian!

# Need to extend this to functions F which are themselves also in-place
_GetGradPass!(Y, F::Function, X) = (Y.=SymbolicPassthrough(F(X), X, :gradient))
_GetJacPass!(Y, F::Function, X) = (Y.=SymbolicPassthrough(F(X), X, :jacobian))
_GetHessPass!(Y, F::Function, X) = (Y.=SymbolicPassthrough(F(X), X, :hessian))
_GetMatrixJacPass!(Y, F::Function, X) = (Y.=SymbolicPassthrough(F(X), X, :matrixjacobian))

function _ConsistencyCheckOld(Fexpr, var, deriv::Symbol)
    if deriv == :jacobian || deriv == :doublejacobian
        @assert (Fexpr isa AbstractVector{<:Num}) && (var isa AbstractVector{<:Num}) "Got $deriv but Fexpr=$(typeof(Fexpr)) and argument=$(typeof(var))."
    elseif deriv == :gradient || deriv == :hessian
        @assert (Fexpr isa Num) && (var isa AbstractVector{<:Num}) "Got $deriv but Fexpr=$(typeof(Fexpr)) and argument=$(typeof(var))."
    elseif deriv == :matrixjacobian
        @assert (Fexpr isa AbstractArray{<:Num}) && (var isa AbstractVector{<:Num}) "Got $deriv but Fexpr=$(typeof(Fexpr)) and argument=$(typeof(var))."
    elseif deriv == :derivative
        @assert (Fexpr isa Num) && (var isa Num) "Got $deriv but Fexpr=$(typeof(Fexpr)) and argument=$(typeof(var))."
    else
        throw("Invalid deriv type: $deriv.")
    end
end

_ConsistencyCheck(Fexpr, var, deriv::Symbol) = _ConsistencyCheck(Fexpr, var, Val(deriv))
_ConsistencyCheck(Fexpr::AbstractVector{<:Num}, var::AbstractVector{<:Num}, ::Union{Val{:jacobian},Val{:doublejacobian}}) = nothing
_ConsistencyCheck(Fexpr::Num, var::AbstractVector{<:Num}, ::Union{Val{:gradient},Val{:hessian}}) = nothing
_ConsistencyCheck(Fexpr::Num, var::Num, ::Val{:derivative}) = nothing
_ConsistencyCheck(Fexpr::AbstractArray{<:Num}, var::AbstractVector{<:Num}, ::Val{:matrixjacobian}) = nothing
function _ConsistencyCheck(Fexpr, var, deriv::Val{T}) where T
    if T ∉ [:derivative, :gradient, :jacobian, :hessian, :doublejacobian, :matrixjacobian]
        throw("Invalid deriv type: $T.")
    else
        throw("Got :$T but Fexpr=$(typeof(Fexpr)) and arg=$(typeof(var)).")
    end
end

"""
Executes symbolic derivative as specified by `deriv::Symbol`.
"""
function SymbolicPassthrough(Fexpr::Union{<:AbstractArray{<:Num},<:Num}, var::Union{<:AbstractVector{<:Num},<:Num}, deriv::Symbol=:jacobian; simplify::Bool=true)
    _ConsistencyCheck(Fexpr, var, deriv)
    SymbolicDoubleJacobian(V::AbstractVector{<:Num}, z::AbstractVector{<:Num}; simplify::Bool=true) = SymbolicMatrixJacobian(Symbolics.jacobian(V,z),z; simplify=simplify)
    function SymbolicMatrixJacobian(M::AbstractArray{<:Num}, z::AbstractVector{<:Num}; simplify::Bool=true)
        reshape(Symbolics.jacobian(vec(M), z; simplify=simplify), size(M)..., length(z))
    end
    if deriv == :doublejacobian
        SymbolicDoubleJacobian(Fexpr, var; simplify=simplify)
    elseif deriv == :matrixjacobian
        SymbolicMatrixJacobian(Fexpr, var; simplify=simplify)
    else
        (@eval Symbolics.$deriv)(Fexpr, var; simplify=simplify)
    end
end
"""
    GetSymbolicDerivative(F::Function, inputdim::Int=GetArgLength(F), deriv::Symbol=:jacobian; timeout::Real=5, inplace::Bool=false, parallel::Bool=false)
Computes symbolic derivatives, including `:jacobian`, `:gradient`, `:hessian` and `:derivative` which are specified via `deriv`.
Special care has to be taken that the correct `inputdim` is specified! Silent errors may occur otherwise.
"""
function GetSymbolicDerivative(F::Function, inputdim::Int=GetArgLength(F), deriv::Symbol=:jacobian; timeout::Real=5, kwargs...)
    @assert deriv ∈ [:derivative, :gradient, :jacobian, :hessian, :doublejacobian, :matrixjacobian]
    @variables x X[1:inputdim]
    var = inputdim > 1 ? X : (try F(1.0); x catch; X end)
    Fexpr = KillAfter(F, var; timeout=timeout)
    # Warning already thrown in KillAfter
    isnothing(Fexpr) && return nothing
    GetSymbolicDerivative(Fexpr, var, deriv; kwargs...)
end
GetSymbolicDerivative(F::Function, deriv::Symbol; kwargs...) = GetSymbolicDerivative(F, GetArgLength(F), deriv; kwargs...)

function GetSymbolicDerivative(Fexpr::Union{<:AbstractArray{<:Num},<:Num}, var::Union{<:AbstractVector{<:Num},<:Num}, deriv::Symbol=:jacobian; simplify::Bool=true, inplace::Bool=false, parallel::Bool=false, kwargs...)
    derivative = SymbolicPassthrough(Fexpr, var, deriv; simplify=simplify)
    Builder(derivative, var; parallel=parallel, inplace=inplace, kwargs...)
end


text:
# %% Color theme definitions
struct CyclicContainer <: AbstractVector{String}
    c::Vector{String}
    n::Int
end
CyclicContainer(c) = CyclicContainer(c, 0)

Base.length(c::CyclicContainer) = length(c.c)
Base.size(c::CyclicContainer) = size(c.c)
Base.iterate(c::CyclicContainer, state=1) = Base.iterate(c.c, state)
Base.getindex(c::CyclicContainer, i) = c.c[(i-1)%length(c.c) + 1]
Base.getindex(c::CyclicContainer, i::AbstractArray) = c.c[i]
function Base.getindex(c::CyclicContainer)
    c.n += 1
    c[c.n]
end
Base.iterate(c::CyclicContainer, i = 1) = iterate(c.c, i)

COLORSCHEME = [
    "#6D44D0",
    "#2CB3BF",
    "#1B1B1B",
    "#DA5210",
    "#03502A",
    "#866373",
]

COLORS = CyclicContainer(COLORSCHEME)
LINESTYLES = CyclicContainer(["-", ":", "--", "-."])

# %% Makie styling
# other styling elements for Makie
set_theme!(;
    palette = (color = COLORSCHEME,),
    fontsize = 22,
    figure_padding = 1,
    resolution = (1000, 500),
    linewidth = 3.0,
)


text:
# Adapted from https://github.com/xofbd/pandas_vs_PostgreSQL/blob/master/create_dataset.py
# (which was released under MIT license: Copyright (c) 2018 Don Bruce Fox)
# See also: https://blog.thedataincubator.com/2018/01/pandas-vs-postgresql/

using MinimumViableTables

using Random

function create_dataset(n = 1000)
    Random.srand(1)

    id = randperm(n)
    section = rand(['A', 'B', 'C', 'D'], n)
    score_1 = rand(n)
    score_2 = rand(n)
    score_3 = rand(n)

    t_a = Table(id = id, section = section, score_1 = score_1, score_2 = score_2)
    t_b = Table(id2 = id, score_3 = score_3)

    return (t_a, t_b)
end


text:
"""
compute_scc(m::Model=get_model(); year::Int = nothing, last_year::Int = 2595), prtp::Float64 = 0.03)

Computes the social cost of CO2 for an emissions pulse in `year` for the provided MimiRICE2010 model.
If no model is provided, the default model from MimiRICE2010.get_model() is used.
Constant discounting is used from the specified pure rate of time preference `prtp`.
"""
function compute_scc(m::Model=get_model(); year::Union{Int, Nothing} = nothing, last_year::Int = model_years[end], prtp::Float64 = 0.015, eta::Float64=1.5)
    year === nothing ? error("Must specify an emission year. Try `compute_scc(m, year=2015)`.") : nothing
    !(last_year in model_years) ? error("Invlaid value of $last_year for last_year. last_year must be within the model's time index $model_years.") : nothing
    !(year in model_years[1]:10:last_year) ? error("Cannot compute the scc for year $year, year must be within the model's time index $(model_years[1]):10:$last_year.") : nothing

    mm = get_marginal_model(m; year = year)

    return _compute_scc(mm, year=year, last_year=last_year, prtp=prtp, eta=eta)
end

"""
compute_scc_mm(m::Model=get_model(); year::Int = nothing, last_year::Int = 2595, prtp::Float64 = 0.03)

Returns a NamedTuple (scc=scc, mm=mm) of the social cost of carbon and the MarginalModel used to compute it.
Computes the social cost of CO2 for an emissions pulse in `year` for the provided MimiRICE2010 model.
If no model is provided, the default model from MimiRICE2010.get_model() is used.
Constant discounting is used from the specified pure rate of time preference `prtp`.
"""
function compute_scc_mm(m::Model=get_model(); year::Union{Int, Nothing} = nothing, last_year::Int = model_years[end], prtp::Float64 = 0.015, eta::Float64=1.5)
    year === nothing ? error("Must specify an emission year. Try `compute_scc_mm(m, year=2015)`.") : nothing
    !(last_year in model_years) ? error("Invlaid value of $last_year for last_year. last_year must be within the model's time index $model_years.") : nothing
    !(year in model_years[1]:10:last_year) ? error("Cannot compute the scc for year $year, year must be within the model's time index $(model_years[1]):10:$last_year.") : nothing

    mm = get_marginal_model(m; year = year)
    scc = _compute_scc(mm; year=year, last_year=last_year, prtp=prtp, eta=eta)

    return (scc = scc, mm = mm)
end

# helper function for computing SCC from a MarginalModel, not to be exported or advertised to users
function _compute_scc(mm::MarginalModel; year::Int, last_year::Int, prtp::Float64, eta::Float64)
    ntimesteps = findfirst(isequal(last_year), model_years)     # Will run through the timestep of the specified last_year
    run(mm, ntimesteps=ntimesteps)

    marginal_damages = -1 * mm[:neteconomy, :C][1:ntimesteps, :] * 10.0^12 * 12/44 # convert from trillion $/ton C to $/ton CO2; multiply by -1 to get positive value for damages
    global_marginal_damages = dropdims(sum(marginal_damages, dims = 2), dims=2)

    global_c = dropdims(sum(mm.base[:neteconomy, :C], dims = 2), dims=2)
    global_pop = dropdims(sum(mm.base[:neteconomy, :l], dims = 2), dims=2)
    global_cpc = 1000 .* global_c ./ global_pop

    year_index = findfirst(isequal(year), model_years)

    df = [zeros(year_index-1)..., ((global_cpc[year_index]/global_cpc[i])^eta * 1/(1+prtp)^(t-year) for (i,t) in enumerate(model_years) if year<=t<=last_year)...]
    scc = sum(df .* marginal_damages * 10)  # currently implemented as a 10year step function; so each timestep of discounted marginal damages is multiplied by 10
    return scc
end

"""
get_marginal_model(m::Model = get_model(); year::Int = nothing)

Creates a Mimi MarginalModel where the provided m is the base model, and the marginal model has additional emissions of CO2 in year `year`.
If no Model m is provided, the default model from MimiDICE2010.get_model() is used as the base model.
"""
function get_marginal_model(m::Model=get_model(); year::Union{Int, Nothing} = nothing)
    year === nothing ? error("Must specify an emission year. Try `get_marginal_model(m, year=2015)`.") : nothing
    !(year in model_years) ? error("Cannot add marginal emissions in $year, year must be within the model's time index $(model_years[1]):10:$last_year.") : nothing

    mm = create_marginal_model(m, 1e10) # Pulse has a value of 1GtC per year for ten years
    add_marginal_emissions!(mm.modified, year)

    return mm
end

"""
Adds a marginal emission component to year m which adds 1Gt of additional C emissions per year for ten years starting in the specified `year`.
"""
function add_marginal_emissions!(m::Model, year::Int)
    add_comp!(m, Mimi.adder, :marginalemission, before=:co2cycle)

    time = Mimi.dimension(m, :time)
    addem = zeros(length(time))
    addem[time[year]] = 1.0     # 1 GtC per year for ten years

    set_param!(m, :marginalemission, :add, addem)
    connect_param!(m, :marginalemission, :input, :emissions, :E)
    connect_param!(m, :co2cycle, :E, :marginalemission, :output)
end

# OLD previously used marginal damages functions

function getmarginal_rice_models(;emissionyear=2005,datafile=joinpath(@__DIR__, "..", "data", "RICE_2010_base_000.xlsm"))

    RICE = get_model()
    run(RICE)

    mm = MarginalModel(RICE)
    m1 = mm.base
    m2 = mm.modified

    add_comp!(m2, Mimi.adder, :marginalemission, before=:co2cycle)

    time = Mimi.dimension(m1, :time)
    addem = zeros(length(time))
    addem[time[emissionyear]] = 1.0

    set_param!(m2,:marginalemission,:add,addem)
    connect_param!(m2,:marginalemission,:input,:emissions,:E)
    connect_param!(m2, :co2cycle,:E,:marginalemission,:output)

    run(m1)
    run(m2)

    return m1, m2
end

# This function returns a matrix of marginal damages per one t of carbon emission in the
# emissionyear parameter year.
function getmarginaldamages_rice(;emissionyear=2005,datafile=joinpath(@__DIR__, "..", "data", "RICE_2010_base_000.xlsm"))
    m1, m2 = getmarginal_rice_models(emissionyear=emissionyear, datafile=datafile)

    run(m1)
    run(m2)

    damage1 = m1[:damages,:DAMAGES]
    damage2 = m2[:damages,:DAMAGES]

    # Calculate the marginal damage between run 1 and 2 for each
    # year/region
    marginaldamage = (damage2.-damage1) .* 1e12 / 1e9 / 10.

    return marginaldamage
end


text:
"""
# module AdventOfCode

- Julia version:
- Author: Paul.Mealor
- Date: 2020-12-02

# Examples

```jldoctest
julia>
```
"""
module AoC

    export lines, exampleLines
    export firstLine, firstExampleLine
    export charMap, exampleCharMap
    export intMap, exampleIntMap
    export ints, exampleInts
    export parseLine

    lines(day) = open(readlines, "src/day" * string(day) * "-input.txt")
    exampleLines(day, n) = open(readlines, "src/day" * string(day) * "-example-" * string(n) * ".txt")

    firstLine(day) = open(readline, "src/day" * string(day) * "-input.txt")
    firstExampleLine(day, n) = open(readline, "src/day" * string(day) * "-example-" * string(n) * ".txt")

    linesToCharMap(lines) = reshape([c for line in lines for c in line], length(lines[1]), length(lines))
    charMap(day) = lines(day) |> linesToCharMap
    exampleCharMap(day, n) = exampleLines(day, n) |> linesToCharMap

    intMap(day) = charMap(day) .|> c -> parse(Int, c)
    exampleIntMap(day, n) = exampleCharMap(day, n) .|> c -> parse(Int, c)

    ints(day) = lines(day) |> ll -> parse.(Int, ll)
    exampleInts(day, n) = exampleLines(day, n) |> ll -> parse.(Int, ll)

    parse(::Type{T}, string::AbstractString) where T<:AbstractString = string
    parse(type, string) = Base.parse(type, string)
    parseLine(line, types...) = split(line, r"\s+"; limit=length(types)) |> tt -> parse.(types, tt)
end


text:
module WindTrigger

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Wind Pattern" => Ahorn.EntityPlacement(
        Maple.WindTrigger,
        "rectangle"
    )
)

function Ahorn.editingOptions(trigger::Maple.WindTrigger)
    return Dict{String, Any}(
        "pattern" => Maple.wind_patterns
    )
end

end


text:
using Test

include("all_examples.jl")

@testset "run_examples.jl" begin
    @testset "$(example)" for example in EXAMPLES
        include(example)
    end
end


text:
using LinearAlgebra, Statistics

function nelder_mead(f, S, ϵ; α=1.0, β=2.0, γ=0.5)
    Δ, y_arr = Inf, f.(S)
    iteractions = 0
    while Δ > ϵ
        p = sortperm(y_arr) # sort lowest to highest
        S, y_arr = S[p], y_arr[p]
        xl, yl = S[1], y_arr[1] # lowest
        xh, yh = S[end], y_arr[end] # highest
        xs, ys = S[end-1], y_arr[end-1] # second-highest
        xm = mean(S[1:end-1]) # centroid
        xr = xm + α*(xm - xh) # reflection point
        yr = f(xr)

        if yr < yl
            xe = xm + β*(xr-xm) # expansion point
            ye = f(xe)
            S[end],y_arr[end] = ye < yr ? (xe, ye) : (xr, yr)
        elseif yr > ys
            if yr ≤ yh
                xh, yh, S[end], y_arr[end] = xr, yr, xr, yr
            end
            xc = xm + γ*(xh - xm) # contraction point
            yc = f(xc)
            if yc > yh
                for i in 2 : length(y_arr)
                    S[i] = (S[i] + xl)/2
                    y_arr[i] = f(S[i])
                end
            else
                S[end], y_arr[end] = xc, yc
            end
        else
            S[end], y_arr[end] = xr, yr
        end

        Δ = std(y_arr, corrected=false)
        iteractions = iteractions + 1
    end
    println("Iteractions: ",iteractions)
    return S,iteractions
end

#valores dos parametros
ka = 9.0
kb = 2.0
La = 10.0
Lb = 10.0

#Caso 1
#F1 = 2.0; F2 = 4.0;

#Caso 2
#F1 = 0.0 ;F2 = 4.0;

#Caso 3
F1 = 0.0; F2 = 400.0;

#tolerancia
tol = 1e-8

#ponto inicial
u0 = [-0.01,0.0]
#u0 = [0.0,0.0]
#numero de pontos iniciais = dimensao do problema + 1
u1 = [u0[1]-0.00001,u0[2]+0.00001]
u2 = [u0[1]+0.00001,u0[2]+0.00001]
u = [u0,u1,u2]

#function
f(u) = 0.5e0 * ka * (sqrt(u[1] ^ 2 + (La - u[2]) ^ 2) - La) ^ 2 + 0.5e0 * kb * (sqrt(u[1] ^ 2 + (Lb + u[2]) ^ 2) - Lb) ^ 2 - F1 * u[1] - F2 * u[2];

#tolerancia
tol = 1e-8

@time u_est,iterac = nelder_mead(f, u, tol)
println("u*: ",mean(u_est))
println("f(u*): ",f(mean(u_est)))


text:
"""
    function get_Pg_thermal_slack(args ...)

Get Slack Generator, locate their buses and returns a matrix with their Name, Power output
and Bus number (sorted by Bus).

# Arguments
- `sys::System`:                            Power system in p.u. (from PowerSystems.jl)
- `res::PowerSimulations.
    OperationsProblemResults`:              Results of the solved OPF for the system
                                             (from PowerSimulations.jl)

"""
function get_Pg_thermal_slack(sys::System, res::PowerSimulations.OperationsProblemResults)

    # Define elements
    gen_loc = 0
    ngen_slack = 0 #number of gen slacks
    gen_thermal_slack_name = "Empty"
    gen_thermal_slack_busnumber = 0
    gens_thermal_slack_Pg = Array{Float64}(undef, 0)
    gens_thermal_slack_names = Array{String}(undef, 0)
    gens_thermal_slack_busnumbers = Array{Int64}(undef, 0)

    # Get all thermal Generators
    gens_thermal = get_components(ThermalStandard, sys)
    all_PGenThermal = get_variables(res)[:P__ThermalStandard] #Optimised PGen

    # Identify Slack Generators and save their information
    for gen_thermal in gens_thermal
        gen_loc = gen_loc +1
        if get_bustype(gen_thermal.bus)== BusTypes.REF
            ngen_slack = ngen_slack + 1;
            gen_thermal_slack_name = gen_thermal.name
            resize!(gens_thermal_slack_Pg, ngen_slack)
            resize!(gens_thermal_slack_names, ngen_slack)
            resize!(gens_thermal_slack_busnumbers, ngen_slack)
            gens_thermal_slack_Pg[ngen_slack] = all_PGenThermal[1, gen_thermal_slack_name]
            gens_thermal_slack_names[ngen_slack] = gen_thermal_slack_name
            gens_thermal_slack_busnumbers[ngen_slack] = gen_thermal.bus.number
        end
    end

    # Matrix of generators with their information sorted
    gens_thermal_slack_Matrix = Array{Any}(undef,(ngen_slack, 3))
    gens_thermal_slack_Matrix[:, 1] = gens_thermal_slack_names
    gens_thermal_slack_Matrix[:, 2] = gens_thermal_slack_Pg
    gens_thermal_slack_Matrix[:, 3] = gens_thermal_slack_busnumbers
    # sort by bus to match with PTDF structure
    gens_thermal_slack_Matrix[sortperm(gens_thermal_slack_Matrix[:, 3]), :]

    return gens_thermal_slack_Matrix
end

max_stars_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_issues_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_forks_repo_name: invenia/MaximumGeneratorProfitTCRDD


text:
function parseall(str)
    pos = firstindex(str)
    exs = []
    while pos <= lastindex(str)
        ex, pos = Meta.parse(str, pos)
        push!(exs, ex)
    end
    if length(exs) == 0
        throw(Meta.ParseError("end of input"))
    elseif length(exs) == 1
        return exs[1]
    else
        return Expr(:block, exs...)
    end
end

"""
    fix_makefile(makefile)

Takes in the path to a Documenter.jl-compatible `make.jl` file and
1. Removes all calls to `deploydocs`.
2. Records the custom `build` kwarg in the `makedocs` call.
3. Makes sure that `format` is set to `Documenter.HTML(...)` in the `makedocs` call.

Return a tuple of `new_make_expr, buildpath`.
"""
function fix_makefile(makefile, documenter_version = v"0.24")
    # default output path:
    buildpath = joinpath(dirname(makefile), "build")
    should_break = false

    ast = parseall(read(makefile, String))
    make_expr = Expr(:block)

    for elem in ast.args
        # skip deploydocs
        Meta.isexpr(elem, :call) && elem.args[1] == :deploydocs && continue
        if Meta.isexpr(elem, :call) && elem.args[1] == :makedocs
            should_break = true

            # rewrite makedoc call to respect our requirements
            new_args = []

            has_fmt = false
            has_sitename = false
            has_linkcheck = false
            has_doctest = false
            html = documenter_version < v"0.21" ? QuoteNode(:html) : :(Documenter.HTML())

            fixkwarg = argument -> begin
                if Meta.isexpr(argument, :kw)
                    name, arg = argument.args
                    # assure that we generate HTML
                    if name == :format
                        has_fmt = true
                        if Meta.isexpr(arg, :call) && arg.args[1] == :(Documenter.HTML) && !(html isa QuoteNode)
                            append!(html.args, arg.args[2:end])
                        end
                        argument.args[2] = html
                    end
                    # filter out root + build dir
                    if name == :build
                        # record custom output dir if we can infer it statically, otherwise just use `build`
                        if arg isa String
                            buildpath = joinpath(dirname(makefile), arg)
                        else
                            argument.args[2] = "build"
                            buildpath = joinpath(dirname(makefile), "build")
                        end
                    end
                    if name == :sitename
                        has_sitename = true
                    end
                    if name == :linkcheck
                        has_linkcheck = true
                        argument.args[2] = false
                    end
                    if name == :doctest
                        has_doctest = true
                        argument.args[2] = false
                    end
                end
            end

            for argument in elem.args
                if Meta.isexpr(argument, :parameters)
                    for argument in argument.args
                        fixkwarg(argument)
                    end
                else
                    fixkwarg(argument)
                end

                push!(new_args, argument)
            end

            if !has_sitename && documenter_version < v"0.20"
                name = basename(dirname(dirname(dirname(makefile))))
                push!(new_args, Expr(:kw, :sitename, name))
            end

            if !has_fmt
                push!(new_args, Expr(:kw, :format, html))
            end

            if !has_linkcheck
                push!(new_args, Expr(:kw, :linkcheck, false))
            end

            if !has_doctest
                push!(new_args, Expr(:kw, :doctest, false))
            end

            # make sure to overwrite `root`:
            push!(new_args,
                Expr(:kw, :root, dirname(makefile)),
            )
            elem = Expr(:call, new_args...)
        end

        push!(make_expr.args, elem)

        # ignore everything after `makedocs` call
        should_break && break
    end

    # it's possible a make.jl-file doesn't have `using Documenter` at all, so let's just unconditionally add that here
    pushfirst!(make_expr.args, Expr(:using, Expr(:., :Documenter)))

    fix_lnns(make_expr, makefile)

    return make_expr, buildpath
end

function fix_lnns(expr, filepath)
    for (i, el) in enumerate(expr.args)
        if el isa LineNumberNode
            expr.args[i] = LineNumberNode(el.line, Symbol(filepath))
        elseif el isa Expr
            fix_lnns(el, filepath)
        end
    end
end


text:
using Base.Iterators: partition

@testset "ConjugateBernoulli" begin
    # underlying distribution
    n = 100000
    θ = 0.55

    Random.seed!(1234)
    data = rand(Bernoulli(θ), n)

    @testset "Model update and posteriros sampling" begin
        model = ConjugateBernoulli(1, 1)
        stats = BetaStatistics(s=sum(data), n=n)
        update!(model, stats)
        posterior_samples = samplepost(model, 100_000)
        @test mean(posterior_samples.θ) ≈ 0.55
    end

    @testset "BetaStatistics constructor" begin
        stats = BetaStatistics(data)
        @test stats.n == n
        @test stats.s / stats.n ≈ mean(data)
    end

    @testset "Online update vs. One shot update" begin
        numsteps = 10
        α = 10
        β = 30
        data_by_step = reshape(data, (10, div(n, numsteps)))
        model_online = ConjugateBernoulli(α, β)
        for i = 1:numsteps
            stats = BetaStatistics(data_by_step[i,:])
            update!(model_online, stats)
        end

        model_oneshot = ConjugateBernoulli(α, β)
        stats = BetaStatistics(data)
        update!(model_oneshot, stats)

        @test model_online.dist == model_oneshot.dist
    end

end

@testset "ConjugateExponential" begin
    # underlying distribution
    n = 100_000
    Random.seed!(1234)
    θ = 25.0
    data = rand(Exponential(θ), n)

    @testset "Posteriors sampling" begin
        model = ConjugateExponential(0.001, 1000)
        stats = GammaStatistics(n=n, x̄=mean(data))
        update!(model, stats)
        posterior_samples = samplepost(model, 100_000)
        @test isapprox(mean(posterior_samples.θ), θ, atol=0.1)
    end

    @testset "GammaStatistics constructor" begin
        stats = GammaStatistics(data)
        @test stats.n == n
        @test isapprox(stats.x̄, 25, atol=0.1)
    end

    @testset "Online update vs. One shot update" begin
        numsteps = 10
        α = 2
        θ = 5
        data_by_step = reshape(data, (10, div(n, numsteps)))
        model_online = ConjugateExponential(α, θ)
        for i = 1:numsteps
            stats = GammaStatistics(data_by_step[i, :])
            update!(model_online, stats)
        end

        model_oneshot = ConjugateExponential(α, θ)
        stats = GammaStatistics(data)
        update!(model_oneshot, stats)

        @test typeof(model_online.dist) === typeof(model_oneshot.dist)
        @test model_online.dist.α ≈ model_oneshot.dist.α
        @test model_online.dist.θ ≈ model_oneshot.dist.θ
    end
end

@testset "ConjugatePoisson" begin
    # underlying distribution
    n = 100_000
    Random.seed!(1234)
    λ = 0.5
    data = rand(Poisson(λ), n)

    @testset "Posteriors sampling" begin
        model = ConjugatePoisson(0.001, 1000)
        stats = GammaStatistics(n=n, x̄=mean(data))
        update!(model, stats)
        posterior_samples = samplepost(model, 100_000)
        @test isapprox(mean(posterior_samples.λ), λ, atol=0.1)
    end

    @testset "GammaStatistics constructor" begin
        stats = GammaStatistics(data)
        @test stats.n == n
        @test isapprox(stats.x̄, λ, atol=0.1)
    end

    @testset "Online update vs. One shot update" begin
        numsteps = 10
        α = 2
        θ = 5
        data_by_step = reshape(data, (10, div(n, numsteps)))
        model_online = ConjugatePoisson(α, θ)
        for i = 1:numsteps
            stats = GammaStatistics(data_by_step[i, :])
            update!(model_online, stats)
        end

        model_oneshot = ConjugatePoisson(α, θ)
        stats = GammaStatistics(data)
        update!(model_oneshot, stats)

        @test typeof(model_online.dist) === typeof(model_oneshot.dist)
        @test model_online.dist.α ≈ model_oneshot.dist.α
        @test model_online.dist.θ ≈ model_oneshot.dist.θ
    end
end

@testset "ConjugateNormal" begin
    # underlying distribution
    n = 100_000
    Random.seed!(1234)
    μ = 13.0
    σ = 2.5
    data = rand(Normal(μ, σ), n)

    @testset "Posteriors sampling" begin
        model = ConjugateNormal(0.0, 1.0, 0.001, 0.001)
        stats = NormalStatistics(data)
        update!(model, stats)
        posterior_samples = samplepost(model, 100_000)
        @test mean(posterior_samples.μ) ≈ μ
        @test std(posterior_samples.μ) ≈ σ / sqrt(n)
    end

    @testset "NormalStatistics constructor" begin
        stats = NormalStatistics(data)
        @test stats.n == n
        @test isapprox(stats.meanx, μ, atol=0.01)
        @test isapprox(stats.sdx, σ, atol=0.01)
    end

    @testset "Online update vs. One shot update" begin
        numsteps = 10
        data_by_step = reshape(data, (10, div(n, numsteps)))
        model_online = ConjugateNormal(0.0, 10.0, 10., 5.0)
        for i = 1:numsteps
            stats = NormalStatistics(data_by_step[i, :])
            update!(model_online, stats)
        end

        model_oneshot = ConjugateNormal(0.0, 10.0, 10., 5.0)
        stats = NormalStatistics(data)
        update!(model_oneshot, stats)

        @test typeof(model_online.dist) === typeof(model_oneshot.dist)
        @test model_online.dist.μ ≈ model_oneshot.dist.μ
        @test model_online.dist.v ≈ model_oneshot.dist.v
        @test model_online.dist.α ≈ model_oneshot.dist.α
        @test model_online.dist.θ ≈ model_oneshot.dist.θ
    end
end

@testset "ConjugateLogNormal" begin
    # underlying distribution
    meanlogx = 7.0
    sdlogx = 0.9
    varlogx = sdlogx * sdlogx
    n = 100000

    Random.seed!(1234)
    data  = rand(LogNormal(meanlogx, sdlogx), n)

    @testset "Posteriros sampling" begin
        model = ConjugateLogNormal(0.0, 1.0, 0.001, 0.001)
        stats = LogNormalStatistics(n=n,
            meanlogx=mean(log.(data)), sdlogx=std(log.(data)))
        update!(model, stats)

        @test model.dist.μ ≈ 7.0
        @test model.dist.v ≈ 0
        @test isapprox(model.dist.α, 50000, rtol=0.01)
        @test isapprox(model.dist.θ, 40525, rtol=0.01)

        posterior_samples = samplepost(model, 100_000)
        @test mean(posterior_samples.μ_logx) ≈ meanlogx
        @test mean(posterior_samples.σ²_logx) ≈ varlogx
        @test isapprox(mean(posterior_samples.μ_x), 1645, rtol=0.01)
        @test isapprox(mean(posterior_samples.σ²_x), 3.4e6, rtol=0.01)
    end

    @testset "LogNormalStatistics constructor" begin
        stats = LogNormalStatistics(data)
        @test stats.n == n
        @test stats.meanlogx ≈ meanlogx
        @test stats.sdlogx ≈ sdlogx
    end

    @testset "Online update vs. One shot update" begin
        numsteps = 10
        data_by_step = reshape(data, (10, div(n, numsteps)))
        model_online = ConjugateLogNormal(0.0, 10.0, 10., 5.0)
        for i = 1:numsteps
            stats = LogNormalStatistics(data_by_step[i, :])
            update!(model_online, stats)
        end

        model_oneshot = ConjugateLogNormal(0.0, 10.0, 10., 5.0)
        stats = LogNormalStatistics(data)
        update!(model_oneshot, stats)

        @test typeof(model_online.dist) === typeof(model_oneshot.dist)
        @test model_online.dist.μ ≈ model_oneshot.dist.μ
        @test model_online.dist.v ≈ model_oneshot.dist.v
        @test model_online.dist.α ≈ model_oneshot.dist.α
        @test model_online.dist.θ ≈ model_oneshot.dist.θ
    end
end


text:
using Plots
using Random
Random.seed!(1)

include_ddp()

# Model
include_model("double_integrator")

function f(model::DoubleIntegratorContinuous, x, u, w)
    [x[2]; (1.0 + w[1]) * u[1]]
end

model = DoubleIntegratorContinuous{Midpoint, FixedTime}(2, 1, 1)
n = model.n
m = model.m

struct FOHModel{I, T} <: Model{I, T}
	n::Vector{Int}
	m::Vector{Int}
	d::Vector{Int}

	model::Model
end

function foh_model(model::Model, T)
	n = [t == 1 ? model.n : model.n + model.m for t = 1:T]
	m = [t == 1 ? 2 * model.m : model.m for t = 1:T-1]
	d = [model.d for t = 1:T-1]

	FOHModel{typeof(model).parameters...}(n, m, d, model)
end

function fd(models::FOHModel, x, u, w, h, t)
	if t == 1
		return [fd(models.model,
					view(x, 1:model.n),
					view(u, 1:model.m),
					w, h, t);
				 view(u, 1:model.m) + h * view(u, model.m .+ (1:model.m))]
	else
		return [fd(models.model,
					view(x, 1:model.n),
					view(x, model.n .+ (1:model.m)),
					w, h, t);
				 view(x, model.n .+ (1:model.m)) + h * view(u, 1:model.m)]
	end
end

# Time
T = 11
h = 0.1

# FOH model
model_foh = foh_model(model, T)

# Initial conditions, controls, disturbances
x1 = [1.0; 0.0]
x_ref = [t == 1 ? zeros(model.n) : zeros(model.n + model.m) for t = 1:T]
xT = zeros(model.n + model.m)
ū = [t == 1 ? 1.0 * rand(2 * model.m) : 1.0 * rand(model.m) for t = 1:T-1]
w = [zeros(model.d) for t = 1:T-1]

# Rollout
x̄ = rollout(model_foh, x1, ū, w, h, T)

# Objective
Q = [(t > 1 ?
	 Diagonal([1.0; 1.0; 1.0e-3])
		: Diagonal([1.0; 1.0])) for t = 1:T]
q = [-2.0 * Q[t] * x_ref[t] for t = 1:T]
R = [t == 1 ? Diagonal(1.0e-3 * ones(2 * model.m)) : Diagonal(1.0e-3 * ones(model.m)) for t = 1:T-1]
r = [t == 1 ? zeros(2 * model.m) : zeros(model.m) for t = 1:T-1]

obj = StageCosts([QuadraticCost(Q[t], q[t],
	t < T ? R[t] : nothing, t < T ? r[t] : nothing) for t = 1:T], T)

function g(obj::StageCosts, x, u, t)
	T = obj.T
    if t < T
		Q = obj.cost[t].Q
		q = obj.cost[t].q
	    R = obj.cost[t].R
		r = obj.cost[t].r
        return x' * Q * x + q' * x + u' * R * u + r' * u
    elseif t == T
		Q = obj.cost[T].Q
		q = obj.cost[T].q
        return x' * Q * x + q' * x
    else
        return 0.0
    end
end

# Constraints
# ul = [-10.0]
# uu = [10.0]
p = [t < T ? 0 : model.n for t = 1:T]
info_t = Dict()#:ul => ul, :uu => uu, :inequality => (1:2 * m))
info_T = Dict(:xT => xT[1:model.n])
con_set = [StageConstraint(p[t], t < T ? info_t : info_T) for t = 1:T]

function c!(c, cons::StageConstraints, x, u, t)
	T = cons.T
	p = cons.con[t].p

	if t < T
		# ul = cons.con[t].info[:ul]
		# uu = cons.con[t].info[:uu]
		# c .= [ul - u; u - uu]
	else
		c .= view(x, 1:model.n) - cons.con[T].info[:xT]
	end
end

prob = problem_data(model_foh, obj, con_set, copy(x̄), copy(ū), w, h, T,
	n = model_foh.n, m = model_foh.m, d = model_foh.d)

# Solve
@time constrained_ddp_solve!(prob,
    max_iter = 100, verbose = true)

x, u = current_trajectory(prob)
x̄, ū = nominal_trajectory(prob)
hcat([x[1:model.n] for t = 1:T]...)

# Visualize
using Plots
plot(hcat([xT[1:model.n] for t = 1:T]...)',
    width = 2.0, color = :black, label = "")
plot!(hcat([x[t][1:model.n] for t = 1:T]...)', color = :magenta, label = "")

u_foh = [u[1][1:model.m], [x[t][model.n .+ (1:model.m)] for t = 2:T]...]
plot(hcat(u_foh...)', linetype = :steppost)

plot(hcat(u..., u[end])', linetype = :steppost)
# Simulate policy
include(joinpath(@__DIR__, "simulate.jl"))

# Model
model_sim = DoubleIntegratorContinuous{RK3, FixedTime}(model.n, model.m, model.d)
x1_sim = copy(x1)
T_sim = 10 * T

# Time
tf = h * (T - 1)
t = range(0, stop = tf, length = T)
t_sim = range(0, stop = tf, length = T_sim)
dt_sim = tf / (T_sim - 1)

# Policy
# K = [K for K in prob.p_data.K]
# plot(vcat(K...))
# # K = [prob.p_data.K[t] for t = 1:T-1]
# K, _ = tvlqr(model, x̄, ū, h, Q, R)
# # K = [-k for k in K]
# K = [-K[t] for t = 1:T-1]
# plot(vcat(K...))

# Simulate
N_sim = 100
x_sim = []
u_sim = []
J_sim = []
Random.seed!(1)
for k = 1:N_sim
	wi_sim = min(0.1, max(-0.1, 1.0e-1 * randn(1)[1]))
	w_sim = [wi_sim for t = 1:T-1]
	println("sim: $k - w = $(wi_sim[1])")

	x_ddp, u_ddp, J_ddp, Jx_ddp, Ju_ddp = simulate_linear_feedback(
		model_sim,
		K,
	    x̄, ū,
		x_ref, u_ref,
		Q, R,
		T_sim, h,
		x1_sim,
		w_sim,
		ul = ul,
		uu = uu)

	push!(x_sim, x_ddp)
	push!(u_sim, u_ddp)
	push!(J_sim, J_ddp)
end

# Visualize
idx = (1:2)
plt = plot(t, hcat(x̄...)[idx, :]',
	width = 2.0, color = :black, label = "",
	xlabel = "time (s)", ylabel = "state",
	title = "double integrator (J_avg = $(round(mean(J_sim), digits = 3)), N_sim = $N_sim)")

for xs in x_sim
	plt = plot!(t_sim, hcat(xs...)[idx, :]',
	    width = 1.0, color = :magenta, label = "")
end
display(plt)

plt = plot(t, hcat(ū..., ū[end])',
	width = 2.0, color = :black, label = "",
	xlabel = "time (s)", ylabel = "control",
	title = "double integrator (J_avg = $(round(mean(J_sim), digits = 3)), N_sim = $N_sim)")

for us in u_sim
	plt = plot!(t_sim, hcat(us..., us[end])',
		width = 1.0, color = :magenta, label = "",
		linetype = :steppost)
end
display(plt)
# u_sim
# plot(vcat(K...))


text:
# This file was generated, do not modify it. # hide
data_nmiss = dropmissing(data, :commissioning_year);


text:
const USE_GPU = false
const VISUALIZE = false
const BENCHMARK = true
include("elastic_wave_3D.jl")

# Grid sizes for the scaling experiment
grid_size = 16 * 2 .^ (1:5)

# Array allocation
T_effs = zeros(5)

# Scaling experiment loop
for it = 1:5
    T_effs[it] = elastic_wave_3D(grid_size[it], 4096.0 / 8^it)
end

# Visualization of Results
gr()
max_bandwith = 47.68 # Theoretical Max Memory Bandwidth of Ryzen 5 5600G
num_threads = Threads.nthreads()

plot(grid_size, max_bandwith .* ones(5), label="Theoretical Max Bandwidth", lw=3)
plot!(grid_size, T_effs, xlabel="grid_size", ylabel="T_eff GB/s", ylims=(0.0, 60.0), title="CPU Scaling Experiment $(num_threads) Thread(s)", label="elastic_wave_3D.jl", lw=3)
png("./docs/img/elastic_wave_3D_scaling_experiment_cpu_$(num_threads)threads")