text:
using QuantEcon: tauchen, MarkovChain, simulate

# ------------------------------------------------------------------- #
# Define the main Arellano Economy type
# ------------------------------------------------------------------- #

"""
Arellano 2008 deals with a small open economy whose government
invests in foreign assets in order to smooth the consumption of
domestic households. Domestic households receive a stochastic
path of income.

##### Fields
* `β::Real`: Time discounting parameter
* `γ::Real`: Risk aversion parameter
* `r::Real`: World interest rate
* `ρ::Real`: Autoregressive coefficient on income process
* `η::Real`: Standard deviation of noise in income process
* `θ::Real`: Probability of re-entering the world financial sector after default
* `ny::Int`: Number of points to use in approximation of income process
* `nB::Int`: Number of points to use in approximation of asset holdings
* `ygrid::Vector{Float64}`: This is the grid used to approximate income process
* `ydefgrid::Vector{Float64}`: When in default get less income than process
  would otherwise dictate
* `Bgrid::Vector{Float64}`: This is grid used to approximate choices of asset
  holdings
* `Π::Array{Float64, 2}`: Transition probabilities between income levels
* `vf::Array{Float64, 2}`: Place to hold value function
* `vd::Array{Float64, 2}`: Place to hold value function when in default
* `vc::Array{Float64, 2}`: Place to hold value function when choosing to
  continue
* `policy::Array{Float64, 2}`: Place to hold asset policy function
* `q::Array{Float64, 2}`: Place to hold prices at different pairs of (y, B')
* `defprob::Array{Float64, 2}`: Place to hold the default probabilities for
  pairs of (y, B')
"""
immutable ArellanoEconomy
    # Model Parameters
    β::Float64
    γ::Float64
    r::Float64
    ρ::Float64
    η::Float64
    θ::Float64

    # Grid Parameters
    ny::Int
    nB::Int
    ygrid::Array{Float64, 1}
    ydefgrid::Array{Float64, 1}
    Bgrid::Array{Float64, 1}
    Π::Array{Float64, 2}

    # Value function
    vf::Array{Float64, 2}
    vd::Array{Float64, 2}
    vc::Array{Float64, 2}
    policy::Array{Float64, 2}
    q::Array{Float64, 2}
    defprob::Array{Float64, 2}
end

"""
This is the default constructor for building an economy as presented
in Arellano 2008.

##### Arguments
* `;β::Real(0.953)`: Time discounting parameter
* `;γ::Real(2.0)`: Risk aversion parameter
* `;r::Real(0.017)`: World interest rate
* `;ρ::Real(0.945)`: Autoregressive coefficient on income process
* `;η::Real(0.025)`: Standard deviation of noise in income process
* `;θ::Real(0.282)`: Probability of re-entering the world financial sector
  after default
* `;ny::Int(21)`: Number of points to use in approximation of income process
* `;nB::Int(251)`: Number of points to use in approximation of asset holdings
"""
function ArellanoEconomy(;β=.953, γ=2., r=0.017, ρ=0.945, η=0.025, θ=0.282,
                          ny=21, nB=251)

    # Create grids
    Bgrid = collect(linspace(-.4, .4, nB))
    mc = tauchen(ny, ρ, η)
    Π = mc.p
    ygrid = exp(mc.state_values)
    ydefgrid = min(.969 * mean(ygrid), ygrid)

    # Define value functions (Notice ordered different than Python to take
    # advantage of column major layout of Julia)
    vf = zeros(nB, ny)
    vd = zeros(1, ny)
    vc = zeros(nB, ny)
    policy = Array(Int, nB, ny)
    q = ones(nB, ny) .* (1 / (1 + r))
    defprob = Array(Float64, nB, ny)

    return ArellanoEconomy(β, γ, r, ρ, η, θ, ny, nB, ygrid, ydefgrid, Bgrid, Π,
                            vf, vd, vc, policy, q, defprob)
end

u(ae::ArellanoEconomy, c) = c^(1 - ae.γ) / (1 - ae.γ)
_unpack(ae::ArellanoEconomy) =
    ae.β, ae.γ, ae.r, ae.ρ, ae.η, ae.θ, ae.ny, ae.nB
_unpackgrids(ae::ArellanoEconomy) =
    ae.ygrid, ae.ydefgrid, ae.Bgrid, ae.Π, ae.vf, ae.vd, ae.vc, ae.policy, ae.q, ae.defprob

# ------------------------------------------------------------------- #
# Write the value function iteration
# ------------------------------------------------------------------- #
"""
This function performs the one step update of the value function for the
Arellano model-- Using current value functions and their expected value,
it updates the value function at every state by solving for the optimal
choice of savings

##### Arguments

* `ae::ArellanoEconomy`: This is the economy we would like to update the
  value functions for
* `EV::Matrix{Float64}`: Expected value function at each state
* `EVd::Matrix{Float64}`: Expected value function of default at each state
* `EVc::Matrix{Float64}`: Expected value function of continuing at each state

##### Notes

* This function updates value functions and policy functions in place.
"""
function one_step_update!(ae::ArellanoEconomy, EV::Matrix{Float64},
                          EVd::Matrix{Float64}, EVc::Matrix{Float64})

    # Unpack stuff
    β, γ, r, ρ, η, θ, ny, nB = _unpack(ae)
    ygrid, ydefgrid, Bgrid, Π, vf, vd, vc, policy, q, defprob = _unpackgrids(ae)
    zero_ind = searchsortedfirst(Bgrid, 0.)

    for iy=1:ny
        y = ae.ygrid[iy]
        ydef = ae.ydefgrid[iy]

        # Value of being in default with income y
        defval = u(ae, ydef) + β*(θ*EVc[zero_ind, iy] + (1-θ)*EVd[1, iy])
        ae.vd[1, iy] = defval

        for ib=1:nB
            B = ae.Bgrid[ib]

            current_max = -1e14
            pol_ind = 0
            for ib_next=1:nB
                c = max(y - ae.q[ib_next, iy]*Bgrid[ib_next] + B, 1e-14)
                m = u(ae, c) + β * EV[ib_next, iy]

                if m > current_max
                    current_max = m
                    pol_ind = ib_next
                end

            end

            # Update value and policy functions
            ae.vc[ib, iy] = current_max
            ae.policy[ib, iy] = pol_ind
            ae.vf[ib, iy] = defval > current_max ? defval: current_max
        end
    end

    Void
end

"""
This function takes the Arellano economy and its value functions and
policy functions and then updates the prices for each (y, B') pair

##### Arguments

* `ae::ArellanoEconomy`: This is the economy we would like to update the
  prices for

##### Notes

* This function updates the prices and default probabilities in place
"""
function compute_prices!(ae::ArellanoEconomy)
    # Unpack parameters
    β, γ, r, ρ, η, θ, ny, nB = _unpack(ae)

    # Create default values with a matching size
    vd_compat = repmat(ae.vd, nB)
    default_states = vd_compat .> ae.vc

    # Update default probabilities and prices
    copy!(ae.defprob, default_states * ae.Π')
    copy!(ae.q, (1 - ae.defprob) / (1 + r))

    Void
end

"""
This performs value function iteration and stores all of the data inside
the ArellanoEconomy type.

##### Arguments

* `ae::ArellanoEconomy`: This is the economy we would like to solve
* `;tol::Float64(1e-8)`: Level of tolerance we would like to achieve
* `;maxit::Int(10000)`: Maximum number of iterations

##### Notes

* This updates all value functions, policy functions, and prices in place.

"""
function vfi!(ae::ArellanoEconomy; tol=1e-8, maxit=10000)

    # Unpack stuff
    β, γ, r, ρ, η, θ, ny, nB = _unpack(ae)
    ygrid, ydefgrid, Bgrid, Π, vf, vd, vc, policy, q, defprob = _unpackgrids(ae)
    Πt = Π'

    # Iteration stuff
    it = 0
    dist = 10.

    # Allocate memory for update
    V_upd = zeros(ae.vf)

    while dist > tol && it < maxit
        it += 1

        # Compute expectations for this iterations
        # (We need Π' because of order value function dimensions)
        copy!(V_upd, ae.vf)
        EV = ae.vf * Πt
        EVd = ae.vd * Πt
        EVc = ae.vc * Πt

        # Update Value Function
        one_step_update!(ae, EV, EVd, EVc)

        # Update prices
        compute_prices!(ae)

        dist = maxabs(V_upd - ae.vf)

        if it%25 == 0
            println("Finished iteration $(it) with dist of $(dist)")
        end
    end

    Void
end

"""
This function simulates the Arellano economy

##### Arguments

* `ae::ArellanoEconomy`: This is the economy we would like to solve
* `capT::Int`: Number of periods to simulate
* `;y_init::Float64(mean(ae.ygrid)`: The level of income we would like to
  start with
* `;B_init::Float64(mean(ae.Bgrid)`: The level of asset holdings we would like
  to start with

##### Returns

* `B_sim_val::Vector{Float64}`: Simulated values of assets
* `y_sim_val::Vector{Float64}`: Simulated values of income
* `q_sim_val::Vector{Float64}`: Simulated values of prices
* `default_status::Vector{Float64}`: Simulated default status
  (true if in default)

##### Notes

* This updates all value functions, policy functions, and prices in place.

"""
function QuantEcon.simulate(ae::ArellanoEconomy, capT::Int=5000;
                            y_init=mean(ae.ygrid), B_init=mean(ae.Bgrid))

    # Get initial indices
    zero_index = searchsortedfirst(ae.Bgrid, 0.)
    y_init_ind = searchsortedfirst(ae.ygrid, y_init)
    B_init_ind = searchsortedfirst(ae.Bgrid, B_init)

    # Create a QE MarkovChain
    mc = MarkovChain(ae.Π)
    y_sim_indices = simulate(mc, capT+1, y_init_ind)

    # Allocate and Fill output
    y_sim_val = Array(Float64, capT+1)
    B_sim_val, q_sim_val = similar(y_sim_val), similar(y_sim_val)
    B_sim_indices = Array(Int, capT+1)
    default_status = fill(false, capT+1)
    B_sim_indices[1], default_status[1] = B_init_ind, false
    y_sim_val[1], B_sim_val[1] = ae.ygrid[y_init_ind], ae.Bgrid[B_init_ind]

    for t=1:capT
        # Get today's indexes
        yi, Bi = y_sim_indices[t], B_sim_indices[t]
        defstat = default_status[t]

        # If you are not in default
        if !defstat
            default_today = ae.vc[Bi, yi] < ae.vd[yi] ? true: false

            if default_today
                # Default values
                default_status[t] = true
                default_status[t+1] = true
                y_sim_val[t] = ae.ydefgrid[y_sim_indices[t]]
                B_sim_indices[t+1] = zero_index
                B_sim_val[t+1] = 0.
                q_sim_val[t] = ae.q[zero_index, y_sim_indices[t]]
            else
                default_status[t] = false
                y_sim_val[t] = ae.ygrid[y_sim_indices[t]]
                B_sim_indices[t+1] = ae.policy[Bi, yi]
                B_sim_val[t+1] = ae.Bgrid[B_sim_indices[t+1]]
                q_sim_val[t] = ae.q[B_sim_indices[t+1], y_sim_indices[t]]
            end

        # If you are in default
        else
            B_sim_indices[t+1] = zero_index
            B_sim_val[t+1] = 0.
            y_sim_val[t] = ae.ydefgrid[y_sim_indices[t]]
            q_sim_val[t] = ae.q[zero_index, y_sim_indices[t]]

            # With probability θ exit default status
            if rand() < ae.θ
                default_status[t+1] = false
            else
                default_status[t+1] = true
            end
        end
    end

    return (y_sim_val[1:capT], B_sim_val[1:capT], q_sim_val[1:capT],
            default_status[1:capT])
end

max_stars_repo_name: gsulis/structural

max_issues_repo_name: gsulis/structural

max_forks_repo_name: gsulis/structural


text:
function TLEPropagatorBuilder(arg0::TLE, arg1::PositionAngle, arg2::jdouble)
    return TLEPropagatorBuilder((TLE, PositionAngle, jdouble), arg0, arg1, arg2)
end

function TLEPropagatorBuilder(arg0::TLE, arg1::PositionAngle, arg2::jdouble, arg3::DataContext)
    return TLEPropagatorBuilder((TLE, PositionAngle, jdouble, DataContext), arg0, arg1, arg2, arg3)
end

function build_propagator(obj::TLEPropagatorBuilder, arg0::Vector{jdouble})
    return jcall(obj, "buildPropagator", Propagator, (Vector{jdouble},), arg0)
end

ext: jl


text:
# ---
# title: 97. Interleaving String
# id: problem97
# author: Indigo
# date: 2021-06-02
# difficulty: Hard
# categories: String, Dynamic Programming
# link: <https://leetcode.com/problems/interleaving-string/description/>
# hidden: true
# ---
#
# Given strings `s1`, `s2`, and `s3`, find whether `s3` is formed by an
# **interleaving** of `s1` and `s2`.
#
# An **interleaving** of two strings `s` and `t` is a configuration where they
# are divided into **non-empty** substrings such that:
#
#   * `s = s1 + s2 + ... + sn`
#   * `t = t1 + t2 + ... + tm`
#   * `|n - m| <= 1`
#   * The **interleaving** is `s1 + t1 + s2 + t2 + s3 + t3 + ...` or `t1 + s1 + t2 + s2 + t3 + s3 + ...`
#
# **Note:** `a + b` is the concatenation of strings `a` and `b`.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/02/interleave.jpg)
#
#
#
#     Input: s1 = "aabcc", s2 = "dbbca", s3 = "aadbbcbcac"
#     Output: true
#
#
# **Example 2:**
#
#
#
#     Input: s1 = "aabcc", s2 = "dbbca", s3 = "aadbbbaccc"
#     Output: false
#
#
# **Example 3:**
#
#
#
#     Input: s1 = "", s2 = "", s3 = ""
#     Output: true
#
#
#
#
# **Constraints:**
#
#   * `0 <= s1.length, s2.length <= 100`
#   * `0 <= s3.length <= 200`
#   * `s1`, `s2`, and `s3` consist of lower-case English letters.
#
#
## @lc code=start
using LeetCode

function is_interleave(s1::String, s2::String, s3::String)
    len1, len2, len3 = length(s1), length(s2), length(s3)
    len1 + len2 != len3 && return false
    dp = OffsetArray(fill(false, len1 + 1, len2 + 1), -1, -1)
    dp[0, 0] = true
    for i in 0:len1, j in 0:len2
        p = i + j
        i > 0 && (dp[i, j] |= dp[i - 1, j] && s1[i] == s3[p])
        j > 0 && (dp[i, j] |= dp[i, j - 1] && s2[j] == s3[p])
    end
    return dp[end, end]
end
## @lc code=end


text:
# Light properties
source_distance = 1m
source_angle    = 17°

# Eye morphology
aperture        = 251μm
morphing_factor = 1.12

## Volumes
### Refractive indices

ri.water                         = 1.334
ri.cornea                        = 1.37
ri.lens                          = 1.42
ri.lens2distal_retina            = 1.35
ri.distal_retina2proximal_retina = 1.35
ri.gap                           = 1.34

### Thicknesses

thick.cornea                        = 23μm
thick.lens                          = 215μm
thick.lens2distal_retina            = 6μm
thick.distal_retina2proximal_retina = 81μm
thick.gap                           = 49μm

## Interfaces
### Radii of the ellipsoids that form the interface

radii.cornea_distal   = (244μm, 227μm)
radii.lens_distal     = (195μm, 213μm)
radii.lens_proximal   = 337μm
radii.distal_retina   = 337μm
radii.proximal_retina = 337μm
radii.mirror          = 417μm

# Figure ranges
min_aperture = 100μm
max_aperture = 300μm
min_distance = 1cm
max_distance = 1m


text:
"""
    BaseConfig(
        serverurl::String,
        num_comm_rounds::Int,
        fraction_clients::Float32,
        num_total_clients::Int
    )

Part of the configuration which are common to any type of configuration.
"""
struct BaseConfig
    # endpoints
    serverurl::String
    register_node::String
    fit_node::String

    # federated learning
    num_comm_rounds::Int
    fraction_clients::Float32
    num_total_clients::Int

    BaseConfig(
        serverurl::String,
        num_comm_rounds::Int,
        fraction_clients::Float32,
        num_total_clients::Int,
    ) = new(
        serverurl,
        "/register",
        "/fit",
        num_comm_rounds,
        fraction_clients,
        num_total_clients
    )
end


text:
# # Example 8
# This example is from the book Decomposition Techniques in Mathematical Programming
# Chapter 7.2, page 281, [url](https://www.springer.com/gp/book/9783540276852)

# Bold points in example: Here, we are using upperonly, loweronly for defining variables. Also, we are using equality in constraints. Using upperonly/loweronly will create variables that will appear in the proper model. By this definition, we might have slightly smaller reformulated models in the algorithm.

# Model of the problem
# First level
# ```math
# \min -x + 4y + z,\\
# \notag s.t.\\
# y + 2x + z \leq 9,\\
# z = 1,\\
# ```
# Second level
# ```math
# \min -x - y + w,\\
# \notag s.t.\\
# y \geq 0,\\
# x + y + w \leq 8,\\
# x \geq 0,\\
# x \leq 4,\\
# x = 1.\\
# ```

using BilevelJuMP
using Ipopt
using JuMP
using Test

model = BilevelModel(Ipopt.Optimizer, mode = BilevelJuMP.ProductMode(1e-9))

# First we need to create all of the variables in the upper and lower problems:

# Upper level variables
@variable(Upper(model), x)
@variable(UpperOnly(model), z)

# Lower level variables
@variable(Lower(model), y)
@variable(LowerOnly(model), w)

# Then we can add the objective and constraints of the upper problem:

# Upper level objecive function
@objective(Upper(model), Min, -x + 4y + z)

# Upper level constraints
@constraint(Upper(model), y + 2x + z <= 9)
@constraint(Upper(model), z == 1)

# Followed by the objective and constraints of the lower problem:

# Lower objective function
@objective(Lower(model), Min, -x - y + w)

# Lower constraints
@constraint(Lower(model),  y >= 0)
@constraint(Lower(model), x + y + w <= 8)
@constraint(Lower(model),  x >= 0)
@constraint(Lower(model),  x <= 4)
@constraint(Lower(model),  w == 1)

# Initial Starting conditions  #src

# Now we can solve the problem and verify the solution again that reported by book

optimize!(model)

# Auto testing
@test value(x) ≈ 1 atol=1e-2
@test value(y) ≈ 6 atol=1e-2
@test value(z) ≈ 1 atol=1e-2
@test value(w) ≈ 1 atol=1e-2


text:
const ci_cache = GPUCompiler.CodeCache()

const __device_properties_lock = ReentrantLock()
const __device_properties = @NamedTuple{cap::VersionNumber, ptx::VersionNumber, exitable::Bool, debuginfo::Bool}[]
function device_properties(dev)
    @lock __device_properties_lock  begin
        if isempty(__device_properties)
            resize!(__device_properties, ndevices())

            # determine compilation properties of each device
            for dev in devices()
                cap = supported_capability(capability(dev))
                ptx = v"6.3"    # we only need 6.2, but NVPTX doesn't support that

                exitable = true
                if cap < v"7"
                    # JuliaGPU/CUDAnative.jl#4
                    # ptxas for old compute capabilities has a bug where it messes up the
                    # synchronization stack in the presence of shared memory and thread-divergent exit.
                    exitable = false
                end
                if !has_nvml() || NVML.driver_version() < v"460"
                    # JuliaGPU/CUDA.jl#431
                    # TODO: tighten this conditional
                    exitable = false
                end

                debuginfo = false

                __device_properties[deviceid(dev)+1] = (; cap, ptx, exitable, debuginfo)
            end
        end
        @inbounds __device_properties[deviceid(dev)+1]
    end
end

function CUDACompilerTarget(dev::CuDevice; kwargs...)
    PTXCompilerTarget(; device_properties(dev)..., kwargs...)
end

struct CUDACompilerParams <: AbstractCompilerParams end

CUDACompilerJob = CompilerJob{PTXCompilerTarget,CUDACompilerParams}

GPUCompiler.runtime_module(@nospecialize(job::CUDACompilerJob)) = CUDA

# filter out functions from libdevice and cudadevrt
GPUCompiler.isintrinsic(@nospecialize(job::CUDACompilerJob), fn::String) =
    invoke(GPUCompiler.isintrinsic,
           Tuple{CompilerJob{PTXCompilerTarget}, typeof(fn)},
           job, fn) ||
    fn == "__nvvm_reflect" || startswith(fn, "cuda")

function GPUCompiler.link_libraries!(@nospecialize(job::CUDACompilerJob), mod::LLVM.Module,
                                     undefined_fns::Vector{String})
    invoke(GPUCompiler.link_libraries!,
           Tuple{CompilerJob{PTXCompilerTarget}, typeof(mod), typeof(undefined_fns)},
           job, mod, undefined_fns)
    link_libdevice!(mod, job.target.cap, undefined_fns)
end

GPUCompiler.ci_cache(@nospecialize(job::CUDACompilerJob)) = ci_cache

GPUCompiler.method_table(@nospecialize(job::CUDACompilerJob)) = method_table


text:
# %% Demonstration of pitfalls/invalid sets
using DrWatson
@quickactivate :FractalDimension # uses DynamicalSystems, PyPlot
include(srcdir("style.jl"))

datas = Symbol[]
labels = String[]

N = Int(1e4)
dt = 1.0

push!(datas, :lorenz96_nonstationary)
push!(labels, "non-stat. Lorenz96")

push!(datas, :roessler_sphere)
push!(labels, "Periodic+3Dsphere")

# push!(datas, :lorenz96_rounded)
# push!(labels, "rounded Lorenz96")

push!(datas, :experimental_data)
push!(labels, "vostok")

push!(datas, :experimental_data)
push!(labels, "nifty50")

qH = 1
qC = 2
Cmethod = "standard" # bueno or standard. Decides εmax for correlation sum.

# Calculate values for H, C
eHs, eCs, Hs, Cs = [Vector{Float64}[] for i in 1:4]
for i in 1:length(datas)
    # These are the parameters that change between the different elements of the plot
    data = datas[i]

    # Here we simply pack all parameters into a dictionary
    # (other parameters are (probably) globals)
    if string(datas[i])[1:3] != "exp"
        params = @strdict N qH qC data dt
    else
        name = labels[i]
        params = @strdict qH qC data name
    end
    if Cmethod ≠ "standard"
        params["Cmethod"] = Cmethod
    end
    params["z"] = 1

    # This is the main call that calculates everything
    output, s = produce_or_load(
        datadir("main"), params, make_C_H;
        prefix = string(data), suffix = "jld2", force = i == 3,
        ignores = ["data"], storepatch = false,
    )
    @unpack eH, eC, H, C = output
    push!(eHs, eH); push!(Hs, H); push!(eCs, eC); push!(Cs, C)
end

# Do the actual plot
legendtitle = "inappropriate data"

fig, axs = mainplot(Hs, Cs, eHs, eCs, labels, legendtitle; qH, qC, tol = 0.25)
wsave(plotsdir("paper", "inappropriate"), fig)

max_stars_repo_name: Datseris/FractalDimension

max_issues_repo_name: Datseris/FractalDimension

max_forks_repo_name: Datseris/FractalDimension


text:
struct MSQNode{T,R}
    next::R
    data::T
    # Node{T,R}() = new{T,R}()
    MSQNode{T,R}(next) where {T,R} = new{T,R}(next)
    MSQNode{T,R}(next, data) where {T,R} = new{T,R}(next, data)
end

const MSQNodeRef{T} = Reagents.Ref{Union{Nothing,MSQNode{T}}}
const MSQHeadRef{T} = Reagents.Ref{MSQNode{T}}

MSQNode(next::MSQNodeRef{T}) where {T} = MSQNode{T,typeof(next)}(next)
MSQNode(next::MSQNodeRef{T}, data) where {T} = MSQNode{T,typeof(next)}(next, data)

struct MSQueue{T,R<:MSQHeadRef{T}}
    head::R
    tail::R
    function MSQueue{T}() where {T}
        node = MSQNode(MSQNodeRef{T}(nothing))
        head = MSQHeadRef{T}(node)
        tail = MSQHeadRef{T}(node)
        return new{T,typeof(head)}(head, tail)
    end
end

nodetype(::MSQueue{<:Any,R}) where {N,R<:Reagents.Ref{N}} = N

Base.eltype(::Type{<:MSQueue{T}}) where {T} = T

Base.isempty(q::MSQueue) = q.head[].next[] === nothing

trypoppingfirst(q::MSQueue{T}) where {T} =
    Reagents.Update(q.head) do node, _
        next = node.next[]
        if next === nothing
            (node, nothing)
        else
            next::nodetype(q)
            (next, Some(next.data))
        end
    end

pushing(q::MSQueue{T}) where {T} =
    Reagents.Computed() do x
        node = MSQNode(MSQNodeRef{T}(nothing), x)
        while true
            tail = q.tail[]
            tail::nodetype(q)
            next = tail.next[]
            if next === nothing  # found the tail
                return Reagents.CAS(tail.next, nothing, node) ⨟ Reagents.PostCommit() do _
                    Reagents.try!(Reagents.CAS(q.tail, tail, node))
                end
            else  # need the fixup
                next::nodetype(q)
                Reagents.try!(Reagents.CAS(q.tail, tail, next))
            end
        end
    end


text:
@testset "maha" begin
    rng = MersenneTwister(123456)
    x = 2 * rand(rng)
    v1 = rand(rng, 3)
    v2 = rand(rng, 3)

    P = rand(rng, 3, 3)
    k = MahalanobisKernel(P)

    @test kappa(k, x) == exp(-x)
    @test k(v1, v2) ≈ exp(-sqmahalanobis(v1, v2, P))
    @test kappa(ExponentialKernel(), x) == kappa(k, x)
    @test repr(k) == "Mahalanobis Kernel (size(P) = $(size(P)))"
    # test_ADs(P -> MahalanobisKernel(P), P)
    @test_broken "Nothing passes (problem with Mahalanobis distance in Distances)"
end


text:
function lastweek(df)
    l_w = union(df.date)[end]
    last_df = filter(row -> row[:date] == l_w, df)
end

# const to_keep = [:name,
#     :rank_zscore_delta_arena_place,
#     :rank_zscore_delta_arena_victories,
#     :rank_zscore_delta_arena_defences,
#     :rank_zscore_delta_total_power,
#     :rank_zscore_delta_maxed_heroes,
#     :rank_zscore_clan_tokens,
#     :rank_zscore_raid_points,
#     :rank_zscore_portal_stones,
#     :rank_zscore_total]
const to_keep = [:name,
    :zscore_delta_arena_victories,
    :zscore_delta_total_power,
    :zscore_delta_maxed_heroes,
    :zscore_clan_tokens,
    :zscore_raid_points,
    :zscore_portal_stones,
    :zscore_total]
const new_names = [:name,
    :arena_victories,
    :total_power,
    :maxed_heroes,
    :clan_tokens,
    :raid_points,
    :portal_stones,
    :total]

function adjust_to_show(df)
    to_show = select(df,to_keep)
    rename!(to_show, new_names)
    round_names = filter(n-> n != :name, new_names)
    for n in round_names
        println(n)
        to_show[!,n] = round.(to_show[!,n];digits = 2)
    end
    sort!(to_show, :total, rev = true)
end


text:
using DebuggerFramework
using Test
using REPL
using TerminalRegressionTests

include("languages/brainstack.jl")
include("languages/calc.jl")

function setup_repl(emuterm)
    repl = REPL.LineEditREPL(emuterm, true)
    repl.interface = REPL.setup_interface(repl)
    repl.specialdisplay = REPL.REPLDisplay(repl)
    repl
end

const thisdir = dirname(@__FILE__)

# Test stepping in/out of frames
# Language: BrainStack
proga = """
1: 2 3 0
2: 0
3: 0
"""

asta = parse(BrainStack.BrainStackAST, proga)
@test BrainStack.interpret(asta) == 2

TerminalRegressionTests.automated_test(
                joinpath(thisdir,"brainstack/simple.multiout"),
               ["si\n", "f 2\n", "si\n", "si\n", "\n", "\n"]) do emuterm
    DebuggerFramework.debug(asta, setup_repl(emuterm), emuterm)
end

# Test test entering language specific prompt
# Language: Calc
prog = """
1 + (9 * (7 / 2))
_1 * 2
_1 + _2
_1 / 3
"""
calc_code = Calc.parseall(prog)

TerminalRegressionTests.automated_test(
        joinpath(thisdir,"calc/lsrp.multiout"),
        [
        # Evaluate the first three expressions
        "s\ns\ns\n",
        # Look at the first backref
        "`_1\n\b",
        # Make sure we can type ` without erroring (#2)
        "abc`\b\b\b\b`_1\n",
        # Exit
        "\x4"]) do emuterm
    DebuggerFramework.debug(calc_code, setup_repl(emuterm), emuterm)
end


text:
using ForwardDiff, LinearAlgebra, NLsolve
using ForwardDiff: Dual, Tag, tagtype, value, partials

# η = ηₐ ∇a + ηᵦ ∇β + ηᵧ∇γ
# we want η to be a vector variable
# here we want ηi ui to be the ith component of η
# we want ηi to be a scalar variable
# ui needs to be a Co- or Contra-variant vector

# Covariant components Aᵢ
# Contravariant components Aⁱ

abstract type Variance end
struct Covariant <: Variance end
struct Contravariant <: Variance end

abstract type AbstractCoordinateTransform{V<:Variance, D} end
const Act = AbstractCoordinateTransform # short hand

struct CoordinateTransform{V<:Variance, D, F} <: AbstractCoordinateTransform{V, D}
  f::F
end
(ct::CoordinateTransform)(u) = ct.f(u)
Base.:\(c::Act, x) = invert(c, x)
Base.length(c::CoordinateTransform{<:Variance, D}) where {D} = D

const ToCartesian = CoordinateTransform{Covariant}
const FromCartesian = CoordinateTransform{Contravariant}
ToCartesian(f::F, D::Integer) where {F} = CoordinateTransform{Covariant, D, F}(f)
FromCartesian(f::F, D::Integer) where {F} = CoordinateTransform{Contravariant, D, F}(f)

abstract type AbstractScalar{V<:Variance, D} <: Function end
struct Scalar{V<:Variance, D, F} <: AbstractScalar{V, D}
  f::F
  c::CoordinateTransform{V, D}
end
(s::Scalar)(x) = s.f(x)

J⃗(c::Act{Covariant}, r, ::Type{Covariant}) = ForwardDiff.jacobian(c, r)
J⃗(c::Act{Covariant}, r, ::Type{Contravariant}) = inv(J⃗(c, r, Covariant))
J⃗(c::Act{Contravariant}, r, ::Type{Covariant}) = ForwardDiff.jacobian(c.f, c \ r)
J⃗(c::Act{Contravariant}, r, ::Type{Contravariant}) = inv(J⃗(c, r, Covariant))
J⃗(c::Act{V}, x) where {V<:Variance} = J⃗(c, x, V)

covariantbasis(c::Act, x) = J⃗(c, x, Covariant)
contravariantbasis(c::Act, x) = J⃗(c, x, Contravariant)

gᵢⱼ(c::Act, x) = (j = J⃗(c, x); return j' * j)
gⁱʲ(c::Act, x) = inv(gᵢⱼ(c, x))

J(f, x) = sqrt(det(gᵢⱼ(f, x)))
jacobi(f, x) = J⃗(f, x)
jacobian(f, x) = J(f, x)

valueise(x) = x
valueise(x::T) where {T<:Dual} = valueise(value(x))
changevalue(x, y, s) = x == s ? y : x
function changevalue(x::AbstractVector, y, s)
  for i in eachindex(x)
    x[i] = changevalue(x[i], y[i], s[i])
  end
  return x
end
function changevalue(x::AbstractVector{<:Dual}, y, s)
  for i in eachindex(x)
    x[i] = changevalue(x[i], y[i], s[i])
  end
  return x
end
function changevalue(x::Dual{<:Tag}, y, s) where {T<:AbstractFloat, D}
  return Dual{tagtype(x)}(changevalue(x.value, y, s), partials(x))
end
function changevalue(x::Dual{<:Tag{<:Function,T}, T, D}, y::T, s::T
                     ) where {T<:AbstractFloat, D}
  return x.value == s ? Dual{tagtype(x)}(y, partials(x)) : x
end
changevalue(x, y) = changevalue(deepcopy(x), y, valueise.(x))
"""
Given coordinate transform f(y) -> x, find y for given x.
"""
function invert(f::F, x::AbstractVector{T}, initial_y=ones(size(x))
    ) where {F, T}
  f!(m, y) = (m .= f(y) .- valueise.(x); return nothing)
  y = NLsolve.nlsolve(f!, initial_y, autodiff=:forward, ftol=1e-12).zero
  # What is the proper way? This only works for a single AD pass.
  return changevalue(x, y)
end

invert(f::F) where {F} = x -> invert(f, x)
function invert(f::Act{V,D}) where {V<:Variance, D}
  f⁻¹(x) = invert(f, x)
  return CoordinateTransform{convert(V),D,typeof(f⁻¹)}(f⁻¹)
end

struct VectorField{V<:Variance, D, W<:Variance, F,
                   C<:AbstractCoordinateTransform{<:Variance, D}}
  f::F
  c::C
end

const ContravariantVector = VectorField{Contravariant}
const CovariantVector = VectorField{Covariant}

J⃗(v::VectorField, x, ::Type{V}) where {V<:Variance} = J⃗(v.c, x, V)
dl(c::AbstractCoordinateTransform, x, ::Type{Covariant}) = sqrt.(diag(gⁱʲ(c, x)))
dl(c::AbstractCoordinateTransform, x, ::Type{Contravariant}) = sqrt.(diag(gᵢⱼ(c, x)))
dl(v::VectorField{V}, x) where {V<:Variance} = dl(v.c, x, V)

function VectorField{V}(f::F, c::C, construct_from_unit_basis=false) where {
    V<:Variance, D, W<:Variance, F, C<:AbstractCoordinateTransform{W, D}}
  g(x) = construct_from_unit_basis ? f(x) ./ dl(c, x, V) : f(x)
  return VectorField{V, D, W, typeof(g), C}(g, c)
end
(v::VectorField)(x) = v.f(x)

(v::VectorField)(x, normalise::Bool) = normalise ? v(x) .* dl(v, x) : v(x)

function ContravariantVector(Aⱼ::CovariantVector)
  Aⁱ(x) = gⁱʲ(Aⱼ, x) * Aⱼ(x)
  return ContravariantVector(Aⁱ, Aⱼ.c)
end
function CovariantVector(Aʲ::ContravariantVector)
  Aᵢ(x) = gᵢⱼ(Aʲ, x) * Aʲ(x)
  return CovariantVector(Aᵢ, Aʲ.c)
end
CovariantVector(v::CovariantVector) = v
ContravariantVector(v::ContravariantVector) = v

import Base: convert
convert(::Type{Contravariant}) = Covariant
convert(::Type{Covariant}) = Contravariant
convert(Aⁱ::ContravariantVector) = CovariantVector(Aⁱ)
convert(Aᵢ::CovariantVector) = ContravariantVector(Aᵢ)
convert(v::VectorField{V,D,W,F,C}, c::C) where {V,D,W,F,C} = v

const Cov = Covariant
const Con = Contravariant
changecoords(a::Act{Cov,D}, b::Act{Cov,D}) where D = r -> a \ b(r)
changecoords(a::Act{Con,D}, b::Act{Cov,D}) where D = r -> a(b(r))
changecoords(a::Act{Cov,D}, b::Act{Con,D}) where D = r -> a \ (b \ r)
changecoords(a::Act{Con,D}, b::Act{Con,D}) where D = r -> a(b \ r)
function convert(v::VectorField{V,D,<:Variance}, c::Act{<:Variance,D}
    ) where {V<:Variance, D}
  return VectorField{V}(r -> v(changecoords(v.c, c)(r), true), c, true)
end

length(v::VectorField{<:Variance, D}) where {D} = D

dV(a) = J(a)
dV(a, b) = J(a, b)

function dS(c::Act{<:Variance, 3})
  function f(x)
    g_ij = gᵢⱼ(c, x)
    g(n, m) = sqrt(g_ij[n, n] * g_ij[m, m] - g[n, m]^2)
    return [g(filter(i -> i != k, 1:3)...) for k ∈ 1:3]
  end
  return CovariantVector(f, c, false)
end

for op ∈ (:J, :gⁱʲ, :gᵢⱼ)
  @eval $op(v::VectorField, x) = $op(v.c, x)
end

function ∇(f::T, c::Act, construct_from_unit_basis=false) where {T}
  return ContravariantVector(x->∇(f, x), c, construct_from_unit_basis)
end
∇(f::T, x) where {T} = ForwardDiff.jacobian(f, x)
∂(f::T, x) where {T} = inv(∇(f, x))

∇ᵢⱼ(v::VectorField, x) = ForwardDiff.jacobian(f, x)
∇ⁱʲ(v::VectorField, x) = inv(ForwardDiff.jacobian(f, x))

for op in (:+, :-)
  Vs = (ContravariantVector, CovariantVector)
  for (V1, V2) ∈ (Vs, reverse(Vs))
    @eval function Base.$(op)(a::$V1, b::$V2)
      @assert a.c == b.c
      return $V1(x->$op(a(x), $V1(b(x))), a.c)
    end
  end
end
import LinearAlgebra: dot, cross
dot(Aⁱ::ContravariantVector, Bᵢ::CovariantVector) = x->dot(Aⁱ(x), Bᵢ(x))
dot(Aᵢ::CovariantVector, Bⁱ::ContravariantVector) = dot(Bⁱ, Aᵢ)
dot(A::T, B::T) where {V, T<:VectorField{V}} = dot(A, convert(B))
function cross(Aⁱ::ContravariantVector, Bⁱ::ContravariantVector)
  @assert Aⁱ.c == Bⁱ.c
  return CovariantVector(x -> cross(Aⁱ(x), Bⁱ(x)) ./ J(Aⁱ, x), Aⁱ.c)
end
function cross(Aᵢ::CovariantVector, Bᵢ::CovariantVector)
  @assert Aᵢ.c == Bᵢ.c
  return ContravariantVector(x -> cross(Aᵢ(x), Bᵢ(x)) ./ J(Aᵢ, x), Aᵢ.c)
end
cross(A::T, B::T) where {V, T<:VectorField{V}} = cross(A, convert(B))

"""
grad(ϕ) = uⁱ (∂ᵢ ϕ)
"""
function grad(ϕ::F, c::Act{Covariant}) where {F}
  return CovariantVector(x->ForwardDiff.gradient(ϕ, x), c, false)
end
function grad(ϕ::F, c::Act{Contravariant}) where {F}
  return CovariantVector(x->ForwardDiff.gradient(ϕ, c \ x), c, false)
end
grad(s::Scalar) = grad(s.f, s.c)

import Base.div
"""
div(A⃗) = ∂ᵢ(J * Aⁱ) / J
"""
div(v::ContravariantVector) = Scalar(x->div(v, x), v.c)
function div(v::ContravariantVector, r)
  # v.c is CoordinateTransform{Covariant, D}
  return sum(diag(ForwardDiff.jacobian(y->J(v, y) * v(y), r))) / J(v, r)
end
div(v::CovariantVector) = div(ContravariantVector(v))
div(v::CovariantVector, x) = div(ContravariantVector(v), x)

∇o(a) = div(a)
∇o(a, b) = div(a, b)

#"""
#curl(A⃗) = ∑ₖ e⃗ₖ (∂ᵢAⱼ - ∂ⱼAᵢ) / J
#"""
#function curl(v::CovariantVector{3}, x::T) where {T}
#  #∂ⱼAᵢ = ∇(v, x, Covariant)
#  ∂ⱼAᵢ = ForwardDiff.jacobian(v, x)
#  ϵijk∂ⱼAᵢ = [∂ⱼAᵢ[3,2] - ∂ⱼAᵢ[2,3], ∂ⱼAᵢ[1,3] - ∂ⱼAᵢ[3,1], ∂ⱼAᵢ[2,1] - ∂ⱼAᵢ[1,2]]
#  return ϵijk∂ⱼAᵢ ./ J(v, x)
#end
#
#curl(v::CovariantVector) = ContravariantVector(x->curl(v, x), v.c)
#curl(f::F, c::Act) where {F} = curl(CovariantVector(f, c, true))
#curl(v::ContravariantVector) = curl(CovariantVector(v))
#curl(v::ContravariantVector, x) = curl(CovariantVector(v), x)

"""
curl(A⃗) = ∑ₖ e⃗ₖ (∂ᵢAⱼ - ∂ⱼAᵢ) / J
"""
function curl(v::CovariantVector{3}) where {T}
  function inner(x)
    ∂ⱼAᵢ = ForwardDiff.jacobian(v, x)
    ϵijk∂ⱼAᵢ = [∂ⱼAᵢ[3,2] - ∂ⱼAᵢ[2,3], ∂ⱼAᵢ[1,3] - ∂ⱼAᵢ[3,1], ∂ⱼAᵢ[2,1] - ∂ⱼAᵢ[1,2]]
    return ϵijk∂ⱼAᵢ ./ J(v, x)
  end
  return ContravariantVector(inner, v.c)
end
curl(v::ContravariantVector) = curl(CovariantVector(v))
curl(v::ContravariantVector, x) = curl(CovariantVector(v))(x)

#curl(v::CovariantVector) = ContravariantVector(x->curl(v, x), v.c)
#curl(f::F, c::Act) where {F} = curl(CovariantVector(f, c, true))

∇x(a) = curl(a)
∇x(a, b) = curl(a, b)

vectorlaplacian(v::CovariantVector) = grad(div(v)) - curl(curl(v))
vectorlaplacian(v::CovariantVector, x) = grad(div(v))(x) - curl(curl(v))(x)

max_stars_repo_name: jwscook/DifferentialGeometries


text:
@doc """
acf{T<:Number}(x::AbstractArray{T,1}, n::Int=1)

Compute the autocorrelation function of a univariate time series
""" ->
function corlag{T<:Number}(x::AbstractArray{T,1}, n::Int=1)
    if n == 0
        return 1.0
    end
    @assert n > 0
    @assert n < size(x,1)
    return cor(x[1:end-n], x[(n+1):end])
end
function acf{T<:Number}(x::AbstractArray{T,1}, maxlag::Int=15; lags::AbstractArray{Int,1}=0:maxlag)
    @assert all(lags .< size(x,1)-1)
    return map((n) -> corlag(x, n), lags)
end

acf(x::TS, maxlag::Int=15; lags::AbstractArray{Int,1}=0:maxlag) = acf(x.values, maxlag; lags=lags)


text:
using Test
using BoundedStreams

# setup test file
file = joinpath(mkpath(abspath(@__FILE__, "..", "data")), "test.dat")

@testset "$(nameof(S)) general" for S in (BoundedInputStream, BoundedOutputStream)
    io = open(file, "w+")
    is = S(io, 100)
    @test isreadable(is) == (S <: BoundedInputStream)
    @test iswritable(is) == (S <: BoundedOutputStream)
    @test skip(is, 5) === is
    @test position(is) == 5
    @test !ismarked(is)
    @test mark(is) == 5
    @test ismarked(is)
    @test unmark(is)
    mark(is)
    @test seek(is, 1) === is
    @test position(is) == 1
    @test reset(is) == 5
    @test position(is) == 5
    @test bytesavailable(is) == 95
    close(is)
    @test bytesavailable(is) == 0
    @test position(io) == 100
    close(io)
end

@testset "BoundedInputStream read" begin
    open(file, "w") do io
        write(io, '0':'9', 'A':'Z', 'a':'z')
    end
    open(file) do io
        bio = BoundedInputStream(io, 10, close=BoundedStreams.CLOSE)
        @test !eof(bio)
        @test String(read(bio, 4)) == "0123"
        @test position(bio) == 4
        @test position(io) == 4
        @test bytesavailable(bio) == 6
        @test String(read(bio, 5)) == "45678"
        @test position(bio) == 9
        @test bytesavailable(bio) == 1
        @test !eof(bio)
        close(bio)
        @test eof(bio)
        @test eof(io)
    end

    open(file) do io
        bio = BoundedInputStream(io, 3, offset=10, close=26)
        @test String(read(bio)) == "ABC"
        @test position(bio) == 3
        @test position(io) == 13
        @test bytesavailable(bio) == 0
        @test String(read(bio)) == ""
        close(bio)
        @test position(bio) == 26
        @test position(io) == 36
        @test !eof(io)
        bio = BoundedInputStream(io, 5)
        @test String(read(bio)) == "abcde"
        close(bio)
        @test !eof(io)
    end

    open(file) do io
        off = 36
        bio = BoundedInputStream(io, 26, offset=off)
        seek(bio, 10)
        @test position(bio) == 10
        @test position(io) == off + 10
        @test String(read(bio, 2)) == "kl"
        seekstart(bio)
        @test position(bio) == 0
        @test position(io) == off
        seekend(bio)
        @test position(bio) == 26
        @test position(io) == off + 26
        seek(io, 10)
        @test position(bio) == 10 - off
        @test eof(bio)
        @test String(read(bio)) == ""
    end

    open(file) do io
        skip(io, 36)
        bio = BoundedInputStream(io, 40, offset=-26, close=50)
        @test position(io) == 10
        @test skip(bio, 5) === bio
        @test position(io) == 15
        @test seek(bio, 2) === bio
        @test_throws ArgumentError skip(bio, 100)
        @test_throws ArgumentError seek(bio, -200)
        seek(io, 62)
        @test close(bio) === nothing
        @test position(bio) == 50
    end

    open(file) do io
        bio = BoundedInputStream(io, 15)
        @test read(bio, Int64) != nothing
        @test_throws EOFError read(bio, Int64)
    end
end

@testset "BoundedOutputStream write" begin
    open(file, "w") do io
        bio = BoundedOutputStream(io, 20, offset=5, close=BoundedStreams.CLOSE)
        @test write(bio, UInt8(20)) == 1
        @test position(bio) == 1
        @test write(bio, Int64(123)) == 8
        @test position(bio) == 9
        @test bytesavailable(bio) == 11
        @test write(bio, "abcdefghij") == 10
        @test bytesavailable(bio) == 1
        @test_throws EOFError write(bio, "kl")
    end
    open(file, "r") do io
        @test skip(io, 5) === io
        @test read(io, UInt8) == 20
        @test read(io, Int64) == 123
        @test read(io, String) == "abcdefghij"
    end
end


text:
# import Particles package
using Particles

# parameters
Δt = 0.1    # time step size
nt = 100    # number of time steps
ns = 5      # number of samples

# random initial conditions
x₀ = [ [rand() .* 2π, rand() .* 2 .- 1] for i in 1:ns ]

# vector field
function pendulum!(ẋ, x)
    ẋ[1] = x[2]
    ẋ[2] = - sin(x[1])
end

# create an Equation instance
equ = Equation(pendulum!, x₀)

# create Simulation instance
sim = Simulation(equ, ExplicitEuler, Δt, nt)

# run simulation
x = sim()

# load Plots package
# using Plots

# set plot ranges
# xlim = (-π, +π)
# ylim = (-3, +3)

# plot energy contour lines
# plt = contour(
#     LinRange(xlim..., 100),
#     LinRange(ylim..., 100),
#     (x,v) -> v^2 / 2  + 1 - cos(x),
#     xlim = xlim,
#     ylim = ylim,
#     title = "Pendulum",
#     legend = false,
#     size = (800, 600)
# )

# plot solutions
# for i in axes(x,2)
#     scatter!(plt, (x[1,i,:] .+ π) .% 2π .- π, x[2,i,:], marker=5)
# end

# save figure to file
# savefig("pendulum.png")

max_stars_repo_name: michakraus/JuliaAtHomeAssignments

max_issues_repo_name: michakraus/JuliaAtHomeAssignments

max_forks_repo_name: michakraus/JuliaAtHomeAssignments


text:
#= Glacier mass balance toy model

Test with ideal data of a hybrid MB model based on neural networks
for the ablation and accumulation components that respect a classic
temperature-index model, mixed with model interpretation using
SINDy (Brunton et al., 2016).

=#

## Environment and packages
cd(@__DIR__)
using Pkg; Pkg.activate("../."); Pkg.instantiate()
using Infiltrator
using Debugger

using Flux, DiffEqFlux
using Flux: @epochs
using Zygote
using Plots
gr()
using Base: @kwdef
using Statistics
using DataDrivenDiffEq
using ModelingToolkit
using LinearAlgebra
# Set a random seed for reproduceable behaviour
using Random

#Random.seed!(2345)

@kwdef mutable struct Hyperparameters
    batchsize::Int = 500    # batch size
    η::Float64 = 0.1   # learning rate
    epochs::Int = 500        # number of epochs
    use_cuda::Bool = true   # use gpu (if cuda available)
end

###############################################
############  FUNCTIONS   #####################
###############################################

# Fake quadratic response to snowfall and temperature
α(P) = P.^(2)
μ(T) = T.^(2)

# Temperature-index equation for point glacier mass balance
function MB(P, T, Up, Ut)
    T_melt = 0
    MB = Up(P) - Ut(max.(T.-T_melt, 0))
    return MB
end

# Toy function to create artificial reference data
function toy_MB(α, μ, P, T)
    T_melt = 0.0
    MB = α(P) - μ(max.(T.-T_melt, 0.0))
    return MB
end

function train()

    # We determine the loss function
    function loss(batch)
        l, l_acc, l_abl = 0.0f0, 0.0f0, 0.0f0
        num = 0
        for (x, y) in batch

            # Make NN predictions
            p_batch = x[1,:]'
            t_batch = x[2,:]'
            pdd_batch = max.(t_batch.-0, 0)
            Ŷp = Up(p_batch)
            Ŷt = Ut(pdd_batch)

            # We evaluate the MB as the combination of Accumulation - Ablation
            w_pc=1000
            l_MB = sqrt(Flux.Losses.mse(MB(p_batch, t_batch, Up, Ut), y; agg=mean))
            l_range_acc = sum((max.((Ŷp/p_batch).-110, 0)).*w_pc)
            l_range_abl = sum((max.((Ŷt/pdd_batch).-110, 0)).*w_pc)

            #l += l_MB
            l += l_MB + l_range_acc + l_range_abl
            l_acc += l_range_acc
            l_abl += l_range_abl
            num +=  size(x, 2)

            # println("Accumulation loss: ", l_range_acc)
            # println("Ablation loss: ", l_range_abl)

        end

        return l/num
    end

    # Container to track the losses
    losses = Float32[]

    # Callback to show the loss during training
    callback(l) = begin
        push!(losses, l)
        if length(losses)%50==0
            println("Current loss after $(length(losses)) iterations: $(losses[end])")
        end
        false
    end

    function hybrid_train!(loss, ps_Up, ps_Ut, data, opt)
        # Retrieve model parameters
        ps_Up = Params(ps_Up)
        ps_Ut = Params(ps_Ut)

        for batch in data
        # back is a method that computes the product of the gradient so far with its argument.
        train_loss_Up, back_Up = Zygote.pullback(() -> loss(batch), ps_Up)
        train_loss_Ut, back_Ut = Zygote.pullback(() -> loss(batch), ps_Ut)
        # Callback to track the training
        callback(train_loss_Up)
        # Apply back() to the correct type of 1.0 to get the gradient of loss.
        gs_Up = back_Up(one(train_loss_Up))
        gs_Ut = back_Ut(one(train_loss_Ut))
        # Insert what ever code you want here that needs gradient.
        # E.g. logging with TensorBoardLogger.jl as histogram so you can see if it is becoming huge.
        Flux.update!(opt, ps_Up, gs_Up)
        Flux.update!(opt, ps_Ut, gs_Ut)
        # Here you might like to check validation set accuracy, and break out to do early stopping.
        end
    end

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

    ######### Define the network  ############
    # We determine the hyperameters for the training
    hyparams = Hyperparameters()

    # Leaky ReLu as activation function
    leakyrelu(x, a=0.01) = max(a*x, x)
    relu_acc(x) = min(max(0, x), 30)
    relu_abl(x) = min(max(0, x), 35)

    # Define the networks 1->5->5->5->1
    global Up = Chain(
        Dense(1,10,initb = Flux.zeros),
        BatchNorm(10, leakyrelu),
        Dense(10,10,initb = Flux.zeros),
        BatchNorm(10, leakyrelu),
        Dense(10,5,initb = Flux.zeros),
        BatchNorm(5, leakyrelu),
        Dense(5,1, relu, initb = Flux.zeros)
    )

    global Ut = Chain(
        Dense(1,10,initb = Flux.zeros),
        BatchNorm(10, leakyrelu),
        Dense(10,10,initb = Flux.zeros),
        BatchNorm(10, leakyrelu),
        Dense(10,5,initb = Flux.zeros),
        BatchNorm(5, leakyrelu),
        Dense(5,1, relu, initb = Flux.zeros)
    )

    # We define an optimizer
    #opt = RMSProp(hyparams.η, 0.95)
    opt = ADAM(hyparams.η)

    # We get the model parameters to be trained
    ps_Up = Flux.params(Up)
    ps_Ut = Flux.params(Ut)

    sqnorm(x) = sum(abs2, x)

    #######  We generate toy data to train the model  ########
    global snowfall_toy = rand(0.0f0:100.0f0, 500)
    global temperature_toy = 100f0*sin.((1:500)./20)
    forcings = hcat(snowfall_toy, temperature_toy)

    MB_toy = toy_MB(α, μ, snowfall_toy, temperature_toy)

    # Plot the toy dataset
    l1 = @layout [a b; c; d]
    p1₁ = plot(1:500, snowfall_toy, label="Snowfall", color="midnightblue")
    p1₂ = plot(1:500, temperature_toy, label="Temperature", color="darkred")
    hline!(p1₂, [0], c="black", label="")
    p1₃ = plot(1:500, MB_toy, label="Mass balance")
    hline!(p1₃, [0], c="black", label="")

    X = vcat(snowfall_toy', temperature_toy')
    Y = collect(MB_toy')

    # Very important:
    # Shuffle batches and use batchsize similar to dataset size for easy training
    data = Flux.Data.DataLoader((X, Y), batchsize=hyparams.batchsize, (X, Y), shuffle=true)

    # We train the mass balance hybrid model
    @epochs hyparams.epochs hybrid_train!(loss, ps_Up, ps_Ut, data, opt)

    # Perform MB simulations with trained model
    global T_melt = 0
    global snow_norm = X[1,:]
    global temp_norm = X[2,:]
    global PDD_toy = max.(temp_norm.-T_melt, 0)
    global Ŷ_MB = MB(snow_norm', temp_norm', Up, Ut)
    global Ŷ_acc = Up(snow_norm')
    global Ŷ_abl = Ut(PDD_toy')

    # Plot comparison with toy data
    p1₄ = plot(1:500, Ŷ_MB', label="Simulated mass balance")
    hline!(p1₄, [0], c="black", label="")
    p1 = plot(p1₁,p1₂,p1₃, p1₄,layout=l1)

    l2 = @layout [a b;c d]
    p2₁ = plot(1:500, α(snowfall_toy')', label="α")
    plot!(1:500, Ŷ_acc', label="Up")
    hline!(p2₁, [0], c="black", label="")
    p2₂ = plot(1:500, α(snowfall_toy')'./snowfall_toy, label="α/P")
    plot!(1:500, Ŷ_acc'./snowfall_toy, label="Up/P")
    hline!(p2₂, [0], c="black", label="")

    p2₃ = plot(1:500, μ(PDD_toy')', label="μ", lw=2)
    plot!(1:500, Ŷ_abl', label="Ut", lw=2)
    hline!(p2₃, [0], c="black", label="")
    p2₄ = plot(1:500, μ(PDD_toy')'./PDD_toy, label="μ/T", lw=2)
    plot!(1:500, Ŷ_abl'./PDD_toy, label="Ut/T", lw=2)
    hline!(p2₄, [0], c="black", label="")

    p2 = plot(p2₁, p2₂, p2₃, p2₄, layout=l2)

    if(!ispath(pwd(),"plots"))
        mkdir("plots")
    end

    display(p1)
    display(p2)
    savefig(p1,joinpath(pwd(),"plots","hybrid_MB_model.png"))
    savefig(p2,joinpath(pwd(),"plots","hybrid_acc_abl.png"))

    ### Let's take a look at the raw learnt dynamics
    X_raw = 1:150
    Ŷp_raw = Up(X_raw')
    Yp_raw = α(X_raw')
    Ŷt_raw = Ut(X_raw')
    Yt_raw = μ(X_raw')

    l3 = @layout [a b]
    p3₁ = plot(X_raw, Ŷp_raw', lw=3, xlabel="Snowfall (P)", ylabel="Accumulation (α(P))",label="Up (inferred function)")
    plot!(X_raw, Yp_raw', l2=3, lw=3, label="α (true function")
    p3₂ = plot(X_raw, Ŷt_raw', lw=3, xlabel="Temperature (T)", ylabel="Ablation (μ(T))",label="Ut (inferred function)")
    plot!(X_raw, Yt_raw', l2=3, lw=3, label="μ (true function")
    p3 = plot(p3₁, p3₂, layout=l3)
    display(p3)
    savefig(p3,joinpath(pwd(),"plots","recovered_functions.png"))
end

###################################################################
########################  SINDy    ################################
###################################################################

function sindy()

    ## Symbolic regression via sparse regression ( SINDy based )
    # Create a Basis
    Xp = Float64.(copy(reshape(snowfall_toy, (1,length((snowfall_toy))))))
    Xt = max.(temperature_toy.-T_melt, 0)
    @variables X_P[1:1]

    # Generate the basis functions, multivariate polynomials up to deg 5
    b = [polynomial_basis(X_P, 5);]
    basis_p = Basis(b, X_P)

    # Create an optimizer for the SINDy problem
    λ, ρ = 9999e-1, 5.0
    opt = ADMM(λ, ρ)

    #opt = STRRidge(0.1)# Create the thresholds which should be used in the search process

    # opt = SR3(9e-1, 1.0)
    # λs = exp10.(-5:0.1:-1)

    # Test on ideal data
    println("Accumulation:")
    println("SINDy on ideal data \n")
    Ψ_p_true = SINDy(Xp, α(snow_norm'), basis_p, opt, maxiter = 1000) # Succeed
    #Ψ_true = SINDy(Float32.(reshape(X_raw, (1, length(X_raw)))), Float32.(Yp_raw), basis_p, opt, maxiter = 100, normalize = true)
    println(Ψ_p_true)
    print_equations(Ψ_p_true, show_parameter=true)
    p̂_p_true = parameters(Ψ_p_true)
    println("Ψ_true Parameter guess : $(p̂_p_true) \n")

    println("SINDy on synthetic range predicted data \n")
    syn_range = 1:100
    Ψ_p_range = SINDy(Float64.(syn_range), Float64.(Up(syn_range')), basis_p, opt, p=[0], maxiter = 20000, denoise=true, normalize=true)
    println(Ψ_p_range)
    print_equations(Ψ_p_range)
    p̂_p_range = parameters(Ψ_p_range)
    println("Ψ First parameter guess : $(p̂_p_range)")

    println("SINDy on predicted data \n")
    Ψ_p = SINDy(Xp, Float64.(Ŷ_acc), basis_p, opt, p=[0], maxiter = 200000, denoise=false, normalize=true, convergence_error = 1e-8)
    println(Ψ_p)
    print_equations(Ψ_p)
    p̂_p = parameters(Ψ_p)
    println("Ψ First parameter guess : $(p̂_p) \n \n")

    println("Ablation: ")
    println("SINDy on ideal data \n")
    Ψ_t_true = SINDy(Xt, μ(temp_norm'), basis_p, opt, maxiter = 1000) # Succeed
    #Ψ_true = SINDy(Float32.(reshape(X_raw, (1, length(X_raw)))), Float32.(Yp_raw), basis_p, opt, maxiter = 100, normalize = true)
    println(Ψ_t_true)
    print_equations(Ψ_t_true, show_parameter=true)
    p̂_t_true = parameters(Ψ_t_true)
    println("Ψ_true Parameter guess : $(p̂_t_true) \n")

    println("SINDy on synthetic range predicted data \n")
    syn_range = 1:100
    Ψ_t_range = SINDy(Float64.(syn_range), Float64.(Ut(syn_range')), basis_p, opt, p=[0], maxiter = 20000, denoise=true, normalize=true)
    println(Ψ_t_range)
    print_equations(Ψ_t_range)
    p̂_t_range = parameters(Ψ_t_range)
    println("Ψ First parameter guess : $(p̂_t_range)")

    println("SINDy on predicted data \n")
    Ψ_t = SINDy(Xt, Float64.(Ŷ_abl), basis_p, opt, p=[0], maxiter = 200000, denoise=false, normalize=true, convergence_error = 1e-8)
    println(Ψ_t)
    print_equations(Ψ_t)
    p̂_t = parameters(Ψ_t)
    println("Ψ First parameter guess : $(p̂_t) \n")
end

###############################################################
###########################  MAIN #############################
###############################################################

train()

#sindy()


text:
module Planting

# package code goes here

include("utils.jl")
include("geometry.jl")

abstract AbstractObject

immutable SimpleObject <: AbstractObject
    radius::Float64
end

radius(o::SimpleObject) = o.radius

function randomneighbor{T <: AbstractPoint}(p::T, o::AbstractObject)
    r² = radius(o)^2
    R = 2 * r²
    θ = 2π*rand()
    r = √(rand() * R + r²)
    T(x(p) + (r * cos(θ)), y(p) + (r * sin(θ)))
end

function nearbyCells(accelGrid::AccelGrid, radius::Real, p::AbstractPoint)
    xs = (x(p) - 2*radius):accelGrid.gridSize:(x(p) + 2*radius)
    ys = (y(p) - 2*radius):accelGrid.gridSize:(y(p) + 2*radius)
    return (xs, ys)
end

function nearbyItems(p::AbstractPoint, radius::Real, accel::AccelGrid)
    result = eltype(accel)[]
    (xs, ys) = nearbyCells(accel, radius, p)
    # for x in xs, y in ys
    for item in accel[xs, ys]
        if !(item in result)
            push!(result, item)
        end
    end
    return result
end

immutable RejectionPoissonDisk
    N::Int
end

function plant(m::RejectionPoissonDisk, poly::AbstractObject, o::AbstractObject; callback::Dict{String,Function} = Dict{String,Function}(), seed = Dates.datetime2epochms(now()))
    srand(seed)
    (lb, ub) = boundingbox(poly)
    grid = AccelGrid(Int,radius(o)/√(2),(lb...), (ub...))
    final_list = eltype(poly)[]
    done = false
    while !done
        i = 0
        while i < m.N
            i += 1
            np = rand(poly)
            neighborhood = nearbyItems(np,radius(o),grid)
            points = final_list[filter!(x->x!=0,neighborhood)]
        end
    end
end

immutable FastPoissonDisk
    k::Int
end

function plant(m::FastPoissonDisk, poly::AbstractPolygon, o::AbstractObject; callback::Dict{String,Function} = Dict{String,Function}(), seed = Dates.datetime2epochms(now()))
# http://www.cs.ubc.ca/~rbridson/docs/bridson-siggraph07-poissondisk.pdf

    srand(seed)
    (lb, ub) = boundingbox(poly)
    println(radius(o)/√(2))
    grid = AccelGrid(Int,radius(o)/√(2),(lb...), (ub...))

    active_list = eltype(poly)[]
    final_list = eltype(poly)[]

    x₀ = rand(poly)
    push!(active_list,x₀)
    push!(final_list,x₀)
    grid[x₀...] = length(final_list)

    haskey(callback, "init") && callback["init"](poly,o,active_list, final_list, grid)

    genPoint = function (p)
        for i in 1:m.k
            np = randomneighbor(p,o)
            if np in poly
                neighborhood = nearbyItems(np,radius(o),grid)

                points = final_list[filter!(x->x!=0,neighborhood)]

                haskey(callback, "test") && callback["test"](p, np, points,active_list,final_list,grid)

                nearby = any(map(points) do point
                    norm(Line(np,point)) <= radius(o)
                end)

                if !nearby
                    return (true,np)
                end
            end
        end
        return (false,nothing)
    end

    while !isempty(active_list)
        randIdx = rand(1:length(active_list))
        point = active_list[randIdx]
        deleteat!(active_list, randIdx)
        (found,newPoint) = genPoint(point)
        if found
            push!(active_list,point)
            push!(active_list,newPoint)
            push!(final_list,newPoint)
            grid[newPoint...] = length(final_list)
        end
        haskey(callback, "step") && callback["step"](active_list, final_list, grid)
    end

    return final_list
end

export AbstractPoint, Point, AbstractPolygon, Polygon, SimplePolygon, AbstractObject, SimpleObject, x, y, radius, RejectionPoissonDisk, FastPoissonDisk
export plant, Line

# if Pkg.installed("Plots") != nothing
# include("plot.jl")
# end

end # module


text:
male_first_name() = executor(data["en"]["faker"]["name"]["male_first_name"])
female_first_name() = executor(data["en"]["faker"]["name"]["female_first_name"])
first_name() = executor(data["en"]["faker"]["name"]["first_name"])
last_name() = executor(data["en"]["faker"]["name"]["last_name"])
prefixe() = executor(data["en"]["faker"]["name"]["prefix"])
job() = executor(data["en"]["faker"]["name"]["title"]["job"])

function seleccion_field(data::Dict, fields::Tuple)
	e = Dict()
	for (k, v) in data
		if (k in fields)
			push!(e, k => v)
		end
	end
	return e
end

function simple_profile(fields...="None")
	fields[1] == "None" && (fields = []);
	data = Dict(
		"username"=> user_name(),
		"name"=> first_name(),
		"sex"=> rand(["M","F"]),
		"address"=> address(),
		"mail"=> free_email(),
		"birthdate"=> date() )

	if length(fields) > 0
		return  seleccion_field(data, fields)
	else
		return data
	end
end

function profile( fields...="None")
	fields[1] == "None" && (fields = []);
	d = Dict(
		"username"=> user_name(),
		"name"=> first_name(),
		"sex"=> rand(["M","F"]),
		"address"=> address(),
		"mail"=> free_email(),
		"birthdate"=> date(),
		"job"=> job(),
		"company"=> company(),
		"residence"=> address(),
		"current_location"=> (latitude(), longitude()),
		"blood_group"=> rand(["A", "B", "AB", "O"])*rand(["+", "-"]),
		"website"=> url()
		);
	if length(fields) > 0
		return  seleccion_field(d, fields)
	else
		return d
	end
end


text:
using CmdStan

ProjDir = @__DIR__
cd(ProjDir)

shared_file = "shared_funcs.stan"

bernoulli_model = "
  functions{

    #include $(shared_file)
    //#include shared_funcs.stan // a comment
    //#include shared_funcs.stan // a comment
    //#include /Users/rob/shared_funcs.stan // a comment

    void model_specific_function(){
        real x = 1.0;
        return;
    }

  }
  data {
    int<lower=1> N;
    int<lower=0,upper=1> y[N];
  }
  parameters {
    real<lower=0,upper=1> theta;
  }
  model {
    model_specific_function();
    theta ~ beta(my_function(),1);
    y ~ bernoulli(theta);
  }
";

observeddata = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

tmpdir = ProjDir*"/tmp"

stanmodel = Stanmodel(name="bernoulli", model=bernoulli_model,
  printsummary=true, tmpdir=tmpdir);

rc, chn, cnames = stan(stanmodel, observeddata, ProjDir);

if rc == 0
  # Describe the results
  show(chn)
  println()

  # Ceate a ChainDataFrame
  summary_df = read_summary(stanmodel)
  summary_df[:theta, [:mean, :ess]]
end


text:
module Sample

export line_sampler, adaptive_line_sampler, sample_multiple_adaptive_lines,
       InputAdaptiveSampler

using LinearAlgebra
using ProgressMeter
using Rotations
using ..Atoms
using ..Calculators
using ..Clusters
using ..UnitConversions

"""
    _line_sampler(cl1::Cluster, cl2::Cluster, u, step, npoints)

Local helper function.

Returns clusters tuple of cluster arrays that hold positions where clusters
are on different distances from each other in a line like fashion.

# Arguments
- `cl1::Cluster` : first cluster and it's initial location
- `cl2::Cluster` : second cluster and it's initial location
- `u`  : (unit) vector descriping  direction on which `cl2` is moved
- `step` : distance between points
- `npoints` : nunber of points sampled
"""
function _line_sampler(cl1::Cluster, cl2::Cluster, u, step, npoints)
    c1 = deepcopy(cl1)
    c2 = deepcopy(cl2)
    out1 = Array{typeof(c1)}(undef, npoints)
    out2 = Array{typeof(c2)}(undef, npoints)
    for i in eachindex(out1)
        out1[i] = deepcopy(c1)
        out2[i] = deepcopy(c2)
        move!(c2, step*u)
    end
    return out1, out2
end

"""
    line_sampler(cl1::Cluster, cl2::Cluster; npoints=10, mindis=2.0, maxdis=9.0)

Samples cluster position in a line like fashion with even intervals.

Returns clusters tuple of cluster arrays that hold positions where clusters
are on different distances from each other in a line like fashion.

Sampling is done by first generating random direction and placing `cl2` in that
direction from `cl2` at `mindis` and then picking `npoints` in even distances
to `maxdis`.
"""
function line_sampler(cl1::Cluster, cl2::Cluster; npoints=10, mindis=2.0, maxdis=9.0)
    c1 = deepcopy(cl1)
    c2 = deepcopy(cl2)

    center_cluster!(c1)
    center_cluster!(c2)

    rotate_randomly!(c1)
    rotate_randomly!(c2)

    step = (maxdis - mindis) / npoints
    u = [1.0, 0.0, 0.0]
    move!(c2, mindis * u )
    while minimum(distances(c1,c2)) < mindis
        move!(c2, step * u)
    end
     return _line_sampler(c1,c2,u,step,npoints)
end

"""
  adaptive_line_sampler(cal::Calculator, cl1::Cluster, cl2::Cluster, max_e=0; unit="cm-1",
               npoints=10, maxdis=9.0, sstep=0.1, startdistance=3.0, pchannel=undef)

Calculates potential on a line type distances that is suitable for visualization.
Closest distance used in calculations is looked by searching closest distance that
has energy lower than `max_e`

# Arguments
- `max_e`           : Point that is closest and has less energy than this will be starting point for a line

# Keywords
- `unit`            : Unit in which `max_e` is given
- `npoints`         : Number of points in potential
- `maxdis`          : Maximum distance in potential calculation
- `sstep`           : Search step size for adaptive search, also defines accuracy of search
- `startdistance`   : Distance from which adaptive seach is started
- `basename="base"` : prefix for temporary files in calculations
- `id=""`           : extra identificaltion to help to do parallel computations
- `pchannel=undef`  : (Remote)Channel where progess information is added
"""
function adaptive_line_sampler(cal::Calculator, cl1::Cluster, cl2::Cluster, max_e=0;
                               unit="cm-1", npoints=10, maxdis=9.0, sstep=0.1, startdistance=3.0,
                               basename="base", id="", pchannel=undef)
    @debug "Starting adaptive_line_sampler"

    # take deepcopys so that originals are not moved
    c1 = deepcopy(cl1)
    c2 = deepcopy(cl2)

    center_cluster!(c1)
    center_cluster!(c2)

    # Rotate so that we get random orientation
    rotate_randomly!(c1)
    rotate_randomly!(c2)

    # We move the clusters in x-direction
    u = [1.0, 0.0, 0.0]

    # Find good initial point for the search
    # TODO make this more robust
    r = 0.5*startdistance
    move!(c2, 0.5*startdistance*u)
    while minimum(distances(c1,c2)) < startdistance
        move!(c2, sstep * u)
        r += sstep
    end

    # Search for the intial point
    # We need to find a point where potential energy is emax
    # To do this we find the most distant point where energy is more than emax
    # and the closest point where energy is less than e
    emax = energy_from(max_e, unit)
    Δe = bsse_corrected_energy(cal, [c1], [c2], basename=basename, id=id, pchannel=pchannel) .- emax
    ctemp = [deepcopy(c2)] # History for points we tried
    fless = false # We have a point where Δe < 0
    fmore = false # We have a point where Δe > 0
    if Δe[1] > 0
        fmore = true
    else
        fless = true
    end
    while fmore == false || fless == false
        if Δe[end] < 0
            move!(c2, -1*sstep*u)
            r -= sstep
        else
            move!(c2, sstep*u)
            r += sstep
        end
        tmp = bsse_corrected_energy(cal, [c1], [c2], basename=basename, id=id, pchannel=pchannel) .- emax
        push!(Δe,  tmp...)
        push!(ctemp, deepcopy(c2))
        if Δe[end] > 0
            fmore = true
        else
            fless = true
        end
    end
    if length(Δe) == 1
        error("Search fail!")
    end

    # Pretare output and perform some checks
    @debug "adaptive_line_sampler found initial point in $(length(e)) steps"
    out2 = [ctemp[Δe.<0][end]] # We take the closest point where e < emax
    eout = [Δe[Δe.<0][end] + emax] # Calculated potential energy
    out1 = [c1]  # cluster1 coordinates
    c2 = deepcopy(out2[1])
    rr = minimum(distances(c1,c2))
    @debug "Minimum distance at start" r rr
    @assert maxdis > rr "You need to increase maximum distance in search"

    # Lets now find the points
    step = (maxdis - rr ) / npoints
    move!(c2,step*u) # We already calculated the first one, so we need to move
    tmp = _line_sampler(c1, c2, u, step, npoints-1)

    # Now we can calculate the potential
    et = bsse_corrected_energy(cal, tmp[1], tmp[2], basename=basename, id=id, pchannel=pchannel)
    push!(eout, et...)
    push!(out1, tmp[1]...)
    push!(out2, tmp[2]...)
    @debug "Minimum distance at end $(minimum(distances(out1[end],out2[end])))"
    return Dict("Energy" => eout, "Points" =>  out1 .+ out2,  "Mindis" =>  minimum(distances(out1[1],out2[1])))
end

"""
    InputAdaptiveSampler

Structure used to help use of adaptive line samplers

# Fields
- `cal::Calculator`   : [`calculator`](@ref) used in calculations
- `cl1::Cluster`      : first cluster
- `cl2::Cluster`      : second cluster
- `nlines`            : number of lines calculated
- `max_e`             : maximum energy for configuration
- `unit=`             : unit for maximum energy ( see [`energy_from`](@ref) for supported units)
- `npoints`           : number of points in line
- `maxdis`            : maximum cluster distance
- `sstep`             : search step used by [`adaptive_line_sampler`](@ref)
- `startdistance`     : distance where search is started
"""
mutable struct InputAdaptiveSampler
    # This is needed for parallelization.
    # pmap needs that all data is send to othe process
    # wich is intended to be done with this struct

    # TODO clean this - it is way too messy now
    cal
    cl1
    cl2
    nlines
    max_e
    unit
    npoints
    maxdis
    sstep
    startdistance
    function InputAdaptiveSampler(cal::Calculator, cl1::Cluster, cl2::Cluster,
                                  nlines, max_e=0;
                                  unit="cm-1", npoints=10, maxdis=9.0, sstep=0.1,
                                  startdistance=3.0)
        new(cal,cl1,cl2,nlines,max_e,unit,npoints,maxdis,sstep,startdistance)
    end
end

"""
sample_multiple_adaptive_lines(cal::Calculator, cl1::Cluster, cl2::Cluster, nlines, max_e=0;
                               unit="cm-1", npoints=10, maxdis=9.0, sstep=0.1, startdistance=3.0,
                               basename="base", id="")

Calculates energy of two given clusters in different distances of each other.
Distances are sampled in line like fashion in random directions with even distances.
Sampling of line is done with [`adaptive_line_sampler`](@ref).

# Arguments
- `cal::Calculator`   : [`calculator`](@ref) used in calculations
- `cl1::Cluster`      : first cluster
- `cl2::Cluster`      : second cluster
- `nlines`            : number of lines calculated
- `max_e=0`           : maximum energy for configuration
- `unit="cm-1"`       : unit for maximum energy ( see [`energy_from`](@ref) for supported units)
- `npoints=10`        : number of points in line
- `maxdis=9.0`        : maximum cluster distance
- `sstep=0.1`         : search step used by [`adaptive_line_sampler`](@ref)
- `startdistance=3.0` : distance where search is started
- `basename="base"` : prefix for temporary files in calculations
- `id=""`           : extra identificaltion to help to do parallel computations
- `pchannel=undef`                : (Remote)Channel where progess information is added

# Returns
[`Dict`](@ref) with keys
* "Energy" : array holding energy
* "Points" : array of [`Cluster`](@ref) representing points where energy was calculated
* "Mindis" : array of minimum distances in each line
"""
function sample_multiple_adaptive_lines(cal::Calculator, cl1::Cluster, cl2::Cluster, nlines, max_e=0;
                               unit="cm-1", npoints=10, maxdis=9.0, sstep=0.1, startdistance=3.0,
                               basename="base", id="", pchannel=undef)
    c1 = deepcopy(cl1)
    c2 = deepcopy(cl2)
    out = Vector{Dict}(undef, nlines)
    rtmp = Float64[]
    sr = startdistance
    for i in 1:nlines
        tmp = adaptive_line_sampler(cal, c1, c2, max_e; unit=unit, npoints=npoints,
                                    maxdis=maxdis, startdistance=sr, basename=basename,
                                     id=id, pchannel=pchannel)
        push!(rtmp, tmp["Mindis"])
        sr = sum(rtmp) / length(rtmp)
        out[i] = tmp
    end
    energy = hcat(map( x -> x["Energy"], out)...)
    points = hcat(map( x -> x["Points"], out)...)
    mindis = map(x -> x["Mindis"], out)
    return Dict("Energy" => energy, "Points" => points, "Mindis" => mindis)
end

"""
sample_multiple_adaptive_lines(inputs::InputAdaptiveSampler; basename="base", id="", pchannel=undef)

Calculates energy of two given clusters in different distances of each other.
Distances are sampled in line like fashion in random directions with even distances.
Sampling of line is done with [`adaptive_line_sampler`](@ref).

# Arguments
- `inputs::InputAdaptiveSampler`  : input information in clean package
- `basename="base"`               : prefix for temporary files in calculations
- `id=""`                         : extra identificaltion to help to do parallel computations
- `pchannel=undef`                : (Remote)Channel where progess information is added

# Returns
[`Dict`](@ref) with keys
* "Energy" : array holding energy
* "Points" : array of [`Cluster`](@ref) representing points where energy was calculated
* "Mindis" : array of minimum distances in each line
"""
function sample_multiple_adaptive_lines(inputs::InputAdaptiveSampler; basename="base", id="", pchannel=undef)
    sample_multiple_adaptive_lines(inputs.cal, inputs.cl1, inputs.cl2,
                 inputs.nlines, inputs.max_e,
                 unit=inputs.unit, npoints=inputs.npoints, maxdis=inputs.maxdis,
                 startdistance=inputs.startdistance, basename=basename, id=id, pchannel=pchannel)
end

end  # module sample


text:
"""
$(SIGNATURES)

Compute the Cusp normal form.

# Arguments
- `F, dF, d2F, d3F`: function `(x, p) -> F(x, p)` and its differentials `(x, p, dx) -> d1F(x, p, dx)`, `(x, p, dx1, dx2) -> d2F(x, p, dx1, dx2)`...
- `pt::Cusp` Cusp bifurcation point
- `ls` linear solver

# Optional arguments
- `δ = 1e-8` used for finite differences
- `verbose` bool to print information
"""
function cuspNormalForm(F, dF, d2F, d3F,
		br::AbstractBranchResult, ind_bif::Int;
		δ = 1e-8,
		nev = length(eigenvalsfrombif(br, ind_bif)),
		Jᵗ = nothing,
		verbose = false,
		ζs = nothing,
		lens = br.lens,
		Teigvec = getvectortype(br),
		scaleζ = norm,
		issymmetric = false)
	@assert getvectortype(br) <: BorderedArray
	@assert br.specialpoint[ind_bif].x isa BorderedArray
	@assert br.specialpoint[ind_bif].type == :cusp "The provided index does not refer to a Cusp Point"

	verbose && println("#"^53*"\n--> Cusp Normal form computation")

	# scalar type
	T = eltype(Teigvec)
	ϵ2 = T(δ)

	# functional
	prob = br.functional

	# linear solvers
	ls = prob.linsolver
	bls = prob.linbdsolver

	# kernel dimension
	N = 1

	# in case nev = 0 (number of unstable eigenvalues), we increase nev to avoid bug
	nev = min(2N, nev)

	# Newton parameters
	options = br.contparams.newtonOptions

	# bifurcation point
	bifpt = br.specialpoint[ind_bif]
	eigRes = br.eig

	# eigenvalue
	if bifpt.ind_ev > 0
		λ = eigRes[bifpt.idx].eigenvals[bifpt.ind_ev]
	else
		λ = rightmost(eigRes[bifpt.idx].eigenvals)[1]
	end

	# parameters for vector field
	p = bifpt.param
	parbif = set(br.params, lens, p)
	parbif = set(parbif, prob.lens, get(bifpt.printsol, prob.lens))

	# jacobian at bifurcation point
	x0 = convert(Teigvec.parameters[1], bifpt.x.u)
	L = dF(x0, parbif)

	# eigenvectors
	# we recompute the eigen-elements if there were not saved during the computation of the branch
	@info "Eigen-elements not saved in the branch. Recomputing them..."
	eigsolver = getsolver(options.eigsolver)
	_λ0, _ev0, _ = eigsolver(L, nev)
	Ivp = sortperm(_λ0, by = abs)
	_λ = _λ0[Ivp]
	if norm(_λ[1:N] .- 0, Inf) > br.contparams.precisionStability
		@warn "We did not find the correct eigenvalues. We found the eigenvalues:\n $(display(( _λ[1:N]))).\n Difference between the eigenvalues:"
		display(_λ[1:N] .- 0)
	end
	ζ = real.(geteigenvector(eigsolver, _ev0, Ivp[1]))
	ζ ./= scaleζ(ζ)

	# extract eigen-elements for adjoint(L), needed to build spectral projector
	if issymmetric
		λstar = br.eig[bifpt.idx].eigenvals[bifpt.ind_ev]
		ζstar = copy(ζ)
	else
		_Jt = isnothing(Jᵗ) ? adjoint(L) : Jᵗ(x0, parbif)
		ζstar, λstar = getAdjointBasis(_Jt, conj(λ), eigsolver; nev = nev, verbose = verbose)
	end

	ζstar = real.(ζstar); λstar = real.(λstar)

	@assert abs(dot(ζ, ζstar)) > 1e-10 "We got ζ⋅ζstar = $((dot(ζ, ζstar))). This dot product should not be zero. Perhaps, you can increase `nev` which is currently $nev."
	ζstar ./= dot(ζ, ζstar)

	# Kuznetsov, Yu. A. “Numerical Normalization Techniques for All Codim 2 Bifurcations of Equilibria in ODE’s.” SIAM Journal on Numerical Analysis 36, no. 4 (January 1, 1999): 1104–24. https://doi.org/10.1137/S0036142998335005.
	# notations from this paper
	B(dx1, dx2) = d2F(x0, parbif, dx1, dx2)
	C(dx1, dx2, dx3) = d3F(x0, parbif, dx1, dx2, dx3)
	q = ζ; p = ζstar

	h2 = B(q, q)
	h2 .= dot(p, h2) .* q .- h2
	H2, = bls(L, q, p, zero(T), h2, zero(T))

	c = dot(p, C(q, q, q)) + 3dot(p, B(q, H2))
	c /= 6

	pt = Cusp(
		x0, parbif,
		(lens, prob.lens),
		ζ, ζstar,
		(c = c, ),
		:none
	)
end

"""
$(SIGNATURES)

Compute the Bogdanov-Takens normal form.

# Arguments
- `F, dF, d2F, d3F`: function `(x, p) -> F(x, p)` and its differentials `(x, p, dx) -> d1F(x, p, dx)`, `(x, p, dx1, dx2) -> d2F(x, p, dx1, dx2)`...
- `pt::BogdanovTakens` BogdanovTakens bifurcation point
- `ls` linear solver

# Optional arguments
- `δ = 1e-8` used for finite differences
- `verbose` bool to print information
- `autodiff = true` only for Bogdanov-Takens point. Whether to use ForwardDiff for the many differentiations that are required to compute the normal form.
- `detailed = true` only for Bogdanov-Takens point. Whether to compute only a simplified normal form.
"""
function bogdanovTakensNormalForm(prob, L, d2F, d3F,
							pt::BogdanovTakens;
							δ = 1e-8,
							verbose = false,
							detailed = true,
							autodiff = true)
	x0 = pt.x0
	parbif = pt.params
	Ty = eltype(x0)

	# for finite differences
	ϵ = convert(Ty, δ)
	ϵ2 = sqrt(ϵ) # this one is for second order differential

	F = prob.F
	J = prob.J

	# linear solvers
	ls = prob.linsolver
	bls = prob.linbdsolver

	lens1, lens2 = pt.lens

	getp(l::Lens) = get(parbif, l)
	setp(l::Lens, p::Number) = set(parbif, l, p)
	setp(p1::Number, p2::Number) = set(set(parbif, lens1, p1), lens2, p2)

	ζ0, ζ1 = pt.ζ
	ζs0, ζs1 = pt.ζstar

	G = [dot(xs, x) for xs in pt.ζstar, x in pt.ζ]
	norm(G-I(2), Inf) > 1e-5 && @warn "G == I(2) is not valid. We built a basis such that G = $G"

	G = [dot(xs, L*x) for xs in pt.ζstar, x in pt.ζ]
	norm(G-[0 1;0 0], Inf) > 1e-5 && @warn "G is not close to the Jordan block of size 2. We built a basis such that G = $G. The norm of the difference is $(norm(G-[0 1;0 0], Inf))"

	# second differential
	R2(dx1, dx2) = d2F(x0, parbif, dx1, dx2) ./2

	# quadratic coefficients
	R20 = R2(ζ0, ζ0)
	a = dot(ζs1, R20)
	b = 2dot(ζs0, R20) + 2dot(ζs1, R2(ζ0, ζ1))

	# return the normal form coefficients
	pt.nf = (; a, b)
	if detailed == false
		return pt
	end

	###########################
	# computation of the unfolding. We follow the procedure described in Al-Hdaibat et al. 2016

	# Al-Hdaibat, B., W. Govaerts, Yu. A. Kuznetsov, and H. G. E. Meijer. “Initialization of Homoclinic Solutions near Bogdanov--Takens Points: Lindstedt--Poincaré Compared with Regular Perturbation Method.” SIAM Journal on Applied Dynamical Systems 15, no. 2 (January 2016): 952–80. https://doi.org/10.1137/15M1017491.
	###########################
	# to have the same notations as in the paper above
	q0 = ζ0 ; q1 = ζ1;
	p0 = ζs0; p1 = ζs1;

	# second differential notations, to be in agreement with Kuznetsov et al.
	B(dx1, dx2) = d2F(x0, parbif, dx1, dx2)
	Ainv(dx) = bls(L, p1, q0, zero(Ty), dx, zero(Ty))

	H2000, = Ainv(2 .* a .* q1 .- B(q0, q0))
	γ = (-2dot(p0, H2000) + 2dot(p0, B(q0, q1)) + dot(p1, B(q1, q1))) / 2
	H2000 .+= γ .* q0

	H1100, = Ainv(b .* q1 .+ H2000 .- B(q0, q1))
	H0200, = Ainv(2 .* H1100 .- B(q1, q1))

	# first order drivatives
	pBq(p, q) = 2 .* (applyJacobian(prob, x0 + ϵ * q, parbif, p, true) .-
					  applyJacobian(prob, x0, parbif, p, true)) ./ ϵ
	A1(q, lens) = (applyJacobian(prob, x0, setp(lens, get(parbif, lens) + ϵ), q) .-
	 				  applyJacobian(prob, x0, parbif, q)) ./ϵ
	pAq(p, q, lens) =  dot(p, A1(q, lens))

	# second order derivative
	p10 = get(parbif, lens1); p20 = get(parbif, lens2);

	if autodiff
		Jp(p, l)  = ForwardDiff.derivative( P -> F(x0, setp(l, P)) , p)
		Jpp(p, l) = ForwardDiff.derivative( P -> Jp(P, l), p)
		Fp(p1, p2)  = F(x0, setp(p1, p2))
		Jp1p2(p1, p2) = ForwardDiff.derivative(P1 -> ForwardDiff.derivative(P2 -> Fp(P1, P2) , p2), p1)

		J2_11 = Jpp(p10, lens1)
		J2_22 = Jpp(p20, lens2)
		J2_12 = Jp1p2(p10, p20)
	else #finite  differences. We need to be carreful here because (1e-8)^2 is really small!!
		J2_11 = (F(x0, setp(lens1, p10 + ϵ2)) .- 2 .* F(x0, setp(lens1, p10)) .+
				    F(x0, setp(lens1, p10 - ϵ2)) ) ./ ϵ2^2

		J2_22 = (F(x0, setp(lens2, p20 + ϵ2)) .- 2 .* F(x0, setp(lens2, p20)) .+
				   F(x0, setp(lens2, p20 - ϵ2)) )./ ϵ2^2

		J2_12 = (F(x0, setp(p10 + ϵ2, p20 + ϵ2)) .- F(x0, setp(lens1, p10 + ϵ2)) .-
													F(x0, setp(lens2, p20 + ϵ2)) .+ F(x0, parbif))./ ϵ2^2
	end

	# build the big matrix of size (n+2) x (n+2) A = [L J1s; A12 A22]
	J1 = lens -> F(x0, setp(lens, get(parbif, lens) + ϵ)) ./ ϵ
	J1s = (J1(lens1), J1(lens2))

	A12_1 = pBq(p1, q0) ./2
	A12_2 = (pBq(p0, q0) .+ pBq(p1, q1)) ./2
	A22 = [[pAq(p1, q0, lens1), pAq(p0, q0, lens1)+pAq(p1, q1, lens1)] [pAq(p1, q0, lens2), pAq(p0, q0, lens2)+pAq(p1, q1, lens2)] ]

	# solving the linear system of size n+2
	# TODO REMOVE THIS HACK FOR MATRIX-FREE
	Anp2 = [L hcat(J1s...); hcat(A12_1, A12_2)' A22]
	c = 3dot(p0, H1100) - dot(p0, B(q1, q1))

	sol = Anp2 \ vcat(q1, dot(p1, B(q1, q1))/2, c)
	H0010 = sol[1:end-2]
	@assert size(H0010) == size(x0)
	K10 = sol[end-1:end]

	sol = Anp2 \ vcat(zero(q1), zero(Ty), one(Ty))
	H0001 = sol[1:end-2]
	@assert size(H0001) == size(x0)
	K11 = sol[end-1:end]

	# computation of K2
	κ1 = dot(p1, B(H0001, H0001))
	κ2 = pAq(p1, H0001, lens1) * K11[1] +
		 pAq(p1, H0001, lens2) * K11[2]
	J2K = @. J2_11 * K11[1]^2 + 2J2_12 * K11[1] * K11[2] + J2_11 * K11[2]^2
	κ3 = dot(p1, J2K)
	K2 = -( κ1 + 2κ2 + κ3 ) .* K10

	# computation of H0002
	h0002 = B(H0001, H0001)
	h0002 .+= A1(H0001, lens1) .* (2K11[1]) .+ A1(H0001, lens2) .* (2K11[2])
	h0002 .+= J2K
	h0002 .+= J1s[1] .* K2[1] .+ J1s[2] .* K2[2]
	H0002, = Ainv(h0002)
	H0002 .*= -1

	# computation of H1001
	h1001 = B(q0, H0001)
	h1001 .+= A1(q0, lens1) .* K11[1] .+ A1(q0, lens2) .* K11[2]
	H1001, = Ainv(h1001)
	H1001 .*= -1

	# computation of H0101
	h0101 = B(q1, H0001)
	h0101 .+= A1(q1, lens1) .* K11[1] .+ A1(q1, lens2) .* K11[2]
	h0101 .-= H1001 .+ q1
	H0101, = Ainv(h0101)
	H0101 .*= -1

	# computation of H3000 and d
	h3000 = d3F(x0, parbif, q0, q0, q0) .+ 3 .* B(q0, H2000) .- (6a) .* H1100
	d = dot(p1, h3000)/6
	h3000 .-= (6d) .* q1
	H3000, = Ainv(h3000)
	H3000 .*= -1

	# computation of e
	e = dot(p1, d3F(x0, parbif, q0, q0, q0)) + 2dot(p1, B(q0, H1100)) + dot(p1, B(q1, H2000))
	e += -2b * dot(p1, H1100) -2a * dot(p1, H0200) - dot(p1, H3000)
	e /= 2

	# computation of H2001 and a1
	B1(q, p, l) = (d2F(x0, setp(l, getp(l) + ϵ), q, p) .- d2F(x0, parbif, q, p)) ./ ϵ
	h2001 = d3F(x0, parbif, q0, q0, H0001) .+ 2 .* B(q0, H1001) .+ B(H0001, H2000)
	h2001 .+= B1(q0, q0, lens1) .* K11[1] .+ B1(q0, q0, lens2) .* K11[2]
	h2001 .+= A1(H2000, lens1)  .* K11[1] .+ A1(H2000, lens2)  .* K11[2]
	h2001 .-= (2a) .* H0101
	a1 = dot(p1, h2001) / 2
	h2001 .-= (2a1) .* q1
	H2001, = Ainv(h2001)
	H2001 .*= -1

	# computation of b1
	b1 = dot(p1, d3F(x0, parbif, q0, q1, H0001)) +
		 dot(p1, B1(q0, q1, lens1)) * K11[1] +
		 dot(p1, B1(q0, q1, lens2)) * K11[2] +
		 dot(p1, B(q1, H1001)) +
		 dot(p1, B(H0001, H1100)) +
		 dot(p1, B(q0, H0101)) +
		 dot(p1, A1(H1100, lens1)) * K11[1] + dot(p1, A1(H1100, lens2)) * K11[2] -
		 b * dot(p1, H0101) - dot(p1, H1100) - dot(p1, H2001)

	verbose && println(pt.nf)
	return @set pt.nfsupp = (; γ, c, K10, K11, K2, d, e, a1, b1, H0001, H0010, H0002, H1001, H2000)
end

function predictor(bt::BogdanovTakens, ::Val{:HopfCurve}, ds::T; verbose = false, ampfactor = T(1)) where T
	# If we write the normal form [y2, β1 + β2 y2 + a y1^2 + b y1 y2]
	# equilibria y2 = 0, 0 = β1 + a y1^2
	# Characteristic polynomial: t^2 + (-x*b - β2)*t - 2*x*a
	# the fold curve is β1 / a < 0 with x± := ±√(-β1/a)v
	# the Hopf curve is 0 = -x*b - β2, -x⋅a > 0
	# ie β2 = -bx with ±b√(-β1/a)
	@unpack a, b = bt.nf
	@unpack K10, K11, K2 = bt.nfsupp
	lens1, lens2 = bt.lens
	p1 = get(bt.params, lens1)
	p2 = get(bt.params, lens2)
	par0 = [p1, p2]
	getx(s) = a > 0 ? -sqrt(abs(s) / a) : sqrt(abs(s) / abs(a))
	function HopfCurve(s)
		# x = getx(s)
		if a > 0
			x = -sqrt(abs(s) / a)
			β1 = -abs(s)
		else
			x = sqrt(abs(s) / abs(a))
			β1 = abs(s)
		end
		β2 = -b * x
		ω = sqrt(-2x*a)
		return (pars = par0 .+ K10 .* β1 .+ K11 .* β2 .+ K2 .* (β2^2/2), ω = ω)
	end
	# compute eigenvector corresponding to the Hopf branch
	function EigenVec(s)
		x = getx(s)
		# the jacobian is [0 1; 2x*a b*X+β2] with b*X+β2 = 0
		A = [0 1; 2x*a 0]
		F = eigen(A)
		ind = findall(imag.(F.values) .> 0)
		hopfvec = F.vectors[:, ind]
		return bt.ζ[1] .* hopfvec[1] .+ bt.ζ[2] .* hopfvec[2]
	end
	function EigenVecAd(s)
		x = getx(s)
		# the jacobian is [0 1; 2x*a b*X+β2] with b*X+β2 = 0
		A = [0 1; 2x*a 0]'
		F = eigen(A)
		ind = findall(imag.(F.values) .< 0)
		hopfvec = F.vectors[:, ind]
		return bt.ζstar[1] .* hopfvec[1] .+ bt.ζstar[2] .* hopfvec[2]
	end
	# compute point on the Hopf curve
	x0 = getx(ds)

	return (hopf = t->HopfCurve(t).pars,
			ω = t->HopfCurve(t).ω,
			EigenVec = EigenVec,
			EigenVecAd = EigenVecAd,
			x0 = t -> getx(t) .* bt.ζ[1])
end

function predictor(bt::BogdanovTakens, ::Val{:FoldCurve}, ds::T; verbose = false, ampfactor = T(1)) where T
	# If we write the normal form [y2, β1 + β2 y2 + a y1^2 + b y1 y2]
	# equilibria y2 = 0, 0 = β1 + a y1^2
	# the fold curve is β1 / a < 0 with x± := ±√(-β1/a)
	# the Hopf curve is 0 = -x*b - β2, x⋅a > 0
	# ie β2 = -bx with ±b√(-β1/a)
	@unpack a, b = bt.nf
	@unpack K10, K11, K2 = bt.nfsupp
	lens1, lens2 = bt.lens
	p1 = get(bt.params, lens1)
	p2 = get(bt.params, lens2)
	par0 = [p1, p2]
	getx(s) = a > 0 ? -sqrt(abs(s) / a) : sqrt(abs(s) / abs(a))
	function FoldCurve(s)
		β1 = 0
		β2 = s
		return par0 .+ K10 .* β1 .+ K11 .* β2 .+ K2 .* (β2^2/2)
	end
	return (fold = FoldCurve,
			EigenVec = t -> (bt.ζ[1]),
			EigenVecAd = t -> (bt.ζstar[2]),
			x0 = t -> getx(t) .* bt.ζ[1])
end

function predictor(bt::BogdanovTakens, ::Val{:HomoclinicCurve}, ds::T; verbose = false, ampfactor = T(1)) where T
	# we follow
	# Al-Hdaibat, B., W. Govaerts, Yu. A. Kuznetsov, and H. G. E. Meijer. “Initialization of Homoclinic Solutions near Bogdanov--Takens Points: Lindstedt--Poincaré Compared with Regular Perturbation Method.” SIAM Journal on Applied Dynamical Systems 15, no. 2 (January 2016): 952–80. https://doi.org/10.1137/15M1017491.

	@unpack a, b = bt.nf
	@unpack K10, K11, K2, b1, e, d, a1 = bt.nfsupp
	@unpack H0001, H0010, H0002, H1001, H2000 = bt.nfsupp

	lens1, lens2 = bt.lens
	p1 = get(bt.params, lens1)
	p2 = get(bt.params, lens2)
	par0 = [p1, p2]

	# formula 63
	τ2 = 4/a * (25/49*b1 - e/b) + 2/(49a^2) * (144/49b^2 - 25b*a1 + 73d)

	# formula 69
	α(ϵ) = @. par0 + (10b*ϵ^2 / (7a)) * K11 + ϵ^4/a * ( -4*K10 + 50b^2/(49a) * K2 + b * τ2 * K11)

	# formula 71
	q0, q1 = bt.ζ

	u0(ξ) = -6sech(ξ)^2 + 2
	v0(ξ) = 12sech(ξ)^2 * tanh(ξ)
	u1(ξ) = 0
	v1(ξ) = -6b/(7a) * tanh(ξ) * v0(ξ)
	u2(ξ) = -3/(49a^2) * (6b^2 - 70b*a1 + 49d) * sech(ξ)^2 - 2(5a1*b + 7d)/(7a^2)

	function xLP(t, ϵ)
		ξ = ϵ * t
		return @. bt.x0 + (ϵ^2/a) * ( (10b/7) * H0001 + u0(ξ) * q0) +
			(ϵ^3/a) * ( v0(ξ) * q1 + u1(ξ) * q0) +
			(ϵ^4/a) * ( -4 * H0010 + 50 * b^2/(49a) * H0002 + b*τ2 * H0001 +
						u2(ξ) * q0 + v1(ξ) * q1 +
						1/(2a) * u0(ξ)^2 * H2000 + 10b/(7a) * u0(ξ) * H1001)
	end

	return (α = α, orbit = xLP)
end

"""
$(SIGNATURES)

Compute the Bogdanov-Takens normal form.

# Arguments
- `F, dF, d2F, d3F`: function `(x, p) -> F(x, p)` and its differentials `(x, p, dx) -> d1F(x, p, dx)`, `(x, p, dx1, dx2) -> d2F(x, p, dx1, dx2)`...
- `br` branch result from a call to [`continuation`](@ref)
- `ind_bif` index of the bifurcation point in `br`
- `options` options for the Newton solver

# Optional arguments
- `Jᵗ` is the jacobian adjoint, used for computation of the eigen-elements of the jacobian adjoint, needed to compute the spectral projector
- `δ = 1e-8` used for finite differences
- `nev = 5` number of eigenvalues to compute to estimate the spectral projector
- `verbose` bool to print information
- `autodiff = true` only for Bogdanov-Takens point. Whether to use ForwardDiff for the many differentiations that are required to compute the normal form.
- `detailed = true` only for Bogdanov-Takens point. Whether to compute only a simplified normal form.
"""
function bogdanovTakensNormalForm(F, dF, d2F, d3F,
		br::AbstractBranchResult, ind_bif::Int;
		δ = 1e-8,
		nev = length(eigenvalsfrombif(br, ind_bif)),
		Jᵗ = nothing,
		verbose = false,
		ζs = nothing,
		lens = br.lens,
		Teigvec = getvectortype(br),
		scaleζ = norm,
		detailed = true,
		autodiff = true)
	@assert getvectortype(br) <: BorderedArray
	@assert br.specialpoint[ind_bif].type == :bt "The provided index does not refer to a Bogdanov-Takens Point"

	# functional
	prob = br.functional
	@assert prob isa AbstractProblemMinimallyAugmented

	# bordered linear solver
	bls = prob.linbdsolver

	# kernel dimension:
	N = 2

	# in case nev = 0 (number of requested eigenvalues), we increase nev to avoid bug
	nev = max(2N, nev)
	verbose && println("#"^53*"\n--> Bogdanov-Takens Normal form computation")

	# Newton parameters
	optionsN = br.contparams.newtonOptions

	# bifurcation point
	bifpt = br.specialpoint[ind_bif]
	eigRes = br.eig

	# # eigenvalue
	# if bifpt.ind_ev > 0
	# 	λ = eigRes[bifpt.idx].eigenvals[bifpt.ind_ev]
	# else
	# 	λ = rightmost(eigRes[bifpt.idx].eigenvals)[1]
	# end

	# parameters for vector field
	p = bifpt.param
	parbif = set(br.params, lens, p)
	parbif = set(parbif, prob.lens, get(bifpt.printsol, prob.lens))

	# jacobian at bifurcation point
	@assert bifpt.x isa BorderedArray "Need to fill THIS !!!!!!!!!!!"
	x0 = convert(Teigvec.parameters[1], bifpt.x.u)
	Ty = eltype(Teigvec)
	L = dF(x0, parbif)

	# and corresponding eigenvectors
	eigsolver = getsolver(optionsN.eigsolver)
	if isnothing(ζs) # do we have a basis for the kernel?
		if haseigenvector(br) == false # are the eigenvector saved in the branch?
			@info "No eigenvector recorded, computing them on the fly"
			# we recompute the eigen-elements if there were not saved during the computation of the branch
			_λ0, _ev, _ = eigsolver(L, nev)
			Ivp = sortperm(_λ0, by = abs)
			_λ = _λ0[Ivp]
			verbose && (println("--> (λs, λs (recomputed)) = "); display(( _λ[1:N])))
			if norm(_λ[1:N] .- 0, Inf) > br.contparams.precisionStability
				@warn "We did not find the correct eigenvalues (see 1st col). We found the eigenvalues displayed in the second column:\n $(display(( _λ[1:N]))).\n Difference between the eigenvalues:"
				display(_λ[1:N] .- 0)
			end
			ζs = [copy(geteigenvector(eigsolver, _ev, ii)) for ii in Ivp[1:N]]
		else
			ζs = [copy(geteigenvector(eigsolver, br.eig[bifpt.idx].eigenvec, ii)) for ii in indev-N+1:indev]
		end
	end
	###########################
	# Construction of the basis (ζ0, ζ1), (ζstar0, ζstar1). We follow the procedure described in Al-Hdaibat et al. 2016 on page 972.

	# Al-Hdaibat, B., W. Govaerts, Yu. A. Kuznetsov, and H. G. E. Meijer. “Initialization of Homoclinic Solutions near Bogdanov--Takens Points: Lindstedt--Poincaré Compared with Regular Perturbation Method.” SIAM Journal on Applied Dynamical Systems 15, no. 2 (January 2016): 952–80. https://doi.org/10.1137/15M1017491.
	###########################
	vext = real.(ζs[1])
	_λstar, _evstar, _ = eigsolver(transpose(L), nev)
	Ivp = sortperm(_λstar, by = abs)
	# in case the prob is HopfMA, we real it
	zerov = real.(prob.zero)
	wext = real.(geteigenvector(eigsolver, _evstar, Ivp[1]))
	q0, = bls(L, wext, vext, zero(Ty), zerov, one(Ty))
	p1, = bls(transpose(L), vext, wext, zero(Ty), zerov, one(Ty))
	q1, = bls(L, p1, q0, zero(Ty), q0, zero(Ty))
	p0, = bls(transpose(L), q0, p1, zero(Ty), p1, zero(Ty))
	# we want
	# A⋅q0 = 0, A⋅q1 = q0
	# At⋅p1 = 0, At⋅p0 = p1

	μ = √(abs(dot(q0, q0)))
	q0 ./= μ
	q1 ./= μ
	q1 .= q1 .- dot(q0, q1) .* q0
	ν = dot(q0, p0)
	p1 ./= ν
	p0 .= p0 .- dot(p0, q1) .* p1
	p0 ./= ν

	pt = BogdanovTakens(
		x0, parbif, (prob.lens, lens),
		(q0, q1), (p0, p1),
		(a = zero(Ty), b = zero(Ty) ),
		(K2 = zero(Ty),),
		:none
	)
	return bogdanovTakensNormalForm(prob, L, d2F, d3F, pt; δ = δ, verbose = verbose, detailed = detailed, autodiff = autodiff)
end
####################################################################################################
function bautinNormalForm(F, dF, d2F, d3F,
		br::AbstractBranchResult, ind_bif::Int;
		δ = 1e-8,
		nev = length(eigenvalsfrombif(br, ind_bif)),
		Jᵗ = nothing,
		verbose = false,
		ζs = nothing,
		lens = br.lens,
		Teigvec = getvectortype(br),
		scaleζ = norm)
	@assert getvectortype(br) <: BorderedArray
	@assert br.specialpoint[ind_bif].type == :gh "The provided index does not refer to a Bautin Point"

	verbose && println("#"^53*"\n--> Bautin Normal form computation")

	# scalar type
	T = eltype(Teigvec)
	ϵ2 = T(δ)

	# functional
	prob = br.functional
	@assert prob isa HopfProblemMinimallyAugmented
	ls = prob.linsolver
	bls = prob.linbdsolver

	# ``kernel'' dimension
	N = 2

	# in case nev = 0 (number of unstable eigenvalues), we increase nev to avoid bug
	nev = max(N, nev)

	# Newton parameters
	optionsN = br.contparams.newtonOptions

	# bifurcation point
	bifpt = br.specialpoint[ind_bif]
	eigRes = br.eig

	# eigenvalue
	ω = abs(bifpt.x.p[2])
	λ = Complex(0, ω)

	# parameter for vector field
	p = bifpt.param
	parbif = set(br.params, lens, p)
	parbif = set(parbif, prob.lens, get(bifpt.printsol, prob.lens))

	# jacobian at bifurcation point
	x0 = convert(Teigvec.parameters[1], bifpt.x.u)
	L = dF(x0, parbif)

	# right eigenvector
	if haseigenvector(br) == false
		# we recompute the eigen-elements if there were not saved during the computation of the branch
		@info "Recomputing eigenvector on the fly"
		_λ, _ev, _ = optionsN.eigsolver.eigsolver(L, nev)
		_ind = argmin(abs.(_λ .- λ))
		@info "The eigenvalue is $(_λ[_ind])"
		@assert abs(_λ[_ind] - λ)<br.contparams.newtonOptions.tol "We did not find the correct eigenvalue $λ. We found $(_λ[_ind])"
		ζ = geteigenvector(optionsN.eigsolver, _ev, _ind)
	else
		ζ = copy(geteigenvector(optionsN.eigsolver ,br.eig[bifpt.idx].eigenvec, bifpt.ind_ev))
	end
	ζ ./= scaleζ(ζ)

	# left eigen-elements
	_Jt = isnothing(Jᵗ) ? adjoint(L) : Jᵗ(x, p)
	ζstar, λstar = getAdjointBasis(_Jt, conj(_λ[_ind]), optionsN.eigsolver.eigsolver; nev = nev, verbose = verbose)

	# check that λstar ≈ conj(λ)
	abs(λ + λstar) > 1e-2 && @warn "We did not find the left eigenvalue for the Hopf point to be very close to the imaginary part, $λ ≈ $(λstar) and $(abs(λ + λstar)) ≈ 0?\n You can perhaps increase the number of computed eigenvalues, the number is nev = $nev"

	# normalise left eigenvector
	ζstar ./= dot(ζ, ζstar)
	@assert dot(ζ, ζstar) ≈ 1

	# parameters for vector field
	p = bifpt.param
	parbif = set(br.params, lens, p)
	parbif = set(parbif, prob.lens, get(bifpt.printsol, prob.lens))

	# second order differential, to be in agreement with Kuznetsov et al.
	B = BilinearMap( (dx1, dx2) -> d2F(x0, parbif, dx1, dx2) )
	C = TrilinearMap((dx1, dx2, dx3) -> d3F(x0, parbif, dx1, dx2, dx3) )

	q0 = ζ; p0 = ζstar
	cq0 = conj.(q0)

	H20, = ls(L, B(q0, q0); a₀ = Complex(0, 2ω), a₁ = -1)
	H11, = ls(L, -B(q0, cq0))
	H30, = ls(L, C(q0, q0, q0) .+ 3 .* B(q0, H20); a₀ = Complex(0, 3ω), a₁ = -1)

	h21 = C(q0, q0, cq0) .+ B(cq0, H20) .+ 2 .* B(q0, H11)
	G21 = dot(p0, h21)
	h21 .= G21 .* q0 .- h21
	H21, = bls(L, q0, p0, zero(T), h21, zero(T); shift = Complex{T}(0, -ω))
	# sol = [L-λ*I q0; p0' 0] \ [h21..., 0]

	# 4-th order coefficient
	d4F(x0, dx1, dx2, dx3, dx4) = (d3F(x0 .+ ϵ2 .* dx4, parbif, dx1, dx2, dx3) .-
								   d3F(x0 .- ϵ2 .* dx4, parbif, dx1, dx2, dx3)) ./(2ϵ2)

	# implement 4th order differential with finite differences
	function D(x0, dx1, dx2, dx3, dx4)
		dx4r = real.(dx4); dx4i = imag.(dx4);
		# C(dx, dx4r) + i * C(dx, dx4i)
		trilin_r = TrilinearMap((_dx1, _dx2, _dx3) -> d4F(x0, _dx1, _dx2, _dx3, dx4r) )
		out1 = trilin_r(dx1, dx2, dx3)
		trilin_i = TrilinearMap((_dx1, _dx2, _dx3) -> d4F(x0, _dx1, _dx2, _dx3, dx4i) )
		out2 = trilin_i(dx1, dx2, dx3)
		return out1 .+ im .* out2
	end

	h31 = D(x0, q0, q0, q0, cq0) .+ 3 .* C(q0, q0, H11) .+ 3 .* C(q0, cq0, H20) .+ 3 .* B(H20, H11)
	h31 .+= B(cq0, H30) .+ 3 .* B(q0, H21) .- (3 * G21) .* H20
	H31, = ls(L, h31; a₀ = Complex(0, 2ω), a₁ = -1)

	h22 = D(x0, q0, q0, cq0, cq0) .+
		4 .* C(q0, cq0, H11) .+ C(cq0, cq0, H20) .+ C(q0, q0, conj.(H20)) .+
		2 .* B(H11, H11) .+ 2 .* B(q0, conj.(H21)) .+ 2 .* B(cq0, H21) .+ B(conj.(H20), H20) .-
		(2G21 + 2conj(G21)) .* H11
	H22, = ls(L, h22)
	H22 .*= -1

	# 5-th order coefficient
	# implement 5th order differential with finite differences
	function E(dx1, dx2, dx3, dx4, dx5)
		dx5r = real.(dx5); dx5i = imag.(dx5);
		out1 = (D(x0 .+ ϵ2 .* dx5r, dx1, dx2, dx3, dx4) .-
			    D(x0 .- ϵ2 .* dx5r, dx1, dx2, dx3, dx4)) ./(2ϵ2)
	    out2 = (D(x0 .+ ϵ2 .* dx5i, dx1, dx2, dx3, dx4) .-
			    D(x0 .- ϵ2 .* dx5i, dx1, dx2, dx3, dx4)) ./(2ϵ2)
		return out1 .+ im .* out2
	end

	G32 = dot(p0, E(q0, q0, q0, cq0, cq0))
	G32 += dot(p0, D(x0, q0, q0, q0, conj.(H20))) +
		  3dot(p0, D(x0, q0, cq0, cq0, H20)) +
		  6dot(p0, D(x0, q0, q0, cq0, H11))

	G32 += dot(p0, C(cq0, cq0, H30)) +
		  3dot(p0, C(q0, q0, conj.(H21))) +
		  6dot(p0, C(q0, cq0, H21)) +
		  3dot(p0, C(q0, conj.(H20), H20)) +
		  6dot(p0, C(q0, H11, H11)) +
		  6dot(p0, C(cq0, H20, H11))

	G32 += 2dot(p0, B(cq0, H31)) +
		   3dot(p0, B(q0, H22)) +
		    dot(p0, B(conj(H20), H30)) +
		   3dot(p0, B(conj(H21), H20)) +
		   6dot(p0, B(H11, H21))

	# second Lyapunov coefficient
	l2 = real(G32) / 12

	pt = Bautin(
		x0, parbif,
		(lens, prob.lens),
		ζ, ζstar,
		(;ω, G21, G32, l2 ),
		:none
	)
end
####################################################################################################


text:
using LinearAlgebra

@testset "Circle" begin
    @testset "CirclePointRadius" begin
        img = zeros(Gray{N0f8}, 10, 10)

        expected = Gray{N0f8}[  0 0 1 1 1 1 1 0 0 0
                                0 1 1 1 1 1 1 1 0 0
                                1 1 1 1 1 1 1 1 1 0
                                1 1 1 1 1 1 1 1 1 0
                                1 1 1 1 1 1 1 1 1 0
                                1 1 1 1 1 1 1 1 1 0
                                1 1 1 1 1 1 1 1 1 0
                                0 1 1 1 1 1 1 1 0 0
                                0 0 1 1 1 1 1 0 0 0
                                0 0 0 0 0 0 0 0 0 0 ]

        res = @inferred draw(img, CirclePointRadius(Point(5,5), 5))
        @test all(expected .== res) == true

        res = @inferred draw(img, CirclePointRadius(5, 5, 5))
        @test all(expected .== res) == true

        res = @inferred draw(img, CirclePointRadius(CartesianIndex(5,5), 5))
        @test all(expected .== res) == true

        res = @inferred draw(img, CirclePointRadius(Point(5,5), 5), Gray{N0f8}(0.5))
        expected = map(i->(i==1 ? Gray{N0f8}(0.5) : Gray{N0f8}(0)), expected)
        @test all(expected .== res) == true

        img = zeros(Gray{N0f8}, 10, 10)
        expected = Gray{N0f8}[ 0 0 0 0 0 0 0 0 0 0
                               0 0 0 1 1 1 1 1 0 0
                               0 0 1 1 1 1 1 1 1 0
                               0 1 1 1 0 0 0 1 1 1
                               0 1 1 0 0 0 0 0 1 1
                               0 1 1 0 0 0 0 0 1 1
                               0 1 1 0 0 0 0 0 1 1
                               0 1 1 1 0 0 0 1 1 1
                               0 0 1 1 1 1 1 1 1 0
                               0 0 0 1 1 1 1 1 0 0 ]

        res = @inferred draw(img, CirclePointRadius(Point(6, 6), 5; thickness = 3, fill = false))
        @test all(expected .== res)

        err = ArgumentError("Thickness 7 should be smaller than 5.")
        @test_throws err CirclePointRadius(CartesianIndex(6, 6), 5; thickness = 7, fill = false)

        @test draw(img, CirclePointRadius(Point(6, 6), 5; thickness = 0, fill = true)) == draw(img, CirclePointRadius(Point(6, 6), 5))

        @test draw(img, CirclePointRadius(CartesianIndex(5,5), 5; thickness = 0, fill = true)) == draw(img, CirclePointRadius(CartesianIndex(5,5), 5))

        @test draw(img, CirclePointRadius(5, 5, 5; thickness = 0, fill = true)) == draw(img, CirclePointRadius(5, 5, 5))

        # expected = CirclePointRadius(Point(6, 6), 5; thickness=UInt8(1), fill=false)
        # res = CirclePointRadius(Point(6, 6), 5; thickness=1, fill=false)
        # @test expected == res
    end

    @testset "CircleThreePoints" begin
        img = zeros(Gray{N0f8}, 10, 10)

        expected = Gray{N0f8}[ 0 0 0 0 0 0 0 0 0 0
                               0 0 0 1 1 1 1 1 0 0
                               0 0 1 1 1 1 1 1 1 0
                               0 1 1 1 1 1 1 1 1 1
                               0 1 1 1 1 1 1 1 1 1
                               0 1 1 1 1 1 1 1 1 1
                               0 1 1 1 1 1 1 1 1 1
                               0 1 1 1 1 1 1 1 1 1
                               0 0 1 1 1 1 1 1 1 0
                               0 0 0 1 1 1 1 1 0 0 ]

        res = @inferred draw(img, CircleThreePoints(Point(1,5), Point(5,1), Point(10,5)))
        @test all(expected .== res)

        res = @inferred draw(img, CircleThreePoints(1,5, 5,1, 10,5))
        @test all(expected .== res)

        res = @inferred draw(img, CircleThreePoints(CartesianIndex(5,1), CartesianIndex(1,5), CartesianIndex(5,10)))
        @test all(expected .== res)

        res = @inferred draw(img, CircleThreePoints(Point(1,5), Point(5,1), Point(10,5)), Gray{N0f8}(0.5))
        expected = map(i->(i==1 ? Gray{N0f8}(0.5) : Gray{N0f8}(0)), expected)
        @test all(expected .== res) == true

        invalid_circle = CircleThreePoints(Point(1,1), Point(2,2), Point(3,3))
        @test_throws ErrorException draw(img, invalid_circle)

        invalid_circle = CircleThreePoints(Point(1,1), Point(5,5), Point(10,1))
        @test_throws ErrorException draw(img, invalid_circle)

    end
end