text:
using DrWatson
quickactivate(@__DIR__)
@show projectdir()
#-----------------------------------------------------------------------------
# load packages
@info "load packages"
using ArgParse
using CovidSurvey
using Random, Dates, Turing
using Base.Threads
using StatsPlots
using PrettyTables
using DataFrames
using BSON

using Memoization
using ReverseDiff
Turing.setrdcache(true)
setadbackend(:reversediff)

plotlyjs()
@info "number of threads available: $(nthreads())"
#-----------------------------------------------------------------------------
ps = (
    warmup = 10,
    model = "cases",
    steps  = 10,
    seed   = nothing,
    observ = "2021-01-13",#"2021-02-06",#"2021-03-25"
    chains = 1,
    preds = "CT",#"CF,CC,CR,CS",
    semipara=false,
)
#----------------------------------------------------------------------------
# load data
@info "load data"
data = National.load_data(
      observationsend  = ps.observ
    , predictors        = ps.preds|> CovidSurvey.parse_predictors
    , rwstep           = 1
    , epidemicstart    = 30
    , numimpute        = 6
    , deathmodel        = ps.model == "cases" ? National.DeathInit(obs_stop="2020-07-01") : National.DeathInit()
    , casemodel         = National.CaseInit(obs_start="2020-06-01")
    # , hospitmodel       = National.HospitInit()
    , seromodel         = National.SeroInit1(delay=0, std=.5)
    , link              = KLogistic(3.)
    , invlink           = KLogit(3.)
    , lockdown          = "2020-03-18"
    , covariates_kwargs = Dict(
        :semiparametric => ps.semipara,
        :fname => projectdir("data","raw_contact_rates.csv"),
        :shift => -1,
        :startdate => "2020-10-01",
        :enddate => nothing))

turing_data = data.turing_data;
## ============================================================================

@info "prepare model"

model1 = National.name2model[ps.model]
!isnothing(ps.seed) && Random.seed!(ps.seed);
m = model1(turing_data, false)
p = m();
## ============================================================================
@info "start sampling"

@time chain = let
    thinning = 10
    if ps.chains > 1
        sample(m, NUTS(ps.warmup, 0.95; max_depth=5), MCMCThreads(), ps.steps + ps.warmup, ps.chains; progress=true, thinning)
    else
        sample(m, NUTS(ps.warmup, 0.95; max_depth=5), ps.steps + ps.warmup; progress=true, thinning) #; max_depth=15
    end
end
## ==========================================================================
@info "save"

fdir = projectdir("reports/tmp", "$(today())")
fname = National.save_results(fdir, ps, data, chain)

## ==========================================================================
@info "post processing"
National.postprocessing(fname; plot_results=true)

## ==========================================================================
@info "post processing"
fdir, ps, ignores = National.parse_fname(fname)
p = National.PostProcessing3(fdir, ps, ignores, [], fname)
National.savechain(p)
## ==========================================================================
@info "plot chain"
plot_results = true
National.plot_chains(p; plot_results)

## ==========================================================================
@info "meanplot"
National.plot_means(p; plot_results)
## ==========================================================================
# perform diagnostics
p = National.skip_warmup(p)
National.diagnostics(p)
## ==========================================================================
@info "make predictions"
gp = National.generate_posterior(p);
## ==========================================================================
@info "plot regions"
National.plot_prediction(p, gp; plot_results=true)
## -----------------------------------------------------------------------------
@info "store reproduction number"
National.save_rt(p, gp)

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
# Common
# ------------------------------------------------------------
function add_regularization_params!(regularization::AbstractRegularizer; kwargs...)
    push!(regularization.parameters, kwargs...)
    return nothing
end

function decision(::AbstractLShapedSolver, regularization::AbstractRegularization)
    return regularization.ξ
end

function objective(::AbstractLShapedSolver, regularization::AbstractRegularization)
    return regularization.data.Q̃
end

function solve_regularized_master!(lshaped::AbstractLShapedSolver, ::AbstractRegularization)
    lshaped.mastersolver(lshaped.mastervector)
    return nothing
end

function gap(lshaped::AbstractLShapedSolver, regularization::AbstractRegularization)
    @unpack θ = lshaped.data
    @unpack Q̃ = regularization.data
    return abs(θ-Q̃)/(abs(Q̃)+1e-10)
end

function process_cut!(lshaped::AbstractLShapedSolver, cut::AbstractHyperPlane, ::AbstractRegularization)
    return nothing
end

function project!(lshaped::AbstractLShapedSolver, ::AbstractRegularization)
    return nothing
end

function add_regularization_params!(regularization::AbstractRegularization; kwargs...)
    push!(regularization.params, kwargs...)
    return nothing
end

ext: jl


text:
"""
    DDVFA.jl

Description:
    Includes all of the structures and logic for running a Distributed Dual-Vigilance Fuzzy ART (DDVFA) module.
"""

# --------------------------------------------------------------------------- #
# OPTIONS
# --------------------------------------------------------------------------- #

"""
    opts_DDVFA()

Distributed Dual Vigilance Fuzzy ART options struct.

# Examples
```julia-repl
julia> my_opts = opts_DDVFA()
```
"""
@with_kw mutable struct opts_DDVFA <: ARTOpts @deftype Float
    # Lower-bound vigilance parameter: [0, 1]
    rho_lb = 0.7; @assert rho_lb >= 0.0 && rho_lb <= 1.0
    # Upper bound vigilance parameter: [0, 1]
    rho_ub = 0.85; @assert rho_ub >= 0.0 && rho_ub <= 1.0
    # Choice parameter: alpha > 0
    alpha = 1e-3; @assert alpha > 0.0
    # Learning parameter: (0, 1]
    beta = 1.0; @assert beta > 0.0 && beta <= 1.0
    # "Pseudo" kernel width: gamma >= 1
    gamma = 3.0; @assert gamma >= 1.0
    # "Reference" gamma for normalization: 0 <= gamma_ref < gamma
    gamma_ref = 1.0; @assert 0.0 <= gamma_ref && gamma_ref < gamma
    # Similarity method (activation and match):
    #   'single', 'average', 'complete', 'median', 'weighted', or 'centroid'
    method::String = "single"
    # Display flag
    display::Bool = true
    # Maximum number of epochs during training
    max_epoch::Int = 1
    # Normalize the threshold by the feature dimension
    gamma_normalization::Bool = true
end # opts_DDVFA

# --------------------------------------------------------------------------- #
# STRUCTS
# --------------------------------------------------------------------------- #

"""
    DDVFA <: ART

Distributed Dual Vigilance Fuzzy ARTMAP module struct.

# Examples
```julia-repl
julia> DDVFA()
DDVFA
    opts: opts_DDVFA
    subopts::opts_FuzzyART
    ...
```
"""
mutable struct DDVFA <: ART
    # Get parameters
    opts::opts_DDVFA
    subopts::opts_FuzzyART
    config::DataConfig

    # Working variables
    threshold::Float
    F2::Vector{FuzzyART}
    labels::IntegerVector
    n_categories::Int
    epoch::Int
    T::Float
    M::Float
end # DDVFA <: ART

# --------------------------------------------------------------------------- #
# CONSTRUCTORS
# --------------------------------------------------------------------------- #

"""
    DDVFA()

Implements a DDVFA learner with default options.

# Examples
```julia-repl
julia> DDVFA()
DDVFA
    opts: opts_DDVFA
    subopts: opts_FuzzyART
    ...
```
"""
function DDVFA()
    opts = opts_DDVFA()
    DDVFA(opts)
end # DDVFA()

"""
    DDVFA(;kwargs...)

Implements a DDVFA learner with keyword arguments.

# Examples
```julia-repl
julia> DDVFA(rho_lb=0.4, rho_ub = 0.75)
DDVFA
    opts: opts_DDVFA
    subopts: opts_FuzzyART
    ...
```
"""
function DDVFA(;kwargs...)
    opts = opts_DDVFA(;kwargs...)
    DDVFA(opts)
end # DDVFA(;kwargs...)

"""
    DDVFA(opts::opts_DDVFA)

Implements a DDVFA learner with specified options.

# Examples
```julia-repl
julia> my_opts = opts_DDVFA()
julia> DDVFA(my_opts)
DDVFA
    opts: opts_DDVFA
    subopts: opts_FuzzyART
    ...
```
"""
function DDVFA(opts::opts_DDVFA)
    # Set the options used for all F2 FuzzyART modules
    subopts = opts_FuzzyART(
        rho=opts.rho_ub,
        gamma=opts.gamma,
        gamma_ref=opts.gamma_ref,
        gamma_normalization=opts.gamma_normalization,
        display=false
    )

    # Construct the DDVFA module
    DDVFA(opts,
          subopts,
          DataConfig(),
          0.0,
          Array{FuzzyART}(undef, 0),
          Array{Int}(undef, 0),
          0,
          0,
          0.0,
          0.0
    )
end # DDVFA(opts::opts_DDVFA)

# --------------------------------------------------------------------------- #
# ALGORITHMIC METHODS
# --------------------------------------------------------------------------- #

"""
    set_threshold!(art::DDVFA)

Sets the vigilance threshold of the DDVFA module as a function of several flags and hyperparameters.
"""
function set_threshold!(art::DDVFA)
    # Gamma match normalization
    if art.opts.gamma_normalization
        # Set the learning threshold as a function of the data dimension
        art.threshold = art.opts.rho_lb*(art.config.dim^art.opts.gamma_ref)
    else
        # Set the learning threshold as simply the vigilance parameter
        art.threshold = art.opts.rho_lb
    end
end # set_threshold!(art::DDVFA)

"""
    train!(art::DDVFA, x::RealMatrix ; y::IntegerVector=Vector{Int}(), preprocessed::Bool=false)

Train the DDVFA model on the data.
"""

"""
    train!(art::DDVFA, x::RealVector ; y::Integer=0, preprocessed::Bool=false)
"""
function train!(art::DDVFA, x::RealVector ; y::Integer=0, preprocessed::Bool=false)
    # Flag for if training in supervised mode
    supervised = !iszero(y)

    # Run the sequential initialization procedure
    sample = init_train!(x, art, preprocessed)

    # Initialization
    if isempty(art.F2)
        # Set the threshold
        set_threshold!(art)
        # Set the first label as either 1 or the first provided label
        y_hat = supervised ? y : 1
        # Create a new category
        create_category(art, sample, y_hat)
        return y_hat
    end

    # Default to mismatch
    mismatch_flag = true

    # Compute the activation for all categories
    T = zeros(art.n_categories)
    for jx = 1:art.n_categories
        activation_match!(art.F2[jx], sample)
        T[jx] = similarity(art.opts.method, art.F2[jx], "T", sample, art.opts.gamma_ref)
    end

    # Compute the match for each category in the order of greatest activation
    index = sortperm(T, rev=true)
    for jx = 1:art.n_categories
        bmu = index[jx]
        M = similarity(art.opts.method, art.F2[bmu], "M", sample, art.opts.gamma_ref)
        # If we got a match, then learn (update the category)
        if M >= art.threshold
            # Update the stored match and activation values
            art.M = M
            art.T = T[bmu]
            # If supervised and the label differs, trigger mismatch
            if supervised && (art.labels[bmu] != y)
                break
            end
            # Update the weights with the sample
            train!(art.F2[bmu], sample, preprocessed=true)
            # Save the output label for the sample
            y_hat = art.labels[bmu]
            # No mismatch
            mismatch_flag = false
            break
        end
    end

    # If we triggered a mismatch
    if mismatch_flag
        # Update the stored match and activation values
        bmu = index[1]
        art.M = similarity(art.opts.method, art.F2[bmu], "M", sample, art.opts.gamma_ref)
        art.T = T[bmu]
        # Get the correct label
        y_hat = supervised ? y : art.n_categories + 1
        create_category(art, sample, y_hat)
    end

    return y_hat
end # train!(art::DDVFA, x::RealVector ; y::Integer=0, preprocessed::Bool=false)

"""
    create_category(art::DDVFA, sample::RealVector, label::Integer)

Create a new category by appending and initializing a new FuzzyART node to F2.
"""
function create_category(art::DDVFA, sample::RealVector, label::Integer)
    # Global Fuzzy ART
    art.n_categories += 1
    push!(art.labels, label)
    # Local Gamma-Normalized Fuzzy ART
    push!(art.F2, FuzzyART(art.subopts, sample, preprocessed=true))
end # function create_category(art::DDVFA, sample::RealVector, label::Integer)

"""
    stopping_conditions(art::DDVFA)

Stopping conditions for Distributed Dual Vigilance Fuzzy ARTMAP.

Returns true if there is no change in weights during the epoch or the maxmimum epochs has been reached.
"""
function stopping_conditions(art::DDVFA)
    # Compute the stopping condition, return a bool
    return art.epoch >= art.opts.max_epoch
end # stopping_conditions(DDVFA)

"""
    similarity(method::String, F2::FuzzyART, field_name::String, sample::RealVector, gamma_ref::RealFP)

Compute the similarity metric depending on method with explicit comparisons
for the field name.
"""
function similarity(method::String, F2::FuzzyART, field_name::String, sample::RealVector, gamma_ref::RealFP)
    @debug "Computing similarity"

    if field_name != "T" && field_name != "M"
        error("Incorrect field name for similarity metric.")
    end
    # Single linkage
    if method == "single"
        if field_name == "T"
            value = maximum(F2.T)
        elseif field_name == "M"
            value = maximum(F2.M)
        end
    # Average linkage
    elseif method == "average"
        if field_name == "T"
            value = mean(F2.T)
        elseif field_name == "M"
            value = mean(F2.M)
        end
    # Complete linkage
    elseif method == "complete"
        if field_name == "T"
            value = minimum(F2.T)
        elseif field_name == "M"
            value = minimum(F2.M)
        end
    # Median linkage
    elseif method == "median"
        if field_name == "T"
            value = median(F2.T)
        elseif field_name == "M"
            value = median(F2.M)
        end
    # Weighted linkage
    elseif method == "weighted"
        if field_name == "T"
            value = F2.T' * (F2.n_instance ./ sum(F2.n_instance))
        elseif field_name == "M"
            value = F2.M' * (F2.n_instance ./ sum(F2.n_instance))
        end
    # Centroid linkage
    elseif method == "centroid"
        # Get the minimum of each weight element, cast to a 1-D vector
        Wc = vec(minimum(F2.W, dims=2))
        T = norm(element_min(sample, Wc), 1) / (F2.opts.alpha + norm(Wc, 1))^F2.opts.gamma
        if field_name == "T"
            value = T
        elseif field_name == "M"
            value = (norm(Wc, 1)^gamma_ref)*T
        end
    else
        error("Invalid/unimplemented similarity method")
    end

    return value
end # similarity(method::String, F2::FuzzyART, field_name::String, sample::RealVector, gamma_ref::RealFP)

"""
    classify(art::DDVFA, x::RealMatrix ; preprocessed::Bool=false, get_bmu::Bool=false)

Predict categories of 'x' using the DDVFA model.

Returns predicted categories 'y_hat.'

# Examples
```julia-repl
julia> my_DDVFA = DDVFA()
DDVFA
    opts: opts_DDVFA
    ...
julia> x, y = load_data()
julia> train!(my_DDVFA, x)
julia> y_hat = classify(my_DDVFA, y)
```
"""
function classify(art::DDVFA, x::RealVector ; preprocessed::Bool=false, get_bmu::Bool=false)
    # Preprocess the data
    sample = init_classify!(x, art, preprocessed)

    # Calculate all global activations
    T = zeros(art.n_categories)
    for jx = 1:art.n_categories
        activation_match!(art.F2[jx], sample)
        T[jx] = similarity(art.opts.method, art.F2[jx], "T", sample, art.opts.gamma_ref)
    end

    # Sort by highest activation
    index = sortperm(T, rev=true)

    # Default to mismatch
    mismatch_flag = true

    # Iterate over the list of activations
    for jx = 1:art.n_categories
        # Get the best-matching unit
        bmu = index[jx]
        # Get the match value of this activation
        M = similarity(art.opts.method, art.F2[bmu], "M", sample, art.opts.gamma_ref)
        # If the match satisfies the threshold criterion, then report that label
        if M >= art.threshold
            # Update the stored match and activation values
            art.M = M
            art.T = T[bmu]
            # Current winner
            y_hat = art.labels[bmu]
            mismatch_flag = false
            break
        end
    end

    # If we did not find a resonant category
    if mismatch_flag
        @debug "Mismatch"
        # Update the stored match and activation values of the best matching unit
        bmu = index[1]
        art.M = similarity(art.opts.method, art.F2[bmu], "M", sample, art.opts.gamma_ref)
        art.T = T[bmu]
        # Report either the best matching unit or the mismatch label -1
        y_hat = get_bmu ? art.labels[bmu] : -1
    end

    return y_hat
end

# --------------------------------------------------------------------------- #
# CONVENIENCE METHODS
# --------------------------------------------------------------------------- #

"""
    get_W(art::DDVFA)

Convenience functio; return a concatenated array of all DDVFA weights.
"""
function get_W(art::DDVFA)
    # Return a concatenated array of the weights
    return [art.F2[kx].W for kx = 1:art.n_categories]
end # get_W(art::DDVFA)

"""
    get_n_weights_vec(art::DDVFA)

Convenience function; return the number of weights in each category as a vector.
"""
function get_n_weights_vec(art::DDVFA)
    return [art.F2[i].n_categories for i = 1:art.n_categories]
end # get_n_weights_vec(art::DDVFA)

"""
    get_n_weights(art::DDVFA)

Convenience function; return the sum total number of weights in the DDVFA module.
"""
function get_n_weights(art::DDVFA)
    # Return the number of weights across all categories
    return sum(get_n_weights_vec(art))
end # get_n_weights(art::DDVFA)


text:
module RobustModels

include("compat.jl")

# Use README as the docstring of the module and doctest README
@doc let path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    replace(read(path, String), "```julia" => "```jldoctest; output = false")
end RobustModels

using Distributions: ccdf, pdf, quantile, Normal, Chisq, TDist, FDist
using SparseArrays: SparseMatrixCSC, spdiagm
using LinearAlgebra: diag, dot, tr, I, UniformScaling, rmul!, lmul!
using Random: AbstractRNG, GLOBAL_RNG
using Printf: @printf, @sprintf
using GLM: Link, canonicallink, FPVector, lm, SparsePredChol, DensePredChol
using StatsBase: mean, mad, ConvergenceException, sample, quantile
using IterativeSolvers: lsqr!, cg!
#using Roots: find_zero, Order1, ConvergenceFailed
#using QuadGK: quadgk
#import Tulip

import Base: ==, show, broadcastable
import GLM:
    dispersion, LinPred, DensePred, ModResp, delbeta!, linpred!, installbeta!, cholpred
import StatsBase: fit, fit!, deviance, nulldeviance, nobs, weights, Weights, confint,
                  dof, dof_residual, loglikelihood, nullloglikelihood, stderror,
                  vcov, residuals, predict, response, islinear, fitted, isfitted,
                  mean, var, std, sem, mean_and_std, mean_and_var
import StatsModels:
    @delegate, @formula, RegressionModel, coef, coeftable, CoefTable,
    leverage, modelmatrix, TableRegressionModel

## Reexports
export coef,
       coeftable,
       confint,
       deviance,
       nulldeviance,
       dof,
       dof_residual,
       loglikelihood,
       nullloglikelihood,
       nobs,
       stderror,
       vcov,
       residuals,
       predict,
       fit,
       fit!,
       response,
       modelmatrix,
       dispersion,
       islinear,
       isfitted,
       fitted,
       weights,
       leverage,
       quantile,
       @formula,
       mean,
       var,
       std,
       sem,
       mean_and_std,
       mean_and_var,
       nothing  # stopper

export LossFunction,
       BoundedLossFunction,
       ConvexLossFunction,
       CompositeLossFunction,
       L2Loss,
       L1Loss,
       HuberLoss,
       L1L2Loss,
       FairLoss,
       LogcoshLoss,
       ArctanLoss,
       CauchyLoss,
       GemanLoss,
       WelschLoss,
       TukeyLoss,
       YohaiZamarLoss,
       AbstractEstimator,
       AbstractQuantileEstimator,
       MEstimator,
       SEstimator,
       MMEstimator,
       TauEstimator,
       GeneralizedQuantileEstimator,
       ExpectileEstimator,
       L2Estimator,
       DensePredCG,
       SparsePredCG,
       RidgePred,
       RobustResp,
       AbstractRobustModel,
       RobustLinearModel,
       QuantileRegression,
       Estimator,
       rlm,
       quantreg,
       loss,
       tuning_constant,
       refit!,
       projectionmatrix,
       workingweights,
       scale,
       tauscale,
       mean_and_sem,
       nothing  # stopper

"""
An estimator needs a cost/loss function for the modified (weighted) least squares
problems of the form:

```math
\\min \\sum_i \\rho\\left(\\dfrac{r_i}{\\hat{\\sigma}}\right)
```
"""
abstract type LossFunction end

"Bounded loss function type for hard rejection of outlier."
abstract type BoundedLossFunction <: LossFunction end

"Convex loss function type with no local minimum solution."
abstract type ConvexLossFunction <: LossFunction end

"""
A robust estimator is a location or scale estimator associated to one (or more) loss function and
solved using Iteratively Reweighted Least Square.
"""
abstract type AbstractEstimator end

"Generalized M-Quantile estimator"
abstract type AbstractQuantileEstimator <: AbstractEstimator end

"""
    AbstractRobustModel

Abstract type for robust models.

`RobustModels.jl` implements one subtype: [`RobustLinearModel`](@ref).
See the documentation for each for more details.
"""
abstract type AbstractRobustModel{T} <: RegressionModel end

abstract type RobustResp{T} <: ModResp end

abstract type AbstractRegularizedPred{T} <: LinPred end

Base.broadcastable(m::T) where {T<:AbstractEstimator} = Ref(m)
Base.broadcastable(m::T) where {T<:LossFunction} = Ref(m)

include("estimators.jl")
include("robustlinearmodel.jl")
include("linpred.jl")
include("linresp.jl")
include("univariate.jl")
include("quantileregression.jl")

end # module


text:
using Quoridor
using Base.Test

Q = Quoridor

# test movement blockage

@test !Q.separates(Q.Wall((6,1), (6,3)), Q.Loc(5,1), Q.Loc(5,2))
@test Q.separates(Q.Wall((4,3), (6,3)), Q.Loc(4,3), Q.Loc(4,2))


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

using Markdown
using InteractiveUtils

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

# ╔═╡ 1e992cf1-7ea6-4c23-a573-3c867c22ada0
md"""
# Scrubbable numbers

Try clicking and dragging the numbers in the text below:
"""

# ╔═╡ 023869e9-236d-426e-8419-d398f1d20c3e
html"""

"""

# ╔═╡ 7777f629-dc29-44a8-9127-8e236c88d1ef
# using PlutoUI

# ╔═╡ d86baba8-1207-466b-9170-fdc5f2616152
html"""

"""

# ╔═╡ 10494af3-3d7d-47ae-8844-3ca139e7708b
md"""
## Examples

We define a variable `x` as a _scrubbable number_ like so:
"""

# ╔═╡ 99a2907f-f748-47bd-92d2-88eb21e8e0f1
md"""
Using interpolation (with `$`), you can write this definition inside Markdown text:
"""

# ╔═╡ 0be135f4-ef36-41d1-8a20-1080ca8919a0
md"""
## Arguments

Besides an initial value, a scrubbable number also has an array of possible values that can be reached.

When you pass a **single number** to `Scrubbable`, this array is automatically created. Have a look at the minimum and maximum allowed values of the following scrubbables:
"""

# ╔═╡ 207d3d2f-2841-4eca-839d-88f80c80c788
md"""
_(Here we created a `Scrubbable` on its own, without binding its value to a variable. Not very useful!)_

You can also **specify the array manually**:
"""

# ╔═╡ a1d5ebd5-dec1-4486-9968-a9d6a539751a
md"""
If no `default` is specified, the middle value is used.
"""

# ╔═╡ 270beb18-8288-48eb-9165-423e12ec8411
md"""
### Formatting

The library [`d3-format`](https://github.com/d3/d3-format) is used to format floating-point numbers. You can specify a **format string** like `".2f"` to be used to format the scrubbable value. Have a look at their [documentation](https://github.com/d3/d3-format) to see more examples.
"""

# ╔═╡ c17663e3-8a3b-4d6b-8296-f8b996739a01
html"""

"""

# ╔═╡ f8b9398e-afe6-4fdd-bc7a-9bb7b7c3cdc2
md"""
# Appendix
## What is a good default range?
"""

# ╔═╡ 4e40dcad-9840-4517-9947-7fcd4dce69cd
function default_range(x::Integer)
	if x == 0
		-10:10
	elseif 0 < x <= 10
		0:10
	elseif -10 <= x < 0
		-10:0
	else
		sort(round.([Int64], (0.0:0.1:2.0) .* x))
	end
end

# ╔═╡ e14b3395-831d-43d6-874e-cfb32f8edd05
up_or_down_one_order_of_magnitude = 10 .^ (-1.0:0.1:1.0)

# ╔═╡ 018035dd-78d7-46d3-bd3c-5c136fa47929
function default_range(x::Real) # not an integer
	if x == 0
		-1.0 : 0.1 : 1.0
	elseif 0 < x < 1
		between_zero_and_x = LinRange(0.0, x, 10)
		between_x_and_one = LinRange(x, 1.0, 10)
		[between_zero_and_x..., between_x_and_one[2:end]...]
	else
		sort(x .* up_or_down_one_order_of_magnitude)
	end
end

# ╔═╡ e0ea8328-ddcc-4130-b2f6-0079f98402ad
md"""
Integer inputs become integer ranges with constant step:
"""

# ╔═╡ e8a6ce9f-9f48-4e31-ad2a-e21d7e0415af
default_range(20)

# ╔═╡ 1e139b79-b703-43bf-a351-1a2d811aea6a
default_range(2000)

# ╔═╡ 1e60fea6-33dc-417b-9928-eb5922b41759
md"""
Floating points get a logarithmic scale across one order of magnitude:
"""

# ╔═╡ 6fe1f0f7-e8a1-441a-a676-ddd4b3b26668
default_range(1.0e6)

# ╔═╡ cfcc5fbb-72f7-4223-a83e-96fe6971b143
default_range(-2.0)

# ╔═╡ 1fa8edc6-9aa2-4082-bdeb-7517d9e2dd71
md"""
Zero becomes a range around zero:
"""

# ╔═╡ 3f1c3fa5-2257-4c3a-aa75-0b3c59a7fcdc
default_range(0)

# ╔═╡ 2554121f-13e6-4b07-9c45-b2ccf154d07d
default_range(0.0)

# ╔═╡ 8bce6b13-9600-49ca-a7ea-ba2f53028f1b
# using HypertextLiteral

# ╔═╡ d17d259b-8379-46a7-ab54-cd2f697ec713
# @htl("""
# 	$(bc)
# 	$(cool)
# 	""")

# ╔═╡ aed5fa58-4fe3-4596-b18d-a76cd98a5a1b
md"""
## Definition
"""

# ╔═╡ 1a649975-9e31-4ae8-8b2c-b615852cfc9d
function skip_as_script(m::Module)
	if isdefined(m, :PlutoForceDisplay)
		return m.PlutoForceDisplay
	else
		isdefined(m, :PlutoRunner) && parentmodule(m) == Main
	end
end

# ╔═╡ 8c1d2b3b-8fa1-4356-8a9d-dff10cd0a336
if skip_as_script(@__MODULE__)
	import Pkg
	Pkg.activate(Base.current_project(@__DIR__))
	Text("Project env active")
end

# ╔═╡ 58a5b0d2-88a6-4a83-bb26-05c05a51716b
import AbstractPlutoDingetjes

# ╔═╡ 9c7ce2da-4ad8-11eb-14cd-cfcc8d2a6bf8
begin
	"""

	An inline number that can be changed by clicking and dragging the mouse.

	# Examples
	```julia
	md\"\"\"
	_If Alice has \$(@bind a Scrubbable(20)) apples,
	and she gives \$(@bind b Scrubbable(3)) apples to Bob..._
	\"\"\"
	```
	```julia
	md\"\"\"
	_...then Alice has **\$(a - b)** apples left._
	\"\"\"
	```

	In the examples above, we give the **initial value** as parameter, and the reader can change it to be lower or higher.

	## Custom range
	Besides an initial value, a scrubbable number also has an array of possible values that can be reached. When you pass a **single number** to `Scrubbable`, this array is automatically created.

	You can also **specify the array manually**:

	```julia
	@bind apples Scrubbable(200 : 300; default=220)
	```

	## Formatting

	The library [`d3-format`](https://github.com/d3/d3-format) is used to format floating-point numbers. You can specify a **format string** like `".2f"` to be used to format the scrubbable value. Have a look at their [documentation](https://github.com/d3/d3-format) to see more examples.

	`@bind money Scrubbable(30e6, format=".0s", prefix="€ ")`

	`@bind coolness Scrubbable(0.80 : 0.01 : 1.00, format=".0%", prefix="you are 🌝 ", suffix=" cool")`
	"""
	Base.@kwdef struct Scrubbable
		values::AbstractVector{<:Real}
		default::Real
		format::Union{AbstractString,Nothing}=nothing
		prefix::AbstractString=""
		suffix::AbstractString=""
        id::String=join(rand('a':'z', 16))
	end
	Scrubbable(range::AbstractVector; kwargs...) = Scrubbable(;values=range, default=range[1 + length(range) ÷ 2], kwargs...)
	Scrubbable(x::Number; kwargs...) = Scrubbable(;values=default_range(x), default=x, kwargs...)

	Base.get(s::Scrubbable) = s.default

	AbstractPlutoDingetjes.Bonds.initial_value(s::Scrubbable) =
		s.default
	AbstractPlutoDingetjes.Bonds.possible_values(s::Scrubbable) =
		s.values
	function AbstractPlutoDingetjes.Bonds.validate_value(s::Scrubbable, val)
		val isa Real && minimum(s.values) - 0.001 <= val <= maximum(s.values) + 0.001
	end

	function Base.show(io::IO, m::MIME"text/html", s::Scrubbable)
		format = if s.format === nothing
			# TODO: auto format
			if eltype(s.values) <: Integer
				""
			else
				".1f"
			end
		else
			String(s.format)
		end

		write(io, """""")
	end

	Scrubbable
end

# ╔═╡ b62db8c0-4352-4d0f-83a2-ac170ef3337a
md"""
_If Alice has $(@bind a Scrubbable(20)) apples,
and she gives $(@bind b Scrubbable(3)) apples to Bob..._
"""

# ╔═╡ 861a570a-3a77-4445-9217-3d38682cbb8c
md"""
_...then Alice has **$(a - b)** apples left._
"""

# ╔═╡ 5d6299e0-fb4a-448f-b87a-9d17b607fd6a
@bind x Scrubbable(5)

# ╔═╡ 188e548a-2d5d-4680-9b0a-7bfb03b24dab
1000 + x

# ╔═╡ e9169074-6fec-43a6-badc-4b20d05050ab
md"""
If Alice has $(@bind num_apples Scrubbable(20)) apples...
"""

# ╔═╡ 65584961-4085-4fd6-8968-6951056f3b1c
num_apples

# ╔═╡ 9e97e893-066c-4fc6-b380-9637994b7dac
(Scrubbable(5), Scrubbable(2000), Scrubbable(30.0))

# ╔═╡ b69fb562-7e8e-4366-a53d-85173afc49be
Scrubbable(200 : 300; default=220)

# ╔═╡ 19549e7c-cf8a-4d79-814b-db89efc91a15
md"""
# I spent $(@bind money Scrubbable(30e6, format=".0s", prefix="€ ")) on Pluto.jl stickers

And it was worth it!

"""

# ╔═╡ 71bb41eb-908c-4b43-881b-5b6726fc6d0a
@bind coolness Scrubbable(0.80 : 0.01 : 1.00, format=".0%", prefix="you are 🌝 ", suffix=" cool")

# ╔═╡ 574564ab-ab15-490e-aff6-6300718ae751
if coolness >= 1
	md"![](https://media.giphy.com/media/GwGXoeb0gm7sc/giphy.gif)"
end

# ╔═╡ b081aa76-f080-4dd0-bcff-4bcc82a1c50a
bc = @bind cool Scrubbable(199.1)

# ╔═╡ 5c78e637-5733-4b0a-8dab-bf1cd9656d11
HTML(join(repr.([MIME"text/html"()], [Scrubbable(1.0) for _ in 1:100])))

# ╔═╡ 1d34fec8-01cb-4bee-8144-d8cc13a87b8b
export Scrubbable

# ╔═╡ Cell order:
# ╟─1e992cf1-7ea6-4c23-a573-3c867c22ada0
# ╟─023869e9-236d-426e-8419-d398f1d20c3e
# ╠═7777f629-dc29-44a8-9127-8e236c88d1ef
# ╠═b62db8c0-4352-4d0f-83a2-ac170ef3337a
# ╠═861a570a-3a77-4445-9217-3d38682cbb8c
# ╟─d86baba8-1207-466b-9170-fdc5f2616152
# ╟─10494af3-3d7d-47ae-8844-3ca139e7708b
# ╠═5d6299e0-fb4a-448f-b87a-9d17b607fd6a
# ╠═188e548a-2d5d-4680-9b0a-7bfb03b24dab
# ╟─99a2907f-f748-47bd-92d2-88eb21e8e0f1
# ╠═e9169074-6fec-43a6-badc-4b20d05050ab
# ╠═65584961-4085-4fd6-8968-6951056f3b1c
# ╟─0be135f4-ef36-41d1-8a20-1080ca8919a0
# ╠═9e97e893-066c-4fc6-b380-9637994b7dac
# ╟─207d3d2f-2841-4eca-839d-88f80c80c788
# ╠═b69fb562-7e8e-4366-a53d-85173afc49be
# ╟─a1d5ebd5-dec1-4486-9968-a9d6a539751a
# ╟─270beb18-8288-48eb-9165-423e12ec8411
# ╟─19549e7c-cf8a-4d79-814b-db89efc91a15
# ╟─71bb41eb-908c-4b43-881b-5b6726fc6d0a
# ╟─574564ab-ab15-490e-aff6-6300718ae751
# ╟─c17663e3-8a3b-4d6b-8296-f8b996739a01
# ╟─f8b9398e-afe6-4fdd-bc7a-9bb7b7c3cdc2
# ╠═4e40dcad-9840-4517-9947-7fcd4dce69cd
# ╠═e14b3395-831d-43d6-874e-cfb32f8edd05
# ╠═018035dd-78d7-46d3-bd3c-5c136fa47929
# ╟─e0ea8328-ddcc-4130-b2f6-0079f98402ad
# ╠═e8a6ce9f-9f48-4e31-ad2a-e21d7e0415af
# ╠═1e139b79-b703-43bf-a351-1a2d811aea6a
# ╟─1e60fea6-33dc-417b-9928-eb5922b41759
# ╠═6fe1f0f7-e8a1-441a-a676-ddd4b3b26668
# ╠═cfcc5fbb-72f7-4223-a83e-96fe6971b143
# ╟─1fa8edc6-9aa2-4082-bdeb-7517d9e2dd71
# ╠═3f1c3fa5-2257-4c3a-aa75-0b3c59a7fcdc
# ╠═2554121f-13e6-4b07-9c45-b2ccf154d07d
# ╠═b081aa76-f080-4dd0-bcff-4bcc82a1c50a
# ╠═8bce6b13-9600-49ca-a7ea-ba2f53028f1b
# ╠═d17d259b-8379-46a7-ab54-cd2f697ec713
# ╠═5c78e637-5733-4b0a-8dab-bf1cd9656d11
# ╟─aed5fa58-4fe3-4596-b18d-a76cd98a5a1b
# ╟─1a649975-9e31-4ae8-8b2c-b615852cfc9d
# ╟─8c1d2b3b-8fa1-4356-8a9d-dff10cd0a336
# ╠═58a5b0d2-88a6-4a83-bb26-05c05a51716b
# ╠═9c7ce2da-4ad8-11eb-14cd-cfcc8d2a6bf8
# ╠═1d34fec8-01cb-4bee-8144-d8cc13a87b8b


text:
using SeqUMAP
using Test

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

max_stars_repo_name: alecpnkw/SeqUMAP

max_issues_repo_name: alecpnkw/SeqUMAP

max_forks_repo_name: alecpnkw/SeqUMAP


text:
module SegmentationTools

using AxisArrays
using ColorSchemes
using Colors
using DataFrames
using DataStructures
using FreeTypeAbstraction
using ImageAxes
using ImageSegmentation
using ImageMorphology
using ImageMetadata
using Images
using KernelDensity
using MappedArrays
using Statistics
using SimpleTraits

include("utils.jl")
include("segment.jl")
include("diagnostics.jl")

export get_localities

end # module


text:
### ##########################################################
### ####
### #### import time series from yahoo or FRED and return DataFrame
### ####
### ##########################################################
###
### ################# yahoo
###
### function yahoo(stock::String, fm::Int, fd::Int, fy::Int, tm::Int, td::Int, ty::Int, period::String)
###
### # take care of yahoo's 0 indexing for month
###   fm-=1
###   tm-=1
###
###   ydata = download("http://ichart.finance.yahoo.com/table.csv?s=$stock&a=$fm&b=$fd&c=$fy&d=$tm&e=$td&f=$ty&g=$period", "y.csv")
###   df    = readtime("y.csv")
###   run(`rm y.csv`)
###
###   return df
### end
###
### ################# fred
###
### function fred(econdata::String, colname::String)
###
###   fdata = download("http://research.stlouisfed.org/fred2/series/$econdata/downloaddata/$econdata.csv", "f.csv")
###   df    = readtime("f.csv")
###   run(`rm f.csv`)
###
###   return df
### end
###
### ##########################################################
### ####
### #### import time series and return Stock object
### ####
### ##########################################################
###
### function yahoo_stock(stock::String, fm::Int, fd::Int, fy::Int, tm::Int, td::Int, ty::Int, period::String)
###   typed_stock = yahoo(stock, fm, fd, fy, tm, td, ty, period)
###
###   res = Stock(stock,
###               IndexedVector(typed_stock["Date"]),
###               typed_stock[["Open", "High", "Low", "Close", "Adj"]],
###               typed_stock["Vol"],
###               .01)
### end
###
### ##########################################################
### ####
### #### import Yahoo option series
### ####
### ##########################################################
###
### function optionseries(ticker::String)
###   # code here to return a DataFrame
### end
###
### function optionseries!(ticker::String)
###   # code here to return an OptionSeries object
### end
###
### ############### DEFAULT ###########################
###
### ## last three years, daily data
### yahoo(stock::String) = yahoo(stock::String,
###                              month(now()),
###                              day(now()),
###                              year(now())-3,
###                              month(now()),
###                              day(now()),
###                              year(now()), "d")
###
### ## last three years, daily data
### yahoo_stock(stock::String) = yahoo(stock::String,
###                              month(now()),
###                              day(now()),
###                              year(now())-3,
###                              month(now()),
###                              day(now()),
###                              year(now()), "d")
###
### ## name the fred column "VALUE"
### fred(econdata::String) = fred(econdata::String, "VALUE")


text:
"""
    PublicKey

Stores a public key generated by a `KeyGenerator` instance.

See also: [`KeyGenerator`](@ref)
"""
mutable struct PublicKey <: SEALObject
  handle::Ptr{Cvoid}

  function PublicKey()
    handleref = Ref{Ptr{Cvoid}}(C_NULL)
    retval = ccall((:PublicKey_Create1, libsealc), Clong,
                   (Ref{Ptr{Cvoid}},),
                   handleref)
    @check_return_value retval
    return PublicKey(handleref[])
  end

  function PublicKey(handle::Ptr{Cvoid})
    object = new(handle)
    finalizer(destroy!, object)
    return object
  end
end

function destroy!(object::PublicKey)
  if isallocated(object)
    @check_return_value ccall((:PublicKey_Destroy, libsealc), Clong, (Ptr{Cvoid},), object)
    sethandle!(object, C_NULL)
  end

  return nothing
end

function parms_id(key::PublicKey)
  parms_id = zeros(UInt64, 4)
  retval = ccall((:PublicKey_ParmsId, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{UInt64}),
                 key, parms_id)
  @check_return_value retval
  return parms_id
end


text:
using Distributed
addprocs(5)

using Test
using HTTP
using JSON

@testset "HTTP" begin
    for f in ["ascii.jl",
              "issue_288.jl",
              "utils.jl",
              "client.jl",
              "sniff.jl",
              "uri.jl",
              "url.jl",
              "cookies.jl",
              "parser.jl",
              "loopback.jl",
              "WebSockets.jl",
              "messages.jl",
              "handlers.jl",
              "server.jl",
              "async.jl",
              "aws4.jl",
              "insert_layers.jl"]

        println("Running $f tests...")
        include(f)
    end
end


text:
using Optim, HTTP, GLM, LinearAlgebra, Random, Statistics, DataFrames, CSV, FreqTables

#Worked with Ahmed and Aleeze

function ps3()

    #1

    url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2021/master/ProblemSets/PS3-gev/nlsw88w.csv"
    df = CSV.read(HTTP.get(url).body, DataFrame)
    X = [df.age df.white df.collgrad]
    Z = hcat(df.elnwage1, df.elnwage2, df.elnwage3, df.elnwage4, df.elnwage5, df.elnwage6, df.elnwage7, df.elnwage8)
    y = df.occupation

    function mlogit(Beta, X, Z, y)
        K = size(X,2)+1
        J = length(unique(y))
        N = length(y)
        bigY = zeros(N,J)
        for j=1:J
            bigY[:,j] = y.==j
        end
        bigBeta = [reshape(Beta,K,J-1) zeros(K)]

        bigZ = zeros(N,J)
        for j=1:J
            bigZ[:,j] = Z[:,j]-Z[:,J]
        end

        num = zeros(N,J)
        dem = zeros(N)
        for j=1:J
            XZ=cat(X,bigZ[:,j],dims=2)
            num[:,j] = exp.(XZ*bigBeta[:,j])
            dem .+= num[:,j]
        end

        P = num./repeat(dem,1,J)

        loglike = -sum( bigY.*log.(P) )

        return loglike
    end

    mlogit_hat_optim = optimize(b-> mlogit(b,X,Z, y), rand(7*(size(X,2)+1)), LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=50))
    println(mlogit_hat_optim.minimizer)

    #2. The coefficient γ can be interpreted as the change in utility with a percent increase in wage. We are able to make this interpretation from the fact that we are dealing with a linear-log model.

    #3.

    function nested_logit(alpha, X, Z, y, nesting_structure)

        beta = alpha[1:end-3]
        lambda = alpha[end-2:end-1]
        gamma = alpha[end]
        K = size(X,2)
        J = length(unique(y))
        N = length(y)
        bigY = zeros(N,J)
        for j=1:J
            bigY[:,j] = y.==j
        end
        bigBeta = [repeat(beta[1:K],1,length(nesting_structure[1])) repeat(beta[K+1:2K],1,length(nesting_structure[2])) zeros(K)]

        T = promote_type(eltype(X),eltype(alpha))
        num   = zeros(T,N,J)
        lidx  = zeros(T,N,J)
        dem   = zeros(T,N)

        for j=1:J
            if j in nesting_structure[1]
                lidx[:,j] = exp.( (X*bigBeta[:,j] .+ (Z[:,j] .- Z[:,J])*gamma)./lambda[1] )
            elseif j in nesting_structure[2]
                lidx[:,j] = exp.( (X*bigBeta[:,j] .+ (Z[:,j] .- Z[:,J])*gamma)./lambda[2] )
            else
                lidx[:,j] = exp.(zeros(N))
            end
        end

        for j=1:J
            if j in nesting_structure[1]
                num[:,j] = lidx[:,j].*sum(lidx[:,nesting_structure[1][:]];dims=2).^(lambda[1]-1)
            elseif j in nesting_structure[2]
                num[:,j] = lidx[:,j].*sum(lidx[:,nesting_structure[2][:]];dims=2).^(lambda[2]-1)
            else
                num[:,j] = lidx[:,j]
            end

            dem .+= num[:,j]
        end

        P = num./repeat(dem,1,J)

        loglike = -sum( bigY.*log.(P) )

        return loglike
    end

    nesting_structure = [[1 2 3], [4 5 6 7]]
    startvals = [2*rand(2*size(X,2)).-1; 1; 1; .1]

    td2 = TwiceDifferentiable(alpha -> nested_logit(alpha, X, Z, y, nesting_structure), startvals; autodiff = :forward)

    α_hat_nlogit = optimize(td2, startvals, LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=50))
    α_hat_nlogit_ad = α_hat_nlogit.minimizer

    H  = Optim.hessian!(td2, α_hat_nlogit_ad)
    α_hat_nlogit_ad_se = sqrt.(diag(inv(H)))

    println([α_hat_nlogit_ad α_hat_nlogit_ad_se])

    return nothing

end

ps3()


text:
export UserMean

"""
    UserMean(data::DataAccessor)

Recommend based on global user mean rating.
"""
struct UserMean <: Recommender
    data::DataAccessor
    scores::AbstractVector

    function UserMean(data::DataAccessor)
        n_user = size(data.R, 1)
        new(data, vector(n_user))
    end
end

isbuilt(recommender::UserMean) = isfilled(recommender.scores)

function build!(recommender::UserMean)
    n_user = size(recommender.data.R, 1)

    for u in 1:n_user
        v = recommender.data.R[u, :]
        recommender.scores[u] = mean(v)
    end
end

function predict(recommender::UserMean, u::Int, i::Int)
    check_build_status(recommender)
    recommender.scores[u]
end


text:
using Test
using UnPack
using MOONS.Grids
using MOONS.Fields

#####
##### Self-extrapolation
#####

extrap_simple!(f::AbstractField, grid::Grid, dim::Int) =
    extrap_simple!(f, grid, Val(dim))

function extrap_simple!(f::AbstractField, grid::Grid, ::Val{1})
    f[ 1 ,:,:] .= 2*f[  2  ,:,:] .- f[  3  ,:,:]
    f[end,:,:] .= 2*f[end-1,:,:] .- f[end-2,:,:]
end
function extrap_simple!(f::AbstractField, grid::Grid, ::Val{2})
    f[:, 1 ,:] .= 2*f[:,  2  ,:] .- f[:,  3  ,:]
    f[:,end,:] .= 2*f[:,end-1,:] .- f[:,end-2,:]
end
function extrap_simple!(f::AbstractField, grid::Grid, ::Val{3})
    f[:,:, 1 ] .= 2*f[:,:,  2  ] .- f[:,:,  3  ]
    f[:,:,end] .= 2*f[:,:,end-1] .- f[:,:,end-2]
end

@testset "Extrapolations - self extrap" begin
    FT = Float64
    n, a, b, β = 10, FT(0.0), FT(1.0), FT(1.01)
    grid = Cube(a, b, n; warpfun=Roberts_both, args=(β,))

    cent = CellCenter(grid)
    corn = CellCorner(grid)
    face = CellFaces(grid)
    edge = CellEdges(grid)
    cent_simple = CellCenter(grid)
    corn_simple = CellCorner(grid)
    face_simple = CellFaces(grid)
    edge_simple = CellEdges(grid)
    @inbounds for local_grid in GridField(grid, cent)
        @unpack ijk,x,y,z = local_grid
        cent[ijk] = cos(3*π*x)* cos(3*π*y)* cos(3*π*z)
        cent_simple[ijk] = cos(3*π*x)* cos(3*π*y)* cos(3*π*z)
    end
    @inbounds for local_grid in GridField(grid, corn)
        @unpack ijk,x,y,z = local_grid
        corn[ijk] = cos(3*π*x)* cos(3*π*y)* cos(3*π*z)
        corn_simple[ijk] = cos(3*π*x)* cos(3*π*y)* cos(3*π*z)
    end
    for dir in 1:3
        @inbounds for local_grid in GridField(grid, face[dir])
            @unpack ijk,x,y,z = local_grid
            face[dir][ijk] = cos(3*π*x)* cos(3*π*y)* cos(3*π*z)
            face_simple[dir][ijk] = cos(3*π*x)* cos(3*π*y)* cos(3*π*z)
        end
        @inbounds for local_grid in GridField(grid, edge[dir])
            @unpack ijk,x,y,z = local_grid
            edge[dir][ijk] = cos(3*π*x)* cos(3*π*y)* cos(3*π*z)
            edge_simple[dir][ijk] = cos(3*π*x)* cos(3*π*y)* cos(3*π*z)
        end
    end
    for dir in 1:3
        extrap!(face[1], grid, dir)
        extrap!(face[2], grid, dir)
        extrap!(face[3], grid, dir)
        extrap!(edge[1], grid, dir)
        extrap!(edge[2], grid, dir)
        extrap!(edge[3], grid, dir)
        extrap!(cent, grid, dir)
        extrap!(corn, grid, dir)
        extrap_simple!(face_simple[1], grid, dir)
        extrap_simple!(face_simple[2], grid, dir)
        extrap_simple!(face_simple[3], grid, dir)
        extrap_simple!(edge_simple[1], grid, dir)
        extrap_simple!(edge_simple[2], grid, dir)
        extrap_simple!(edge_simple[3], grid, dir)
        extrap_simple!(cent_simple, grid, dir)
        extrap_simple!(corn_simple, grid, dir)
        @test all(face_simple[1] .≈ face[1])
        @test all(face_simple[2] .≈ face[2])
        @test all(face_simple[3] .≈ face[3])
        @test all(edge_simple[1] .≈ edge[1])
        @test all(edge_simple[2] .≈ edge[2])
        @test all(edge_simple[3] .≈ edge[3])
        @test all(cent_simple .≈ cent)
        @test all(corn_simple .≈ corn)
    end
end

#####
##### Staggered 2-field extrapolation
#####

function extrap_x_simple!(f::AbstractField, g::AbstractField,
    grid::Grid, f_dl::Center1D, g_dl::Vertex1D)
    f[ 1 ,:,:] .= 2*g[  2  ,:,:] .- f[  2  ,:,:]
    f[end,:,:] .= 2*g[end-1,:,:] .- f[end-1,:,:]
end
function extrap_x_simple!(f::AbstractField, g::AbstractField,
    grid::Grid, f_dl::Vertex1D, g_dl::Center1D)
    f[ 1 ,:,:] .= 2*g[ 1 ,:,:] .- f[  2  ,:,:]
    f[end,:,:] .= 2*g[end,:,:] .- f[end-1,:,:]
end

@testset "Extrapolations - Staggered 2-field" begin
    FT = Float64
    n, a, b, β = 10, FT(0.0), FT(1.0), FT(1.01)
    grid = Cube(a, b, n; warpfun=Roberts_both, args=(β,))

    cent = CellCenter(grid)
    corn = CellCorner(grid)
    face = CellFace(grid, 1)
    edge = CellEdge(grid, 1)
    cent_simple = CellCenter(grid)
    corn_simple = CellCorner(grid)
    face_simple = CellFace(grid, 1)
    edge_simple = CellEdge(grid, 1)
    cent_out = CellCenter(grid)
    corn_out = CellCorner(grid)
    face_out = CellFace(grid, 1)
    edge_out = CellEdge(grid, 1)
    @inbounds for local_grid in GridField(grid, cent)
        @unpack ijk,x,y,z = local_grid
        cent[ijk] = cos(3*π*x)* cos(3*π*y)* cos(3*π*z)
        cent_simple[ijk] = cos(3*π*x)* cos(3*π*y)* cos(3*π*z)
        cent_out[ijk] = cos(3*π*x)* cos(3*π*y)* cos(3*π*z)
    end
    @inbounds for local_grid in GridField(grid, corn)
        @unpack ijk,x,y,z = local_grid
        corn[ijk] = cos(3*π*x)* cos(3*π*y)* cos(3*π*z)
        corn_simple[ijk] = cos(3*π*x)* cos(3*π*y)* cos(3*π*z)
        corn_out[ijk] = cos(3*π*x)* cos(3*π*y)* cos(3*π*z)
    end
    @inbounds for local_grid in GridField(grid, face)
        @unpack ijk,x,y,z = local_grid
        face[ijk] = cos(3*π*x)* cos(3*π*y)* cos(3*π*z)
        face_simple[ijk] = cos(3*π*x)* cos(3*π*y)* cos(3*π*z)
        face_out[ijk] = cos(3*π*x)* cos(3*π*y)* cos(3*π*z)
    end
    @inbounds for local_grid in GridField(grid, edge)
        @unpack ijk,x,y,z = local_grid
        edge[ijk] = cos(3*π*x)* cos(3*π*y)* cos(3*π*z)
        edge_simple[ijk] = cos(3*π*x)* cos(3*π*y)* cos(3*π*z)
        edge_out[ijk] = cos(3*π*x)* cos(3*π*y)* cos(3*π*z)
    end

    extrap_x_simple!(cent_simple, face, grid, Center1D(), Vertex1D())
    extrap_x_simple!(face_simple, cent, grid, Vertex1D(), Center1D())
    extrap_x_simple!(corn_simple, edge, grid, Vertex1D(), Center1D())
    extrap_x_simple!(edge_simple, corn, grid, Center1D(), Vertex1D())

    extrap!(cent_out, face, grid, 1)
    extrap!(face_out, cent, grid, 1)
    extrap!(corn_out, edge, grid, 1)
    extrap!(edge_out, corn, grid, 1)

    @test all(cent_out .≈ cent_simple)
    @test all(face_out .≈ face_simple)
    @test all(corn_out .≈ corn_simple)
    @test all(edge_out .≈ edge_simple)
end

@testset "Extrapolations - test throws" begin
    FT = Float64
    n, a, b = 10, FT(0.0), FT(1.0)
    grid = Cube(a, b, n)
    face = CellFaces(grid)
    edge = CellEdges(grid)
    for dir in 1:3
        @test_throws MethodError extrap!(face[dir], edge[dir], grid, dir)
        @test_throws MethodError extrap!(edge[dir], face[dir], grid, dir)
    end
end

@testset "Extrapolations - Staggered 2-field convergence" begin
    FT = Float64
    a, b = FT(0.0), FT(1.0)
    nr = 2 .^ (3,4,5,6)
    Δh = zeros(length(nr))
    err_c2v = zeros(length(nr))
    err_v2c = zeros(length(nr))
    for (k,n) in enumerate(2 .^ (3,4,5,6))
        grid = Cube(a, b, n)
        Δh[k] = grid.c[1].n.Δh[1]
        cent = CellCenter(grid)
        corn = CellCorner(grid)
        face = CellFace(grid, 1)
        edge = CellEdge(grid, 1)
        cent_exact = CellCenter(grid)
        corn_exact = CellCorner(grid)
        face_exact = CellFace(grid, 1)
        edge_exact = CellEdge(grid, 1)
        @inbounds for local_grid in GridField(grid,face)
            @unpack ijk,x,y,z = local_grid
            face[ijk] = sin(3*π*x)* sin(3*π*y)* sin(3*π*z)
            face_exact[ijk] = sin(3*π*x)* sin(3*π*y)* sin(3*π*z)
        end
        @inbounds for local_grid in GridField(grid,cent)
            @unpack ijk,x,y,z = local_grid
            cent[ijk] = sin(3*π*x)* sin(3*π*y)* sin(3*π*z)
            cent_exact[ijk] = sin(3*π*x)* sin(3*π*y)* sin(3*π*z)
        end
        extrap!(cent, face, grid, 1)
        extrap!(face, cent, grid, 1)
        err_c2v[k] = norm(cent .- cent_exact) / length(cent)
        err_v2c[k] = norm(face .- face_exact) / length(face)
    end
    conv_c2v = convergence_rate(err_c2v, Δh)
    conv_v2c = convergence_rate(err_v2c, Δh)
    @test conv_c2v[1] ≈ 4.093446045776282
    @test conv_c2v[2] ≈ 2.6665448437905694
    @test conv_c2v[3] ≈ 3.1919109381432618

    @test conv_v2c[1] ≈ 4.5665684649141545
    @test conv_v2c[2] ≈ 4.719101505009674
    @test conv_v2c[3] ≈ 4.848188053541373
end


text:
export AbstractUpdateRule, TikTakRule
export new_guess, weight_parameter, weight_sequence

"""
	$(SIGNATURES)

Abstract type for generating new guesses from the history of already computed points and from the set of starting points provided by the user.
"""
abstract type AbstractUpdateRule end

"""
	$(SIGNATURES)

Generate a sequence of weights for `jMax` points. Useful for tuning the parameters of the update rule.
"""
function weight_sequence(ur :: AbstractUpdateRule, jMax :: Integer)
    return [weight_parameter(ur, j)  for j = 1 : jMax]
end

"""
	$(SIGNATURES)

For generating new points. Defaults based on

Arnoud, Antoine, Fatih Guvenen, and Tatjana Kleineberg. 2019. “Benchmarking Global Optimizers.” Working Paper 26340. National Bureau of Economic Research. https://doi.org/10.3386/w26340.
"""
Base.@kwdef struct TikTakRule <: AbstractUpdateRule
    nPoints :: Int
    θ_min :: Float64 = 0.1
    θ_max :: Float64 = 0.995
    θ_pow :: Float64 = 0.5
    # Bounds for elements of the guess vector.
    lb :: Float64 = 1.0
    ub :: Float64 = 2.0
end

"""
	$(SIGNATURES)

Make a new guess vector from the history of points computed so far and from the set of user provided starting points. Implements `TikTak`.
"""
function new_guess(ur :: TikTakRule, h :: History, j :: Int)
    if j == 1
        guessV = first(h.guessesV);
    else
        wt = weight_parameter(ur, j);
        pBest = best_point(h);
        guessV = wt .* pBest.solnV .+ (1 - wt) .* h.guessesV[j];
    end
    clamp!(guessV, ur.lb, ur.ub);
    return guessV
end

"""
	$(SIGNATURES)

`TikTak` weight parameter.
"""
function weight_parameter(ur :: TikTakRule, j :: Integer)
    return clamp((j / ur.nPoints) ^ ur.θ_pow,
        ur.θ_min, ur.θ_max)
end

# ---------------


text:
draw(PNG("negative-landscape.png", 140, 100),
     Axes(negative=true),
     Paper("lightgrey"),
     # default Pen is 0.02.  use 0.04 here so consistent with bottom
     # left image which is twice as dark because the user coordinates
     # are half as large.
     Pen(0.04)) do
    move(0, 0)
    line(1, 0)
    move(0, 0)
    line(0, 1)
end


text:
include("models.jl")

function test_bernoulli()
    model_bernoulli = bernoulli_example()
    graph_bernoulli = trackdependencies(model_bernoulli)
    @testdependencies(model_bernoulli, w, p, x)
    cond_bernoulli_p = StaticConditional(model_bernoulli, :p)
    @test_nothrow sample(model_bernoulli, Gibbs(cond_bernoulli_p, MH(:w)), 500)

    # Analytic tests
    w = graph_bernoulli[4].value
    p = graph_bernoulli[6].value
    x = graph_bernoulli[2].value
    p_1 = w[1] * pdf(Bernoulli(0.3), x)
    p_2 = w[2] * pdf(Bernoulli(0.7), x)
    Z = p_1 + p_2

    # 𝓅(p | w, x) ∝ 𝓅(p | w) * 𝓅(x | p)
    analytic_conditional = DiscreteNonParametric([0.3, 0.7], [p_1 / Z, p_2 / Z])
    @info "Bernoulli analytic" analytic_conditional
    θ = AutoGibbs.sampled_values(graph_bernoulli)

    local calculated_conditional
    @test_nothrow calculated_conditional = conditionals(graph_bernoulli, @varname(p))[@varname(p)]
    @info "Bernoulli calculated" calculated_conditional(θ)

    @test issimilar(calculated_conditional(θ), analytic_conditional)
end

function test_gmm()
    model_gmm = gmm_example()
    graph_gmm = trackdependencies(model_gmm)
    @testdependencies(model_gmm, μ, w, z, x[1], x[2], x[3])
    cond_gmm_z = StaticConditional(model_gmm, :z)
    @test_nothrow sample(model_gmm, Gibbs(cond_gmm_z, MH(:w, :μ)), 500)
    @test_nothrow sample(model_gmm, Gibbs(cond_gmm_z, HMC(0.01, 10, :w, :μ)), 500)

    # Analytic tests
    x = graph_gmm[2].value
    μ = graph_gmm[7].value
    w = graph_gmm[10].value
    z = graph_gmm[14].value
    p_1 = w[1] .* pdf.(Normal(μ[1], s2_gmm), x)
    p_2 = w[2] .* pdf.(Normal(μ[2], s2_gmm), x)
    (Z1, Z2, Z3) = p_1 .+ p_2

    # 𝓅(zᵢ | μ, w, x, z₋ᵢ) ∝ 𝓅(zᵢ | w) * 𝓅(xᵢ | zᵢ, μ)
    analytic_conditionals = [@varname(z[1]) => Categorical([p_1[1], p_2[1]] ./ Z1),
                             @varname(z[2]) => Categorical([p_1[2], p_2[2]] ./ Z2),
                             @varname(z[3]) => Categorical([p_1[3], p_2[3]] ./ Z3)]
    @info "GMM analytic" analytic_conditionals
    θ = AutoGibbs.sampled_values(graph_gmm)

    local calculated_conditional
    @test_nothrow calculated_conditional = conditionals(graph_gmm, @varname(z))[@varname(z)]
    @info "GMM calculated" calculated_conditional(θ)

    @test issimilar(calculated_conditional(θ), Product([D for (vn, D) in analytic_conditionals]))
end

function test_gmm_loopy()
    model_gmm_loopy = gmm_loopy_example()
    graph_gmm_loopy = trackdependencies(model_gmm_loopy)
    @testdependencies(model_gmm_loopy, μ[1], μ[2], w, z[1], z[2], z[3], x[1], x[2], x[3])
    cond_gmm_loopy_z = StaticConditional(model_gmm_loopy, :z)
    @test_nothrow sample(model_gmm_loopy, Gibbs(cond_gmm_loopy_z, MH(:w, :μ)), 500)
    @test_nothrow sample(model_gmm_loopy, Gibbs(cond_gmm_loopy_z, HMC(0.01, 10, :w, :μ)), 500)

    # Analytic tests
    x = graph_gmm_loopy[2].value
    μ = [graph_gmm_loopy[19].value, graph_gmm_loopy[28].value]
    w = graph_gmm_loopy[32].value
    z = [graph_gmm_loopy[44].value, graph_gmm_loopy[63].value, graph_gmm_loopy[82].value]
    p_1 = w[1] .* pdf.(Normal(μ[1], s2_gmm), x)
    p_2 = w[2] .* pdf.(Normal(μ[2], s2_gmm), x)
    (Z1, Z2, Z3) = p_1 .+ p_2

    # 𝓅(zᵢ | μ, w, x, z₋ᵢ) ∝ 𝓅(zᵢ | w) * 𝓅(xᵢ | zᵢ, μ)
    analytic_conditionals = [@varname(z[1]) => Categorical([p_1[1], p_2[1]] ./ Z1),
                             @varname(z[2]) => Categorical([p_1[2], p_2[2]] ./ Z2),
                             @varname(z[3]) => Categorical([p_1[3], p_2[3]] ./ Z3)]
    @info "Loopy GMM analytic" analytic_conditionals
    θ = AutoGibbs.sampled_values(graph_gmm_loopy)

    local calculated_conditionals
    @test_nothrow calculated_conditionals = conditionals(graph_gmm_loopy, @varname(z))
    @info "Loopy GMM calculated" Dict(vn => cond(θ) for (vn, cond) in calculated_conditionals)

    for (vn, analytic_conditional) in analytic_conditionals
        @test issimilar(calculated_conditionals[vn](θ), analytic_conditional)
    end
end

function test_gmm_shifted()
    model_gmm_shifted = gmm_shifted_example()
    graph_gmm_shifted = trackdependencies(model_gmm_shifted)
    @testdependencies(model_gmm_shifted, μ[1], μ[2], w, z[1], z[2], z[3], x[1], x[2], x[3])
    cond_gmm_shifted_z = StaticConditional(model_gmm_shifted, :z)
    @test_nothrow sample(model_gmm_shifted, Gibbs(cond_gmm_shifted_z, MH(:w, :μ)), 500)
    @test_nothrow sample(model_gmm_shifted, Gibbs(cond_gmm_shifted_z, HMC(0.01, 10, :w, :μ)), 500)
end

function test_hmm()
    model_hmm = hmm_example()
    graph_hmm = trackdependencies(model_hmm)
    @testdependencies(model_hmm, T[1], T[2], m[1], m[2], s[1], s[2], s[3], x[1], x[2], x[3])
    cond_hmm_s = StaticConditional(model_hmm, :s)
    @test_nothrow sample(model_hmm, Gibbs(cond_hmm_s, MH(:m, :T)), 500)
    @test_nothrow sample(model_hmm, Gibbs(cond_hmm_s, HMC(0.01, 10, :m, :T)), 500)

    # Analytic tests
    x = graph_hmm[2].value
    T = [graph_hmm[24].value, graph_hmm[40].value]
    m = [graph_hmm[30].value, graph_hmm[46].value]
    s1 = graph_hmm[50].value
    s2 = graph_hmm[70].value
    s3 = graph_hmm[91].value
    D_obs_1 = Normal(m[1], s2_hmm)
    D_obs_2 = Normal(m[2], s2_hmm)
    p_s1_1 = pdf(Categorical(2), 1) * pdf(Categorical(T[1]), s2) * pdf(D_obs_1, x[1])
    p_s1_2 = pdf(Categorical(2), 2) * pdf(Categorical(T[2]), s2) * pdf(D_obs_2, x[1])
    p_s2_1 = pdf(Categorical(T[s1]), 1) * pdf(Categorical(T[1]), s3) * pdf(D_obs_1, x[2])
    p_s2_2 = pdf(Categorical(T[s1]), 2) * pdf(Categorical(T[2]), s3) * pdf(D_obs_2, x[2])
    p_s3_1 = pdf(Categorical(T[s2]), 1) * pdf(D_obs_1, x[3])
    p_s3_2 = pdf(Categorical(T[s2]), 2) * pdf(D_obs_2, x[3])
    Z_1 = p_s1_1 + p_s1_2
    Z_2 = p_s2_1 + p_s2_2
    Z_3 = p_s3_1 + p_s3_2

    # 𝓅(s₁ | T, m, s₋₁, x) ∝ 𝓅(s₁) 𝓅(s₂ | s₁, T) 𝓅(x₁ | s₁, m)
    #  𝓅(sᵢ | T, m, s₋ᵢ, x) ∝  𝓅(sᵢ | sᵢ₋₁, T) 𝓅(sᵢ₊₁ | sᵢ, T) 𝓅(xᵢ | sᵢ, m) (for i ≥ 2)
    analytic_conditionals = [@varname(s[1]) => Categorical([p_s1_1, p_s1_2] ./ Z_1),
                             @varname(s[2]) => Categorical([p_s2_1, p_s2_2] ./ Z_2),
                             @varname(s[3]) => Categorical([p_s3_1, p_s3_2] ./ Z_3)]
    θ = AutoGibbs.sampled_values(graph_hmm)
    @info "HMM analytic" analytic_conditionals

    local calculated_conditionals
    @test_nothrow calculated_conditionals = conditionals(graph_hmm, @varname(s))
    @info "HMM calculated" Dict(vn => cond(θ) for (vn, cond) in calculated_conditionals)

    for (vn, analytic_conditional) in analytic_conditionals
        # @show vn => probs(calculated_conditionals[vn]), probs(analytic_conditional)
        @test issimilar(calculated_conditionals[vn](θ), analytic_conditional)
    end
end

function test_imm_stick()
    model_imm_stick = imm_stick_example()
    graph_imm_stick = trackdependencies(model_imm_stick)

    # we leave out the μs, because there might be 1--3 of them
    @testdependencies(model_imm_stick,
                      v, # 10 - 1 sticks
                      μ, # 10 cluster centers
                      # z, # 9 data points
                      y[1], y[2], y[3], y[4], y[5], y[6], y[7], y[8], y[9])
    cond_imm_stick = StaticConditional(model_imm_stick, :z)
    @test_nothrow sample(model_imm_stick, Gibbs(cond_imm_stick, MH(:μ, :v)), 500)
    @test_nothrow sample(model_imm_stick, Gibbs(cond_imm_stick, HMC(0.01, 10, :μ, :v)), 500)

    # Analytic tests
    y = graph_imm_stick[2].value
    α = graph_imm_stick[4].value
    K = graph_imm_stick[6].value
    N = graph_imm_stick[7].value
    v = graph_imm_stick[12].value
    # z = graph_imm_stick[16].value
    z = [v.value for v in values(graph_imm_stick.statements)
         if v isa AutoGibbs.Assumption && DynamicPPL.subsumes(@varname(z), v.vn)]
    # μ = graph_imm_stick[19].value
    μ = (v.value for v in values(graph_imm_stick.statements)
         if v isa AutoGibbs.Assumption && DynamicPPL.subsumes(@varname(μ), v.vn)) |> first

    D_w = Categorical(stickbreak(v))
    analytic_conditionals = map(1:N) do n
        p̃ = [exp(logpdf(D_w, z) + logpdf(Normal(μ[z], 1.0), y[n])) for z in support(D_w)]
        @varname(z[n]) => Categorical(p̃ ./ sum(p̃))
    end
    # analytic_conditional = Product([D for (vn, D) in analytic_conditionals])
    # @info "stick-breaking IMM analytic conditionals" analytic_conditionals

    θ = AutoGibbs.sampled_values(graph_imm_stick)

    local calculated_conditionals
    @test_nothrow calculated_conditionals = conditionals(graph_imm_stick, @varname(z))
    # @info "stick-breaking IMM calculated conditionals" Dict(
        # vn => cond(θ) for (vn, cond) in calculated_conditionals)

    for (vn, analytic_conditional) in analytic_conditionals
        # @show vn => probs(calculated_conditionals[vn]), probs(analytic_conditional)
        @test issimilar(calculated_conditionals[vn](θ), analytic_conditional)
    end
end

function test_changepoint()
    model_changepoint = changepoint([1.1, 0.9, 0.2])
    graph_changepoint = trackdependencies(model_changepoint)
    @testdependencies(model_changepoint, λ₁, λ₂, τ, y[1], y[2], y[3])
    @test_nothrow sample(model_changepoint, Gibbs(AutoConditional(:τ), MH(:λ₁, :λ₂)), 20)
end

# function test_reverseps()
#     model_reverse_deps = reverse_deps([0.1, -0.2])
#     graph_reverse_deps = trackdependencies(model_reverse_deps)
#     @testdependencies(model_reverse_deps, m[1], m[2], x)
# end

##########################################################################
#########################################################################
### TEST TOGGLES

test_bernoulli()
test_gmm()
test_gmm_loopy()
test_gmm_shifted()
test_hmm()
test_imm_stick()
test_changepoint()


text:
abstract type AbstractDEBCore end

@columns struct DEBCore{MoMo,GMo} <: AbstractDEBCore
    y_V_E::MoMo  | 0.7   | _        | Beta(2.0, 2.0)  | [0.0, 1.0]   | _     | "From reserve to structure"
    y_E_EC::MoMo | 0.7   | _        | Gamma(2.0, 2.0) | [1e-6, 1.0]  | false | "From C-reserve to reserve, using nitrate"
    y_E_EN::MoMo | 30.0  | _        | Gamma(2.0, 2.0) | [1.0, 50.0]  | false | "From N-reserve to reserve"
    n_N_V::MoMo  | 0.03  | _        | Gamma(2.0, 2.0) | [0.0, 0.1]   | _     | "N/C in structure"
    n_N_E::MoMo  | 0.025 | _        | Gamma(2.0, 2.0) | [0.0, 0.1]   | _     | "N/C in reserve"
    w_V::GMo     | 25.0  | g*mol^-1 | Gamma(2.0, 2.0) | [15.0, 40.0] | _     | "Mol-weight of shoot structure"
    # w_N::GMo     | 25.0  | g*mol^-1 | Gamma(2.0, 2.0) | [15.0, 40.0] | _     | "Mol-weight of shoot N-reserve"
    # w_C::GMo   | 25.0  | g*mol^-1 | Gamma(2.0, 2.0) | [12.0, 40.0] | _     | "Mol-weight of shoot C-reserve"
    # w_E::GMo   | 25.0  | g*mol^-1 | Gamma(2.0, 2.0) | [15.0, 40.0] | _     | "Mol-weight of shoot reserve"
end

abstract type AbstractParams end

" Model parameters that vary between organs "
@selectable @flattenable @default_kw struct Params{As,Sh,Al,Ma,Tr,Re,Ge,Pr} <: AbstractParams
    # Field               | Default                | _     | Selectable Types
    name::Symbol          | :organ                 | false | _
    rate_formula          | FZeroRate()            | _     | _
    assimilation_pars::As | ConstantCAssim()       | _     | Union{Nothing,AbstractAssim}
    shape_pars::Sh        | Plantmorph()           | _     | Union{Nothing,AbstractShape}
    allometry_pars::Al    | nothing                | _     | Union{Nothing,AbstractAllometry}
    maturity_pars::Ma     | nothing                | _     | Union{Nothing,AbstractMaturity}
    trans_pars::Tr        | nothing                | _     | Union{Nothing,AbstractTranslocation}
    rejection_pars::Re    | LosslessRejection()    | _     | AbstractRejection
    germination_pars::Ge  | ThresholdGermination() | _     | Union{Nothing,AbstractGermination}
    production_pars::Pr   | Production()           | _     | Union{Nothing,AbstractProduction}
end

#    %   W   Rel n atoms
# C  45  12  30857
# H  6   1   58600
# O  45  16  27000
# N  2   14   1028
# K  1   39    246
# X  1   40    300ish

" Model parameters shared between organs "
@selectable @udefault_kw struct SharedParams{SU,Co,Fe,Te,Ca,Mt}
    # Field               | Default                   | Selectable Types
    su_pars::SU           | ParallelComplementarySU() | AbstractSynthesizingUnit
    core_pars::Co         | DEBCore()                 | _
    feedback_pars::Fe     | nothing                   | Union{Nothing,AbstractStateFeedback}
    tempcorr_pars::Te     | nothing                   | Union{Nothing,AbstractTemperatureCorrection}
    catabolism_pars::Ca   | CatabolismCN()            | AbstractCatabolism
    maintenance_pars::Mt  | Maintenance()             | AbstractMaintenance
end

###########################################################################################
# Variables

" Model variables "
@plottable @units @udefault_kw struct Vars{F,MoMoD,C,M,T}
    shape::F             | [0.0]   | _                 | _
    rate::MoMoD          | [0.0]   | mol*mol^-1*d^-1   | _
    θE::F                | [0.0]   | _                 | _
    temp::C              | [25.0]  | K                 | _
    tempcorrection::F    | [1.0]   | _                 | _
    height::M            | [0.0]   | m                 | _
    t::T                 | [1]     | _                 | false
end

build_vars(vars, time) = begin
    len = length(time)
    len == length(vars.rate) && return vars

    fields = []
    for fname in fieldnames(typeof(vars))
        ft = fieldtype(typeof(vars), fname)
        if ft <: AbstractArray
            push!(fields, fill(getfield(vars, fname)[1], len))
        else
            push!(fields, getfield(vars, fname))
        end
    end
    typeof(vars).name.wrapper(fields...)
end

###########################################################################################
# Organs and Organisms

abstract type AbstractOrgan end

" Basic model components. For a plants, organs might be roots, stem and leaves "
struct Organ{P,S,V,F,F1} <: AbstractOrgan
    params::P
    shared::S
    vars::V
    J::F
    J1::F1
end

" update vars and flux to record for time t "
define_organs(o, t) = define_organs(o.params, o.shared, o.records, t)
define_organs(params::Tuple{P,Vararg}, shared, records::Tuple{R,Vararg}, t) where {P,R} = begin
    rec = records[1]
    t = calc_tstep(rec.vars, t)
    rec.vars.t[1] = t
    vJ = view(rec.J, :, :, t)
    vJ1 = view(rec.J1, :, :, t)
    J = LArray{Tuple{STATE,TRANS}}(vJ)
    J1 = LArray{Tuple{STATE1,TRANS1}}(vJ1)

    organ = Organ(params[1], shared, rec.vars, J, J1)
    (organ, define_organs(tail(params), shared, tail(records), t)...)
end
define_organs(params::Tuple{}, shared, records::Tuple{}, t) = ()

"Records of mutable variables and flux for ploting and analysis"
@plottable struct Records{V,F,F1}
    vars::V | true
    J::F    | false
    J1::F1  | false
end
"Constructor for records. Arrays use the length of the current timespan"
Records(params, vars, time, val, typ) = begin
    vars = build_vars(vars, time)
    J = build_J(val, typ, time)
    J1 = build_J1(val, typ, time)
    Records(vars, J, J1)
end

build_J(one_flux, T, time) = build_flux(one_flux, T, STATE, TRANS, time)

build_J1(one_flux, T, time) = build_flux(one_flux, T, STATE1, TRANS1, time)

build_flux(one_flux, T, x, y, time) = zeros(typeof(one_flux), length(x), length(y), length(time))

abstract type AbstractOrganism end

"An organism, made up of organs"
@flattenable mutable struct Plant{P,S,R,E,ES,D} <: AbstractOrganism
    params::P             | true
    shared::S             | true
    records::R            | false
    environment::E        | false
    environment_start::ES | false
    dead::D               | false
end

"Outer construtor for defaults"
Plant(; params = (ShootParamsCN(), RootParamsCN()),
        vars = (Vars(), Vars()),
        shared = SharedParams(),
        records = nothing,
        environment = nothing,
        time = 0hr:1hr:8760hr,
        environment_start = setindex!(Array{typeof(1hr),0}(undef), 1hr),
        dead = Ref(false)
      ) = begin
    if records == nothing
        records = []
        for i = 1:length(params)
            push!(records, Records(params[i], vars[i], time, 1.0mol/hr, typeof(1.0mol/hr)))
        end
        records = (records...,)
    end
    Plant(params, shared, records, environment, environment_start, dead)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "DeldirBuilder"
version = v"0.1.15"

# Collection of sources required to build DeldirBuilder
sources = [
    "https://cran.r-project.org/src/contrib/Archive/deldir/deldir_0.1-9.tar.gz" =>
    "090aba6747efc77424a44bf4aaa229dfc01fff7752720689bb70cd861df61f6a",
    "https://cran.r-project.org/src/contrib/Archive/deldir/deldir_0.1-15.tar.gz" =>
    "571011c208829f47ecd9f92a19fd94a8eb59de5a2645ab8c62e73926ade30710",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/deldir/src

for f in *.f; do
    ${FC} -fPIC -O2 -pipe -g -c "${f}" -o "$(basename "${f}" .f).o"
done

if [[ ${target} == *-mingw32 ]]; then
    libdir="bin"
else
    libdir="lib"
fi

mkdir -p "${prefix}/${libdir}"
${CC} -shared -o ${prefix}/${libdir}/libdeldir.${dlext} *.o
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:powerpc64le, libc=:glibc),
    Linux(:i686, libc=:musl),
    Linux(:x86_64, libc=:musl),
    Linux(:aarch64, libc=:musl),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf),
    MacOS(:x86_64),
    Windows(:i686),
    Windows(:x86_64)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libdeldir", :libdeldir)
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: robertdj/DeldirBuilder

max_issues_repo_name: robertdj/DeldirBuilder

max_forks_repo_name: robertdj/DeldirBuilder


text:
using DifferentialEquations
using Plots, LaTeXStrings
using Markdown
using Printf

# Equations
"""Single pendulum with friction"""
spwf = (θ, ω, g, L, μ) -> -μ * ω - g / L * sin(θ)

"""Double pednulum (withou friction)"""
function dp(θ₁, θ₂, ω₁, ω₂, m₁, m₂, g, L₁, L₂)
    dω₁ = (
        -g*m₂*(sin(2*θ₂-θ₁)-sin(θ₁)) - L₁*m₂*ω₁^2*sin(2*θ₂-2*θ₁)
        - 2*L₂*m₂*ω₂^2*sin(θ₂-θ₁) + 2*g*m₁*sin(θ₁)
        ) / (L₁*m₂*(cos(2*(θ₂-θ₁))-1) - 2*L₁*m₁)
    dω₂ = (
        L₂*m₂*ω₂^2*sin(2*θ₂-2*θ₁) + 2*L₁*m₂*ω₁^2*sin(θ₂-θ₁) + 2*L₁*m₁*ω₁^2*sin(θ₂-θ₁) +
        g*m₂*(sin(θ₂-2*θ₁)+sin(θ₂)) + g*m₁*(sin(θ₂-2*θ₁)+sin(θ₂))
        ) / (L₂*m₂*(cos(2*(θ₂-θ₁))-1) - 2*L₂*m₁)
    return (dω₁, dω₂)
end

"""Double pednulum with friction"""
function dpwf(θ₁, θ₂, ω₁, ω₂, m₁, m₂, g, L₁, L₂, μ₁, μ₂)
    dω₁ = (
        L₁*ω₁*(cos(2*(θ₂-θ₁))-1)*μ₂ - 2*L₁*ω₁*μ₁ + g*m₂*(sin(2*θ₂-θ₁)-sin(θ₁)) +
        L₁*m₂*ω₁^2*sin(2*θ₂-2*θ₁) + 2*L₂*m₂*ω₂^2*sin(θ₂-θ₁) - 2*g*m₁*sin(θ₁)
        ) / (2*L₁*m₁ - L₁*m₂*(cos(2*(θ₂-θ₁))-1))
    dω₂ = (
        L₂*m₂*ω₂*(cos(2*(θ₂-θ₁))-1)*μ₂ - 2*L₁*m₁*ω₁*cos(θ₂-θ₁)*μ₂ - 2*L₂*m₁*ω₂*μ₂ +
        2*L₁*m₂*ω₁*cos(θ₂-θ₁)*μ₁ - L₂*m₂^2*ω₂^2*sin(2*θ₂-2*θ₁) - 2*L₁*m₂^2*ω₁^2*sin(θ₂-θ₁) -
        2*L₁*m₁*m₂*ω₁^2*sin(θ₂-θ₁) - g*m₂^2*(sin(θ₂-2*θ₁)+sin(θ₂)) - g*m₁*m₂*(sin(θ₂-2*θ₁)+sin(θ₂))
        ) / (2*L₂*m₁*m₂ - L₂*m₂^2*(cos(2*(θ₂-θ₁))-1))
    return (dω₁, dω₂)
end

# Parameters
g = 9.8
θ₁_, θ₂_ = [0 0]
ω₁_, ω₂_ = [4 0]
m₁, m₂ = [1 4]
L₁, L₂ = [5 10]
μ₁, μ₂ = 0.1 .* [m₁ m₂]
u₀ = [θ₁_, ω₁_]
u₀_double = [θ₁_, θ₂_, ω₁_, ω₂_]
params = [g, L₁, μ₁]
params_double = [m₁, m₂, g, L₁, L₂]
params_double_wf = [m₁, m₂, g, L₁, L₂, μ₁, μ₂]
tspan = (0.0, 100.0)

# Helper functions
x = (L, θ) -> L * sin(θ)
y = (L, θ) -> - L * cos(θ)

function vec_to_list(vecvec, dim)
    vecvec[dim]
end

function pendulum_problem!(du, u, params, t)
    g, L, μ = params
    θ, ω = u
    du[1] = ω
    du[2] = spwf(u..., params...)
end

# Single pendulum solution
prob = ODEProblem(pendulum_problem!, u₀, tspan, params)
sol_single_wf = solve(prob, reltol=1e-5, saveat=10e-3)
plot(sol_single_wf, linewidth=1.5, layout=(2,1), label=["θ [rad]" "ω [rad/s]"], xaxis="t")

# Vector field
θ_bounds = 6*pi
ω_bounds = 2*pi
step = [ pi/5, pi/5]
kh = 0.45
mesh = [[i, j] for i=-θ_bounds:step[1]:θ_bounds for j=-ω_bounds:step[2]:ω_bounds]
slopes = [spwf(point..., params...) for point in mesh]
θf, ωf = @. (vec_to_list(mesh, 1), vec_to_list(mesh, 2))
lens = @. 1 / sqrt(ωf ^ 2 + slopes ^ 2) * kh

θ, ω = @. (vec_to_list(sol_single_wf.u, 1), vec_to_list(sol_single_wf.u, 2));

xticks_str = ["0", "π/2", "π", "3π/2", "2π", "5π/2", "3π", "7π/2", "", "9π/2", "5π"]
fieldp = quiver(θf, ωf, quiver=(lens .* ωf, lens .* slopes), line=(1.2, 0.2))
plot!(fieldp, θ, ω, line=(2), label="θ'(t) = θ(t)", legend=:topright)
plot!(fieldp, framestyle=:origin, xticks = (collect(0:pi/2:5*pi), xticks_str), aspect_ratio=1, dpi=300)
plot!(fieldp,  xaxis=("θ", (-1*pi/2,11*pi/2)), yaxis=("ω = θ'", (-3*pi/3,4*pi/3)), size=(1000, 400))

# Single pendulum animation
@userplot pendulumanim
@recipe function fpend!(oa::pendulumanim)
    x, y, m, title_time = oa.args
    title --> "Pendulum plot :: [$(@sprintf("%2.1f", title_time))%]"
    xlims --> (-L₁-1, L₁+1)
    ylims --> (-L₁-1, L₁+1)
    markersize --> 5 * m
    markercolor --> "black"
    seriestype --> :scatter
    legend --> :none
    framestyle --> :grid
    aspect_ratio --> 1
    dpi --> 100
    [x], [y]
end

function pendulum_example()

    # time frame
    t = 0:0.02:50
    k = 4
    n = length(t) / k |> floor |> Int

    # animation
    anim = @animate for i in 1:n
        x_ball = x(L₁, θ[k*i])
        y_ball = y(L₁, θ[k*i])
        spend = plot([0, x_ball], [0, y_ball])
        pendulumanim!(spend, x_ball, y_ball, m₁, i/n * 100)
    end

    gif(anim, "animations/single_pendulum_with_friction.gif", fps=50)
end

# Double pendulum solution
function double_pendulum_problem!(du, u, params, t)
    m₁, m₂, g, L₁, L₂ = params
    θ₁, θ₂, ω₁, ω₂ = u
    du[1] = ω₁  # du[3]
    du[2] = ω₂  # du[4]
    du[3], du[4] = dp(u..., params...)
end

prob_double = ODEProblem(double_pendulum_problem!, u₀_double, tspan, params_double)
sol_double = solve(prob_double, reltol=1e-5, saveat=10e-3)
double_anim = plot(sol_double,
                   linewidth=1.5,
                   layout=(2,2),
                   label=["θ₁ [rad]" "θ₂ [rad]" "ω₁ [rad/s]" "ω₂ [rad/s]"],
                   xaxis="t")

# Double pednulum with friction solution
function double_pendulum_wf!(du, u, params, t)
    m₁, m₂, g, L₁, L₂, μ₁, μ₂ = params
    θ₁, θ₂, ω₁, ω₂ = u
    du[1] = ω₁  # du[3]
    du[2] = ω₂  # du[4]
    du[3], du[4] = dpwf(u..., params...)
end

prob_double_wf = ODEProblem(double_pendulum_wf!, u₀_double, tspan, params_double_wf)
sol_double_wf = solve(prob_double_wf, reltol=1e-5, saveat=10e-3)

double_wf_anim = plot(sol_double_wf,
                      linewidth=1.5, layout=(2,2),
                      label=["θ₁ [rad]" "θ₂ [rad]" "ω₁ [rad/s]" "ω₂ [rad/s]"],
                      xaxis="t")

# Double pendulum animation
@userplot doublependulumanim
@recipe function ffdpend!(oa::doublependulumanim)
    x, y, m, title_time = oa.args
    title --> "Double pendulum plot :: [$(@sprintf("%2.1f", title_time))%]"
    xaxis --> ("x", (-L₁-L₂-3*m, L₁+L₂+3*m))
    yaxis --> ("y", (-L₁-L₂-3*m, L₁+L₂+3*m))
    markersize --> 3 * m
    markercolor --> "black"
    seriestype --> :scatter
    legend --> :none
    framestyle --> :grid
    aspect_ratio --> 1
    dpi --> 100
    [x], [y]
end

function double_pendulum_animation(solution, filename)

    θ₁, θ₂ = @. (vec_to_list(solution.u, 1), vec_to_list(solution.u, 2))

    # time frame
    t = 0:0.02:100
    k = 4
    n = length(t) / k |> floor |> Int

    # animation
    anim = @animate for i in 1:n
        x_first, y_first = x(L₁, θ₁[k*i]), y(L₁, θ₁[k*i])
        x_second, y_second = x(L₂, θ₂[k*i]), y(L₂, θ₂[k*i])
        plot([0, x_first], [0, y_first])
        plot!([x_first, x_first + x_second], [y_first, y_first + y_second])
        doublependulumanim!(x_first, y_first, m₁, i/n * 100)
        doublependulumanim!(x_first + x_second, y_first + y_second, m₂, i/n * 100)
    end

    gif(anim, "animations/$filename.gif", fps=50)

end

pendulum_example()
double_pendulum_animation(sol_double, "double_pendulum_no_friction")
double_pendulum_animation(sol_double_wf, "double_pendulum_with_friction")

max_stars_repo_name: FrightenedFox/double-pendulum

max_issues_repo_name: FrightenedFox/double-pendulum

max_forks_repo_name: FrightenedFox/double-pendulum


text:
using Random

include("load_data.jl")

function surrogate_randomize_per_train(
        ets::encounter_trains{I,T},
        seed::Int
    ) where {I,T}
    rng = MersenneTwister(seed)
    times_surrogate = similar(trains(ets))
    for (i,train) in enumerate(trains(ets))
        times_surrogate[i] = sort(rand(rng, 0:timestep(ets):duration(ets), length(train)))
    end
    return encounter_trains(times_surrogate, ids(ets), duration(ets), timestep(ets))
end

function surrogate_randomize_all(
        ets::encounter_trains{I,T},
        seed::Int
    ) where {I,T}
    rng = MersenneTwister(seed)

    #create empty trains
    times_surrogate = [T[] for i in 1:length(ets)]
    for i in 1:sum(length.(trains(ets)))
        rand_train = rand(rng, 1:length(ets))
        rand_time  = rand(rng, 0:timestep(ets):duration(ets))
        push!(times_surrogate[rand_train], rand_time)
    end
    return encounter_trains(times_surrogate, ids(ets), duration(ets), timestep(ets))
end


text:
module ch1

import StatsBase; sb = StatsBase

# Page 4
count = 600_000
sample = [1:6...]

x = sb.sample(sample, count, replace=true)
P = sum(x.==3)/count

@show P

# Page 5
wins = zeros(Int32, 1, count)
wins .= 0
for ind = 1:length(wins)
    d = sb.sample(sample, 2)
    if sum(d) == 7 || sum(d) == 11
        wins[ind] = 1
    end
end
wins_count = sum(wins)

@show wins_count
@show wins_count/count

# Page 7
# PDF: Probability Density Function
using Plotly
y = [0,1]
p = [1,1]
line = scatter(x=y, y=p, mode="lines") # "lines+markers"
layout = Layout(
    title="PDF",
    # mirror for the bounding box
    xaxis=attr(title="y", linecolor="black", mirror=true),
    yaxis=attr(title="p(y)", range=[0,1.1], linecolor="black", mirror=true)
)
pt = plot([line], layout)
#display(pt)

# Page 10
import RDatasets; rdata = RDatasets
import CSV
import DataFrames; DF = DataFrames
import Distributions; Dist=Distributions
#using StringEncodings

if !isfile("discoveries.csv")
    download("https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/master/csv/datasets/discoveries.csv",
    "discoveries.csv")
end
#med = rdata.dataset("datasets", "discoveries")
# frame = df.DataFrame()
#frame = CSV.read("discoveries.csv", df)
df = CSV.File(open(read, "discoveries.csv")) |> DF.DataFrame
@show names(df)
#@show frame

# Histogram
hist = histogram(x=df.value, histnorm="probability")
# Poisson distribution
d = Dist.Poisson(3.1)
x = range(0,12,length=13)
line = scatter(x=x,y=Dist.pdf.(d,x))
# Plot
layout = Layout(title="Discoveries")
pt = plot([hist, line], layout)
#display(pt)

# 1.3.1 The Binomial Distribution

# Page 14
# θ is a parameter
# Θ is set of θ
# Binomial

# Page 16
d = Dist.Binomial(3,0.1)
x = range(0,3,length=4)
scat = scatter(x=x,y=Dist.pdf.(d,x), mode="markers",
    #opacity=0.1
    marker=attr(size=15, opacity=1, color="transparent",
        line=attr(color="black", width=2)
    )
)
layout = Layout(title="N=3, p=0.1")
pt = plot([scat], layout)
#display(pt)
savefig(pt, "binomial.png", format="png")
savefig(pt, "binomial.html", format="html")

scat_list = Vector{AbstractTrace}()
d = Dist.Binomial(3,0.1)
x = range(0,3,length=4)
scat = scatter(x=x,y=Dist.pdf.(d,x), mode="markers",
    #opacity=0.1
    marker=attr(size=15, opacity=1, color="transparent",
        line=attr(color="black", width=2)
    )
)
push!(scat_list,scat)

d = Dist.Binomial(3,0.5)
scat = scatter(x=x,y=Dist.pdf.(d,x), mode="markers",
    xaxis="x2", yaxis="y2",
    #opacity=0.1
    marker=attr(size=15, opacity=1, color="transparent",
        line=attr(color="black", width=2)
    )
)
push!(scat_list, scat)
# show how to plot multiple traces in one subplot
line = scatter(x=x, y=Dist.pdf.(d,x), mode="lines",
    xaxis="x2", yaxis="y2",
    line=attr(color="black", width=2)
)
push!(scat_list, line)

d = Dist.Binomial(3,0.9)
scat = scatter(x=x,y=Dist.pdf.(d,x), mode="markers",
    xaxis="x3", yaxis="y3",
    #opacity=0.1
    marker=attr(size=15, opacity=1, color="transparent",
        line=attr(color="black", width=2)
    )
)
push!(scat_list, scat)

x = range(0,30,length=31)
d = Dist.Binomial(30,0.1)
scat = scatter(x=x,y=Dist.pdf.(d,x), mode="markers",
    xaxis="x4", yaxis="y4",
    #opacity=0.1
    marker=attr(size=10, opacity=1, color="transparent",
        line=attr(color="black", width=2)
    )
)
push!(scat_list, scat)

x = range(0,30,length=31)
d = Dist.Binomial(30,0.5)
scat = scatter(x=x,y=Dist.pdf.(d,x), mode="markers",
    xaxis="x5", yaxis="y5",
    #opacity=0.1
    marker=attr(size=10, opacity=1, color="transparent",
        line=attr(color="black", width=2)
    )
)
push!(scat_list, scat)

x = range(0,30,length=31)
d = Dist.Binomial(30,0.9)
scat = scatter(x=x,y=Dist.pdf.(d,x), mode="markers",
    xaxis="x6", yaxis="y6",
    #opacity=0.1
    marker=attr(size=10, opacity=1, color="transparent",
        line=attr(color="black", width=2)
    )
)
push!(scat_list, scat)

layout = Layout(title="N=3, p=0.1",
    grid=attr(rows=2, columns=3, pattern="independent")
    #yaxis2=attr(range=[0,0.4])
    )
pt = plot(scat_list, layout)
#display(pt)

# Example 1.3, page 15
# The hard way
# P[X=k | theta] = binomial(N,k) * theta^k * (1-theta)^(N-k)
# P[X >= 1] =
Ans = 0.0
theta = 2/9
N = 4
for ind = 1:N
    global Ans
    P = binomial(N,ind) * theta^ind * (1-theta)^(N-ind)
    Ans += P
end
println("Example 1.3, Ans = $Ans")

# The easy way
dist = Dist.Binomial(4,2/9)
Ans = 1 - Dist.cdf(dist,0)
println("Example 1.3, Ans = $Ans")

# The R way
using RCall
R"x <- 1 - dbinom(0,4,2/9)"
Ans = @rget x
println("Example 1.3, the R-way, Ans=$Ans")

# 1.3.2 The Poisson Distribution, page 17
x = range(0,100,length=101)
d = Dist.Poisson(60)
scat = scatter(x=x,y=Dist.pdf.(d,x), mode="lines")
cdf = scatter(x=x,y=Dist.cdf.(d,x), mode="lines", yaxis="y2")
layout = Layout(title="lambda=3",
    yaxis=attr(title="PDF"),
    yaxis2=attr(side="right", overlaying="y")
)
pt = plot([scat,cdf], layout)
display(pt)

# Page 19
# Need to run from terminal to pop R-plot
R"lam <- seq(0, 10, length=50)"
R"y <- dpois(3, lam)"
R"plot(lam, y, xlab=\"lambda\", ylab=\"P[x=3]\", type=\"l\")"

end
nothing #suppress printing of Main...


text:
# Plot recipe plotting an expression as a function
# expression may have 1 (ex1) or 2 (ex2) free variables
# plot of ex:𝑹 → 𝐑: `plot(ex1, a, b)`
# parametric plot: `plot(ex1, ex1′, a, b)`
# surface plot: `surface(xs, ys, ex2)`
# contourplot: `contour(xs, ys, ex2)`.
# not plot([ex1, ex1′], a, b); use plot(ex1, a, b); plot!(ex1′)
@recipe f(::Type{T}, v::T) where {T<:Num} = build_function(v, Symbolics.get_variables(v)...; expression=Val{false})


text:
# Capistran M., Moreles M., Lara B.
# Parameter Estimation of Some Epidemic Models. The Case of Recurrent Epidemics Caused by Respiratory Syncytial Virus
# doi.org/10.1007/s11538-009-9429-3
# Equations (7)-(11)
using Logging

using StructuralIdentifiability

logger = Logging.SimpleLogger(stdout, Logging.Info)
global_logger(logger)

ode = @ODEmodel(
  s'(t) = mu - mu * s(t) - b0 * (1 + b1 * x1(t)) * i(t) * s(t) + g * r(t),
  i'(t) = b0 * (1 + b1 * x1(t)) * i(t) * s(t) - (nu + mu) * i(t),
  r'(t) = nu * i(t) - (mu + g) * r(t),
  x1'(t) = -M * x2(t),
  x2'(t) = M * x1(t),
  y1(t) = i(t),
  y2(t) = r(t)
)

@time println(assess_identifiability(ode))


text:
module Flux

# Zero Flux Given

using Base: tail
using Statistics, Random, LinearAlgebra
using Zygote, MacroTools, Juno, Reexport
using MacroTools: @forward
@reexport using NNlib
using Zygote: Params, @adjoint, gradient, pullback, @nograd

export gradient

export Chain, Dense, Maxout, RNN, LSTM, GRU, SamePad, Conv, CrossCor, ConvTranspose,
       GlobalMaxPool, GlobalMeanPool, MaxPool, MeanPool, flatten,
       DepthwiseConv, Dropout, AlphaDropout, LayerNorm, BatchNorm, InstanceNorm, GroupNorm,
       SkipConnection, params, fmap, cpu, gpu, f32, f64, testmode!, trainmode!

include("optimise/Optimise.jl")
using .Optimise
using .Optimise: @epochs
export Descent, ADAM, Momentum, Nesterov, RMSProp,
  ADAGrad, AdaMax, ADADelta, AMSGrad, NADAM,
  ADAMW, RADAM, InvDecay, ExpDecay, WeightDecay,
  ClipValue, ClipNorm

using CuArrays
const use_cuda = Ref(false)

include("utils.jl")
include("zeros.jl")
include("onehot.jl")
include("functor.jl")

include("layers/stateless.jl")
include("layers/basic.jl")
include("layers/conv.jl")
include("layers/recurrent.jl")
include("layers/normalise.jl")

include("data/Data.jl")

include("deprecations.jl")

include("cuda/cuda.jl")

function __init__()
  use_cuda[] = CuArrays.functional() # Can be overridden after load with `Flux.use_cuda[] = false`
  if CuArrays.functional()
    if !CuArrays.has_cudnn()
      @warn "CuArrays.jl found cuda, but did not find libcudnn. Some functionality will not be available."
    end
  end
end

end # module