text:
"""
    getARXregressor(y::AbstractVector,u::AbstractVecOrMat, na, nb; inputdelay = ones(Int, size(nb)))
Returns a shortened output signal `y` and a regressor matrix `A` such that the least-squares ARX model estimate of order `na,nb` is `y\\A`
Return a regressor matrix used to fit an ARX model on, e.g., the form
`A(z)y = B(z)f(u)`
with output `y` and input `u` where the order of autoregression is `na`,
the order of input moving average is `nb` and an optional input delay `inputdelay`. Caution, changing the input delay changes the order to `nb + inputdelay - 1`. An `inputdelay = 0` results in a direct term.
# Example
Here we test the model with the Function `f(u) = √(|u|)`
```julia
A     = [1,2*0.7*1,1] # A(z) coeffs
B     = [10,5] # B(z) coeffs
u     = randn(100) # Simulate 100 time steps with Gaussian input
y     = filt(B,A,u)
yr,A  = getARXregressor(y,u,3,2) # We assume that we know the system order 3,2
x     = A\\yr # Estimate model polynomials
plot([yr A*x], lab=["Signal" "Prediction"])
```
For nonlinear ARX-models, see [BasisFunctionExpansions.jl](https://github.com/baggepinnen/BasisFunctionExpansions.jl/). See also `arx`
"""
function getARXregressor(y::AbstractVector, u::AbstractVecOrMat, na, nb; inputdelay = ones(Int, size(nb)))
    length(nb) == size(u, 2) ||
        throw(ArgumentError("Length of nb must equal number of input signals"))
    size(nb) == size(inputdelay) || throw(ArgumentError("inputdelay has to have the same structure as nb"))
    m = max(na, maximum(nb .+ inputdelay .- 1)) + 1 # Start of yr
    @assert m >= 1
    n = length(y) - m + 1 # Final length of yr
    @assert n <= length(y)
    A = toeplitz(y[m:m+n-1], y[m:-1:m-na])
    @assert size(A, 2) == na + 1
    y = A[:, 1] # extract yr
    A = A[:, 2:end]
    for i = 1:length(nb)
        nb[i] <= 0 && continue
        s = m - inputdelay[i]
        A = [A toeplitz(u[s:s+n-1, i], u[s:-1:s-(nb[i])+1, i])]
    end
    return y, A
end
getARXregressor(d::AbstractIdData, na, nb; inputdelay = ones(Int, size(nb))) =
    getARXregressor(time1(output(d)), time1(input(d)), na, nb; inputdelay = inputdelay)

function getARregressor(dy::AbstractIdData, na)
    noutputs(dy) == 1 || throw(ArgumentError("Only 1d time series supported"))
    y = time1(output(dy))
    getARregressor(vec(y), na)
end

function getARregressor(y::AbstractVector, na)
    m = na + 1 # Start of yr
    n = length(y) - m + 1 # Final length of yr
    A = toeplitz(y[m:m+n-1], y[m:-1:1])
    @assert size(A, 2) == na + 1
    y = A[:, 1] # extract yr
    A = A[:, 2:end]
    return y, A
end

@userplot Find_na
"""
    find_na(y::AbstractVector,n::Int)
Plots the RMSE and AIC For model orders up to `n`. Useful for model selection
"""
find_na
@recipe function find_na(p::Find_na)
    y, n = p.args[1:2]
    error = zeros(n, 2)
    for i = 1:n
        yt, A = getARXregressor(y, 0y, i, 0)
        e = yt - A * (A \ yt)
        error[i, 1] = rms(e)
        error[i, 2] = aic(e, i)
    end
    layout --> 2
    title --> ["RMS error" "AIC"]
    seriestype --> :scatter
    @series begin
        error
    end
end

@userplot Find_nanb
"""
    find_nanb(d::InputOutputData,na,nb)
Plots the RMSE and AIC For model orders up to `n`. Useful for model selection
"""
find_nanb
@recipe function find_nanb(p::Find_nanb; logrms = false)
    d, na, nb = p.args[1:3]
    y, u = time1(output(d)), time1(input(d))
    error = zeros(na, nb, 2)
    for i = 1:na, j = 1:nb
        yt, A = getARXregressor(y, u, i, j)
        e = yt - A * (A \ yt)
        error[i, j, 1] = logrms ? log10.(rms(e)) : rms(e)
        error[i, j, 2] = aic(e, i + j)
    end
    layout --> 2
    seriestype --> :heatmap
    xticks := (1:nb, 1:nb)
    yticks := (1:na, 1:na)
    yguide := "na"
    xguide := "nb"
    @series begin
        title := "RMS error"
        subplot := 1
        error[:, :, 1]
    end
    @series begin
        title := "AIC"
        subplot := 2
        error[:, :, 2]
    end
end

"""
    Gtf = arx(d::AbstractIdData, na, nb; inputdelay = ones(Int, size(nb)), λ = 0, estimator=\\, stochastic=false)

Fit a transfer Function to data using an ARX model and equation error minimization.

- `nb` and `na` are the number of coefficients of the numerator and denominator polynomials.
Input delay can be added via `inputdelay = d`, which corresponds to an additional delay of `z^-d`.
An `inputdelay = 0` results in a direct term.
The highest order of the B polynomial is given by `nb + inputdelay - 1`.  `λ > 0` can be provided for L₂ regularization.
`estimator` defaults to \\ (least squares), alternatives are `estimator = tls` for total least-squares estimation.
`arx(Δt,yn,u,na,nb, estimator=wtls_estimator(y,na,nb)` is potentially more robust in the presence of
heavy measurement noise. The number of free parameters is `na+nb`
- `stochastic`: if true, returns a transfer function with uncertain parameters represented by `MonteCarloMeasurements.Particles`.

Supports MISO estimation by supplying an iddata with a matrix `u`, with nb = [nb₁, nb₂...] and optional inputdelay = [d₁, d₂...]
"""
function arx(d::AbstractIdData, na, nb; inputdelay = ones(Int, size(nb)), λ = 0, estimator = \, stochastic = false)
    y, u, h = time1(output(d)), time1(input(d)), sampletime(d)
    @assert size(y, 2) == 1 "arx only supports single output."
    # all(nb .<= na) || throw(DomainError(nb,"nb must be <= na"))
    na >= 0 || throw(ArgumentError("na must be positive"))
    size(nb) == size(inputdelay) || throw(ArgumentError("inputdelay has to have the same structure as nb"))
    y_train, A = getARXregressor(vec(y), u, na, nb, inputdelay = inputdelay)
    w = ls(A, y_train, λ, estimator)
    a, b = params2poly2(w, na, nb, inputdelay = inputdelay)
    model = tf(b, a, h)
    if stochastic
        local Σ
        try
            Σ = parameter_covariance(y_train, A, w, λ)
        catch
            return minreal(model)
        end
        return minreal(TransferFunction(Particles, model, Σ))
    end

    return minreal(model)
end

"""
    residuals(ARX::TransferFunction, d::InputOutputData)

Calculates the residuals `v = Ay - Bu` of an ARX process and InputOutputData d. The length of the returned residuals is `length(d) - max(na, nb)`
# Example:
```jldoctest
julia> ARX = tf(1, [1, -1], 1)
TransferFunction{Discrete{Int64}, ControlSystems.SisoRational{Int64}}
  1
-----
z - 1

Sample Time: 1 (seconds)
Discrete-time transfer function model

julia> u = 1:5
1:5

julia> y = lsim(ARX, u, 1:5)[1][:]
5-element Vector{Float64}:
  0.0
  1.0
  3.0
  6.0
 10.0

julia> d = iddata(y, u)
InputOutput data of length 5 with 1 outputs and 1 inputs

julia> residuals(ARX, d)
4-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
```
"""
function residuals(ARX::TransferFunction, d::InputOutputData)
    size(ARX, 2) == ninputs(d) || throw(DomainError(d, "number of inputs $(ninputs(d)) does not match ARX model (expects $(size(ARX, 2)) inputs)"))

    w, a, b, inputdelay = params(ARX)
    na = length(a)
    nb = map(length, vec(b))

    y, A = getARXregressor(d, na, nb, inputdelay = inputdelay)
    ypred = A * w
    v = y .- ypred
    return v
end

"""
    predict(ARX::TransferFunction, d::InputOutputData)

One step ahead prediction for an ARX process.  The length of the returned prediction is `length(d) - max(na, nb)`
# Example:
```jldoctest
julia> predict(tf(1, [1, -1], 1), iddata(1:10, 1:10))
9-element Vector{Int64}:
  2
  4
  6
  8
 10
 12
 14
 16
 18
```
"""
function predict(ARX::TransferFunction, d::InputOutputData)
    size(ARX, 2) == ninputs(d) || throw(DomainError(d, "number of inputs $(ninputs(d)) does not match ARX model (expects $(size(ARX, 2)) inputs)"))

    w, a, b, inputdelay = params(ARX)
    na = length(a)
    nb = map(length, vec(b))

    y, A = getARXregressor(d, na, nb, inputdelay = inputdelay)
    ypred = A * w
    return ypred
end

"""
    ar(d::AbstractIdData, na; λ=0, estimator=\\, scaleB=false, stochastic=false)

Estimate an AR transfer function `G = 1/A`, the AR process is defined as `A(z⁻¹)y(t) = e(t)`

# Arguments:
- `d`: IdData, see [`iddata`](@ref)
- `na`: order of the model
- `λ`: `λ > 0` can be provided for L₂ regularization
- `estimator`: e.g. `\\,tls,irls,rtls`
- `scaleB`: Whether or not to scale the numerator using the variance of the prediction error.
- `stochastic`: if true, returns a transfer function with uncertain parameters represented by `MonteCarloMeasurements.Particles`.

# Example
```jldoctest
julia> N = 10000
10000

julia> e = [-0.2; zeros(N-1)] # noise e
10000-element Vector{Float64}:
[...]

julia> G = tf([1, 0], [1, -0.9], 1) # AR transfer function
TransferFunction{Discrete{Int64}, ControlSystems.SisoRational{Float64}}
   1.0z
----------
1.0z - 0.9

Sample Time: 1 (seconds)
Discrete-time transfer function model

julia> y = lsim(G, e, 1:N)[1][:] # Get output of AR transfer function from input noise e
10000-element Vector{Float64}:
[...]

julia> Gest = ar(iddata(y), 1) # Estimate AR transfer function from output y
TransferFunction{Discrete{Float64}, ControlSystems.SisoRational{Float64}}
          1.0z
-------------------------
1.0z - 0.8999999999999998

Sample Time: 1.0 (seconds)
Discrete-time transfer function model

julia> G ≈ Gest # Test if estimation was correct
true

julia> eest = lsim(1/Gest, y, 1:N)[1][:] # recover the input noise e from output y and estimated transfer function Gest
10000-element Vector{Float64}:
[...]

julia> isapprox(eest, e, atol = eps()) # input noise correct recovered
true
```
"""
function ar(d::AbstractIdData, na; λ = 0, estimator = \, scaleB = false, stochastic = false)
    noutputs(d) == 1 || throw(ArgumentError("Only 1d time series supported"))
    y = vec(time1(output(d)))
    na >= 1 || throw(ArgumentError("na must be positive"))
    y_train, A = getARregressor(y, na)
    w = ls(A, y_train, λ, estimator)
    a, b = params2poly(w, na)
    if scaleB
        b = √mean(abs2, y_train - A * w) * b
    end
    model = tf(b, vec(a), sampletime(d))
    if stochastic
        Σ = parameter_covariance(y_train, A, w, λ)
        return TransferFunction(Particles, model, Σ)
    end
    model
end

"""
    G, H, e = arxar(d::InputOutputData, na::Int, nb::Union{Int, Vector{Int}}, nd::Int)

Estimate the ARXAR model `Ay = Bu + v`, where `v = He` and `H = 1/D`, using generalized least-squares method. For more information see Söderström - Convergence properties of the generalized least squares identification method, 1974.

# Arguments:
- `d`: iddata
- `na`: order of A
- `nb`: number of coefficients in B, the order is determined by `nb + inputdelay - 1`. In MISO estimation it takes the form nb = [nb₁, nb₂...].
- `nd`: order of D

# Keyword Arguments:
- `H = nothing`: prior knowledge about the AR noise model
- `inputdelay = ones(Int, size(nb))`: optional delay of input, inputdelay = 0 results in a direct term, takes the form inputdelay = [d₁, d₂...] in MISO estimation
- `λ = 0`: `λ > 0` can be provided for L₂ regularization
- `estimator = \\`: e.g. `\\,tls,irls,rtls`, the latter three require `using TotalLeastSquares`
- `δmin = 10e-4`: Minimal change in the power of e, that specifies convergence.
- `iterations = 10`: maximum number of iterations.
- `verbose = false`: if true, more information is printed

# Example:
```
julia> N = 500
500

julia> sim(G, u) = lsim(G, u, 1:N)[1][:]
sim (generic function with 1 method)

julia> A = tf([1, -0.8], [1, 0], 1)
TransferFunction{Discrete{Int64}, ControlSystems.SisoRational{Float64}}
1.0z - 0.8
----------
   1.0z

Sample Time: 1 (seconds)
Discrete-time transfer function model

julia> B = tf([0, 1], [1, 0], 1)
TransferFunction{Discrete{Int64}, ControlSystems.SisoRational{Int64}}
1
-
z

Sample Time: 1 (seconds)
Discrete-time transfer function model

julia> G = minreal(B / A)
TransferFunction{Discrete{Int64}, ControlSystems.SisoRational{Float64}}
   1.0
----------
1.0z - 0.8

Sample Time: 1 (seconds)
Discrete-time transfer function model

julia> D = tf([1, 0.7], [1, 0], 1)
TransferFunction{Discrete{Int64}, ControlSystems.SisoRational{Float64}}
1.0z + 0.7
----------
   1.0z

Sample Time: 1 (seconds)
Discrete-time transfer function model

julia> H = 1 / D
TransferFunction{Discrete{Int64}, ControlSystems.SisoRational{Float64}}
   1.0z
----------
1.0z + 0.7

Sample Time: 1 (seconds)
Discrete-time transfer function model

julia> u, e = randn(N), randn(N)
[...]

julia> y, v = sim(G, u), sim(H * (1/A), e) # simulate process
[...]

julia> d = iddata(y.+ v, u, 1)
InputOutput data of length 500 with 1 outputs and 1 inputs

julia> na, nb , nd = 1, 1, 1
(1, 1, 1)

julia> Gest, Hest, res = arxar(d, na, nb, nd)
(G = TransferFunction{Discrete{Int64}, ControlSystems.SisoRational{Float64}}
   0.9987917259291642
-------------------------
1.0z - 0.7937837464682017

Sample Time: 1 (seconds)
Discrete-time transfer function model, H = TransferFunction{Discrete{Int64}, ControlSystems.SisoRational{Float64}}
          1.0z
-------------------------
1.0z + 0.7019519225937721

Sample Time: 1 (seconds)
Discrete-time transfer function model, e = [...]
```
"""
function arxar(d::InputOutputData, na::Int, nb::Union{Int, AbstractVector{Int}}, nd::Int;
    H::Union{TransferFunction, Nothing} = nothing,
    inputdelay      = ones(Int, size(nb)),
    δmin::Real      = 0.001,
    iterations::Int = 10,
    estimator       = \,
    verbose::Bool   = false,
    λ::Real         = 0
)
    # Input Checking
    na >= 0 || throw(ArgumentError("na($na) must be positive"))
    all(nb .>= 0 )|| throw(ArgumentError("nb($nb) must be positive"))
    nd >= 1 || throw(ArgumentError("nd($nd) must be positive"))
    iterations >= 1 || throw(DomainError("iterations($iterations) must be >0"))
	δmin > 0 || throw(ArgumentError("δmin($δmin) must be positive"))
    ninputs(d) == length(nb) || throw(ArgumentError("Length of nb ($(length(nb))) must equal number of input signals ($(ninputs(d)))"))

    Ts = sampletime(d)

    # 1. initialize H, GF and v to bring them to scope
    if H === nothing
        H = tf([0], [1], Ts)
        iter = 0 # initialization
    else
        iter = 1 # initial noisemodel is known, H is present
    end
    GF = tf([0], [1], Ts)
    v  = 0

    # Iterate
    eOld    = 0
    δ       = δmin + 1
    timeVec = timevec(d)
    sim(G,u) = lsim(G, u, timeVec)[1][:]
    while iter <= iterations && δ >= δmin
        # Filter input/output according to errormodel H, after initialization
        if iter > 0
            yF = sim(1/H, output(d)')
            if ninputs(d) == 1
                uF = sim(1/H, input(d)')
            else
                uF = fill(1.0, size(d.u))
                for i in 1:ninputs(d)
                    uF[i, :] = sim(1/H, d.u[i,:])
                end
            end
        else
            # pass unfiltered signal for initialization
            yF = copy(output(d))
            uF = copy(input(d))
        end
		dF = iddata(yF, uF, d.Ts)

		# 2. fit arx model
        GF = arx(dF, na, nb, estimator = estimator, λ = λ, inputdelay = inputdelay)

        # 3. Evaluate residuals
        v = residuals(GF, d)

        # 4. check if converged
        e = var(v)
        if eOld != 0
            δ = abs((eOld - e) / e)
        end
        verbose && println("iter: $iter, δ: $δ, e: $e")
        eOld = e

        # 5. estimate new noise model from residuals
        dH = iddata(v, d.Ts)
        H = ar(dH, nd, estimator = estimator)

        iter += 1
    end

    # total residuals e
    e = lsim(1/H, v, timevec(v, Ts))[1][:]
    return (G = GF, H = H, e = e)
end

function reversal_ls(A, y)
    n = size(A, 2)
    J = I(n)[:, end:-1:1]
    R = A'A
    R .= 1 / 2 .* (R + J * R' * J)
    R \ (A'y)
end

function ls(A, y, λ = 0, estimator = \)
    if λ == 0
        w = estimator(A, y)
    else
        w = estimator([A; λ * I], [y; zeros(size(A, 2))])
    end
    w
end

"""
G, Gn = plr(d::AbstractIdData,na,nb,nc; initial_order = 20)

Perform pseudo-linear regression to estimate a model on the form
`Ay = Bu + Cw`
The residual sequence is estimated by first estimating a high-order arx model, whereafter the estimated residual sequence is included in a second estimation problem. The return values are the estimated system model, and the estimated noise model. `G` and `Gn` will always have the same denominator polynomial.

`armax` is an alias for `plr`. See also [`pem`](@ref), [`ar`](@ref), [`arx`](@ref)
"""
function plr(d::AbstractIdData, na, nb, nc; initial_order = 20, method = :ls)
    y, u, h = time1(output(d)), time1(input(d)), sampletime(d)
    all(nb .<= na) || throw(DomainError(nb, "nb must be <= na"))
    na >= 1 || throw(ArgumentError("na must be positive"))
    # na -= 1
    y_train, A = getARXregressor(y, u, initial_order, initial_order)
    w1 = A \ y_train
    yhat = A * w1
    ehat = yhat - y_train
    ΔN = length(y) - length(ehat)
    y_train, A =
        getARXregressor(y[ΔN+1:end-1], [u[ΔN+1:end-1, :] ehat[1:end-1]], na, [nb; nc])
    if method == :wtls
        w = wtls_estimator(y[ΔN+1:end-1], na, nb + nc, [zeros(nb); ones(nc)])(A, y_train)
    elseif method == :rtls
        w = rtls(A, y_train)
    elseif method == :ls
        w = A \ y_train
    elseif method isa Function
        w = method(A, y_train)
    end
    a, b = params2poly(w, na, nb)
    model = tf(b, a, h)
    c = w[na+sum(nb)+1:end]
    noise_model = tf(c, a, h)
    model, noise_model
end

const armax = plr

"""
    model = arma(d::AbstractIdData, na, nc; initial_order=20, method=:ls)

Estimate a Autoregressive Moving Average model with `na` coefficients in the denominator and `nc` coefficients in the numerator.
Returns the model and the estimated noise sequence driving the system.

# Arguments:
- `d`: iddata
- `initial_order`: An initial AR model of this order is used to estimate the residuals
- `estimator`: A function `(A,y)->minimizeₓ(Ax-y)` default is `\` but another option is `wtls_estimator(1:length(y)-initial_order,na,nc,ones(nc))`

See also [`estimate_residuals`](@ref)
"""
function arma(d::AbstractIdData, na, nc; initial_order = 20, estimator = \)
    y, h = time1(output(d)), sampletime(d)
    all(nc .<= na) || throw(DomainError(nb, "nc must be <= na"))
    na >= 1 || throw(ArgumentError("na must be positive"))
    # na -= 1
    y_train, A = getARregressor(y, initial_order)
    # w1 = method isa Function ? method(A,y_train) : tls(A,y_train)
    w1 = A \ y_train # The above seems like it should work but tests indicate that it makes performance worse, maybe exactly because it removes too much correlation in the residuals
    yhat = A * vec(w1)
    ehat = y_train - yhat
    σ = √mean(abs2, ehat)
    ΔN = length(y) - length(ehat)
    y2 = @view y[ΔN:end-1]
    y_train, A = getARXregressor(y2, ehat, na, nc)
    w = estimator(A, y_train)
    a, b = params2poly(w, na, nc)
    # b[1] = 1
    b = σ * b
    model = tf(b, a, h)
    model
end

"""
    tfest(
        data::FRD,
        p0,
        link = log ∘ abs;
        freq_weight = sqrt(data.w[1]*data.w[end]),
        opt = BFGS(),
        opts = Optim.Options(
            store_trace       = true,
            show_trace        = true,
            show_every        = 1,
            iterations        = 100,
            allow_f_increases = false,
            time_limit        = 100,
            x_tol             = 0,
            f_tol             = 0,
            g_tol             = 1e-8,
            f_calls_limit     = 0,
            g_calls_limit     = 0,
        ),
    )

Fit a parametric transfer function to frequency-domain data.

# Arguments:
- `data`: An `FRD` onbject with frequency domain data.
- `p0`: Initial parameter guess. Can be a `NamedTuple` or `ComponentVector` with fields `b,a` specifying numerator and denominator as they appear in the call to `tf`, i.e., `(b = [1.0], a = [1.0,1.0,1.0])`. Can also be an instace of `TransferFunction`.
- `link`: By default, phase information is discarded in the fitting. To include phase, change to `link = log`.
- `freq_weight`: Apply weighting with the inverse frequency. The value determines the cutoff frequency before which the weight is constant, after which the weight decreases linearly. Defaults to the geometric mean of the smallest and largest frequency.
- `opt`: The Optim optimizer to use.
- `opts`: `Optim.Options` controlling the solver options.
"""
function tfest(
    data::FRD,
    p0::Union{NamedTuple, ComponentArray},
    link = log ∘ abs;
    freq_weight = sqrt(data.w[2] * data.w[end]),
    opt = BFGS(),
    opts = Optim.Options(
        store_trace       = true,
        show_trace        = true,
        show_every        = 1,
        iterations        = 10000,
        allow_f_increases = false,
        time_limit        = 100,
        x_tol             = 0,
        f_tol             = 0,
        g_tol             = 1e-8,
        f_calls_limit     = 0,
        g_calls_limit     = 0,
    ),
)

    ladr = @. link(data.r)
    if freq_weight > 0
        wv = 1 ./ (data.w .+ freq_weight)
    end

    function loss(p)
        a, b = p.a, p.b
        G = tf(b, a)
        mag = vec(freqresp(G, data.w))
        @. mag = link(mag) - ladr
        if freq_weight > 0
            mag .*= wv
        end
        mean(abs2, mag)
    end

    res = Optim.optimize(loss, ComponentVector(p0), opt, opts, autodiff = :forward)

    tf(res.minimizer.b, res.minimizer.a)
end

tfest(data, G::LTISystem, args...; kwargs...) = tfest(tfest(data)[1], G, args...; kwargs...)

function tfest(data::FRD, G::LTISystem, args...; kwargs...)
    ControlSystems.issiso(G) || throw(ArgumentError("Can only fit SISO model to FRD"))
    ControlSystems.isdiscrete(G) && throw(DomainError("Continuous-time model expected"))
    b, a = numvec(G)[], denvec(G)[]
    tfest(data, (; b = b, a = a), args...; kwargs...)
end

function tfest(data::FRD, basis::AbstractStateSpace;
    freq_weight = 1 ./ (data.w .+ data.w[2]),
    opt = BFGS(),
    metric::M = abs2,
    opts = Optim.Options(
        store_trace       = true,
        show_trace        = true,
        show_every        = 50,
        iterations        = 1000000,
        allow_f_increases = false,
        time_limit        = 100,
        x_tol             = 1e-5,
        f_tol             = 0,
        g_tol             = 1e-8,
        f_calls_limit     = 0,
        g_calls_limit     = 0,
    ),
) where M
    ω      = data.w
    Fs     = basis_responses(basis, ω, inverse=false)
    p0     = randn(length(Fs))
    resp_P = abs2.(data.r) .|> log
    Gmat   = reduce(hcat, Fs)
    Gθ     = similar(Gmat, size(Gmat, 1))

    function loss(p)
        mul!(Gθ, Gmat, p)
        c = zero(eltype(p))
        @inbounds for i in eachindex(ω)
            fp = log(abs2(Gθ[i]))
            c += metric(fp - resp_P[i]) * freq_weight[i]
        end
        c/length(ω)
    end

    res = Optim.optimize(
        loss,
        p0,
        opt,
        opts,
        # autodiff=:forward
    )
    F = sum_basis(basis, res.minimizer)
    if dcgain(F)[] < 0
        F = -F
        res.minimizer .= .- res.minimizer
    end
    F, res
end

"""
    arma_ssa(d::AbstractIdData, na, nc; L=nothing, estimator=\\, robust=false)

DOCSTRING

# Arguments:
- `d`: iddata
- `na`: number of denominator parameters
- `nc`: number of numerator parameters
- `L`: length of the lag-embedding used to separate signal and noise. `nothing` corresponds to automatic selection.
- `estimator`: The function to solve the least squares problem. Examples `\\,tls,irls,rtls`.
- `robust`: Use robust PCA to be resistant to outliers.
"""
function arma_ssa(d::AbstractIdData, na, nc; L = nothing, estimator = \, robust = false)
    y, h = time1(output(d)), sampletime(d)
    all(nc .<= na) || throw(DomainError(nc, "nc must be <= na"))
    na >= 1 || throw(ArgumentError("na must be positive"))

    L === nothing && (L = min(length(y) ÷ 2, 2000))
    H = hankel(y, L)
    if robust
        H, _, s = rpca(H)
        # yhat = unhankel(H)
    else
        s = svd(H)
    end
    yinds = 1:L-L÷4
    noiseinds = L-L÷4+1:L
    yhat = unhankel(s.U[:, yinds] * Diagonal(s.S[yinds]) * s.Vt[yinds, :])
    ehat = unhankel(s.U[:, noiseinds] * Diagonal(s.S[noiseinds]) * s.Vt[noiseinds, :])
    dhat = iddata(yhat, ehat, d.Ts)

    arma(dhat, na, nc, initial_order = L ÷ 5)
end

"""
    estimate_residuals(model, y)

Estimate the residuals driving the dynamics of an ARMA model.
"""
function estimate_residuals(model, y)
    y = time1(output(y))
    eest = zeros(length(y))
    na = length(denvec(model)[1, 1]) - 1
    nc = length(numvec(model)[1, 1])
    w = params(model)[1]
    for i = na:length(y)-1
        ϕ = [y[i:-1:i-na+1]; 0; eest[i-1:-1:i-nc+1]]
        yh = w'ϕ
        eest[i] = y[i+1] - yh
        ϕ[na+1] = eest[i]
        yh = w'ϕ
        eest[i] += y[i+1] - yh
    end
    eest
end

# Helper constructor to make a MISO system after MISO arx estimation
function ControlSystems.tf(
    b::AbstractVector{<:AbstractVector{<:Number}},
    a::AbstractVector{<:Number},
    h,
)
    tfs = map(b) do b
        tf(b, a, h)
    end
    hcat(tfs...)
end

"""
wtls_estimator(y,na,nb, σu=0)
Create an estimator function for estimation of arx models in the presence of measurement noise. If the noise variance on the input `σu` (model errors) is known, this can be specified for increased accuracy.
"""
function wtls_estimator(y, na, nb, σu = 0)
    y = output(y)
    rowQ = [diagm(0 => [ones(na); σu .* ones(nb); 1]) for i = 1:length(y)-(na)]
    Qaa, Qay, Qyy = rowcovariance(rowQ)
    (A, y) -> wtls(A, y, Qaa, Qay, Qyy)
end

"""
    a,b = params2poly(params,na,nb; inputdelay = zeros(Int, size(nb)))
Used to get numerator and denominator polynomials after arx fitting
"""
function params2poly(w, na, nb; inputdelay = zeros(Int, size(nb)))
    maxb = maximum(nb .+ inputdelay)
    a = [1; -w[1:na]]
    a = [a; zeros(max(0, maxb - na))] # if nb > na
    w = w[na+1:end]
    b = map(1:length(nb)) do i
        b = w[1:nb[i]]
        w = w[nb[i]+1:end]
        b = [zeros(inputdelay[i]); b; zeros(maxb - inputdelay[i] - nb[i])] # compensate for different nbs and delay
        b
    end
    return a, b
end
"""
    a,b = params2poly2(params,na,nb; inputdelay = ones(Int, size(nb)))
Used to get numerator and denominator polynomials after arx fitting. Updated version supporting a direct term.
"""
function params2poly2(w, na, nb; inputdelay = ones(Int, size(nb)))
    maxb = maximum(nb .+ inputdelay)
    a = [1; -w[1:na]]
    a = [a; zeros(max(0, maxb - na -1))] # if nb > na
    w = w[na+1:end]
    b = map(1:length(nb)) do i
        b = w[1:nb[i]]
        w = w[nb[i]+1:end]
        b = [zeros(inputdelay[i]); b; zeros(max(na+1, maxb) - inputdelay[i] - nb[i])] # compensate for different nbs and delay, as well as na > nb
        b
    end
    return a, b
end

function params2poly(w, na)
    a = [1; -w]
    b = [1; zeros(na)]
    return a, b
end

# """
# w, a, b = params(G::TransferFunction)
# w = [a;b]
# """
# function params(G::TransferFunction)
#     am, bm = -denvec(G)[1][2:end], numvec(G)[1]
#     wm     = [am; bm]
#     wm, am, bm
# end
"""
w, a, b, inputdelay = params(G::TransferFunction)
w = [a; vcat(b...)]

retrieve na and nb with:
na = length(a)
nb = map(length, vec(b))
"""
function params(G::TransferFunction)
    ams = vec(denvec(G))
    am = -ams[1][2:end]
    filter!(!iszero, am) # assumption, that no coefficient actually iszero

    bm = vec(numvec(G))
    inputdelay = length.(ams) .- length.(bm)
    filter!.(!iszero, bm)

    wm = [am; vcat(bm...) ]
    wm, am, bm, inputdelay
end

"""
    Σ = parameter_covariance(y_train, A, w, λ=0)
"""
function parameter_covariance(y_train, A, w, λ = 0)
    σ² = var(y_train .- A * w)
    iATA = if λ == 0
        inv(A'A)
    else
        ATA = A'A
        ATAλ = ATA + λ * I
        ATAλ \ ATA / ATAλ
    end
    iATA = (iATA + iATA') / 2
    Σ = σ² * iATA + sqrt(eps()) * Matrix(LinearAlgebra.I, size(iATA))
end

"""
TransferFunction(T::Type{<:AbstractParticles}, G::TransferFunction, Σ, N=500)

Create a `TransferFunction` where the coefficients are `Particles` from [`MonteCarloMeasurements.jl`](https://github.com/baggepinnen/MonteCarloMeasurements.jl) that can represent uncertainty.
# Example
```julia
using MonteCarloMeasurements
G       = ar(d,2,stochastic=true)
w       = exp10.(LinRange(-3,log10(π/Δt),100))
mag     = bode(Glsp,w)[1][:]
errorbarplot(w,mag,0.01; yscale=:log10, xscale=:log10, layout=3, subplot=1, lab="ls")

See full example [here](https://github.com/baggepinnen/MonteCarloMeasurements.jl/blob/master/examples/controlsystems.jl)
```
"""
function ControlSystems.TransferFunction(
    T::Type{<:MonteCarloMeasurements.AbstractParticles},
    G::TransferFunction,
    Σ::AbstractMatrix,
    N = 500,
)
    wm, am, bm, inputdelay = params(G)
    na = length(am)
    nb = map(length, vec(bm))
    if length(nb) == 1 ## SISO
        b = bm[1]
        isAR = b[1] == 1 && all(b[2:end] .== 0)
    else
        isAR = false # MISO -> ARX
    end

    if isAR && size(Σ, 1) < length(wm)
        p = T(N, MvNormal(wm[1:end-nb[1]], Σ))
        a, b = params2poly(p, na)
    else
        p = T(N, MvNormal(wm, Σ))
        a, b = params2poly2(p, na, nb, inputdelay = inputdelay)
    end
    arxtf = tf(b, a, G.Ts)
end

# function ControlSystems.TransferFunction(T::Type{<:MonteCarloMeasurements.AbstractParticles}, G::TransferFunction, p::AbstractMatrix)
#       wm, am, bm = params(G)
#       na,nb      = length(am), length(bm)
#       p          = T(p .+ wm')
#       a,b        = params2poly(p,na,nb)
#       arxtf      = tf(b,a,G.Ts)
# end

function DSP.filt(tf::ControlSystems.TransferFunction, y)
    b, a = numvec(tf)[], denvec(tf)[]
    filt(b, a, y)
end


text:
# Robust optimization
# 1. Function for MILP with Box Uncertainty and Budget of Uncertainty
# 2. Fucntion for Two-Stage Stochastic LP with Box Uncertainty

"""
    milpBoxBudget(; num_x, num_y, vec_min_y, vec_max_y, vec_c, vec_f, vec_b, mat_a,
            mat_b, mat_a_bar, mat_a_hat, mat_b_bar, vec_b_bar, vec_gammaCap)

Function for MILP with Box Uncertainty and Budget of Uncertainty
"""
function milp_box_budget(;
    num_x,
    num_y,
    vec_min_y,
    vec_max_y,
    vec_c,
    vec_f,
    vec_b,
    mat_a,
    mat_b,
    mat_a_bar,
    mat_a_hat,
    mat_b_bar,
    vec_b_bar,
    vec_gammaCap,
)
    m = length(mat_a_bar[:,1])
    n = length(mat_a_bar[1,:])
    for i = 1: m
        if vec_gammaCap[i] > n
            @error "vec_gammaCap[$i] is greater than m in nominal and hat data!"
        end
    end
    @info "Robust MILP with box uncertainty and budget."
    @debug "Input: vec_gammaCap = $vec_gammaCap."
    #
    model = Model(with_optimizer(GLPK.Optimizer))
    @variable(model, vec_y[1: num_y], Int)
    @variable(model, vec_x[1: num_x] >= 0)
    @objective(model, Min, (transpose(vec_c) * vec_x + transpose(vec_f) * vec_y)[1])
    @constraint(model, vec_y[1: num_y] .<= vec_max_y)
    @constraint(model, vec_y[1: num_y] .>= vec_min_y)
    @constraint(model, mat_a * vec_x + mat_b * vec_y .>= vec_b)
    # Transformation of Box Uncertainty.
    @variable(model, vec_lambda[1: m] >= 0)
    @variable(model, mat_mu[1: m, 1: num_x] >= 0)
    @variable(model, vec_z[1: num_x] >= 0)
    @constraint(model, - vec_z .<= vec_x)
    @constraint(model, vec_x .<= vec_z)
    for i = 1: m
        # sum() of variables without coefficients can be used directly
        @constraint(model, transpose(mat_a_bar[i, :]) * vec_x + vec_gammaCap[i] * vec_lambda[i] +
            sum(mat_mu[i, :]) + transpose(mat_b_bar[i, :]) * vec_y <= vec_b_bar[i])
        @constraint(model, vec_lambda[i] .+ hcat(mat_mu[i, :]) .>=  hcat(mat_a_hat[i, :]) .* vec_z[:])
    end
    optimize!(model)
    obj_result = objective_value(model)
    vec_result_y = value(vec_y)
    println("Result: vec_y = $vec_result_y.")
    vec_result_x = value(vec_x)
    println("Result: vec_x = $vec_result_x.")
    vec_result_z = value(vec_z)
    mat_result_mu = value(mat_mu)
    vec_result_lambda = value(vec_lambda)

    println("---------------------------------------------------------\n",
            "   2/2. Nominal Ending\n",
            "---------------------------------------------------------\n")

    return (
        obj_result,
        vec_result_y,
        vec_result_x,
        vec_result_z,
        mat_result_mu,
        vec_result_lambda,
    )
end

"""
Function for LP with Box Uncertainty and Budget of Uncertainty
"""
function lpBoxBudget(;
    num_x,
    vec_c,
    vec_b,
    mat_a,
    mat_a_bar,
    mat_a_hat,
    vec_b_bar,
    vec_gammaCap
)
    m = length(mat_a_bar[:,1])
    n = length(mat_a_bar[1,:])
    for i = 1: m
        if vec_gammaCap[i] > n
            println("Error: vec_gammaCap[$i] is greater than m in nominal and hat data!")
        end
    end
    println("---------------------------------------------------------\n",
            "    1/2. Robust MILP with Box Uncertainty and Budget\n",
            "---------------------------------------------------------\n")
    println("Input: vec_gammaCap = $vec_gammaCap.")
    #
    model = Model(with_optimizer(GLPK.Optimizer))
    @variable(model,
    vec_x[1: num_x] >= 0)
    @objective(model, Min, (transpose(vec_c) * vec_x)[1])
    # Transformation of Box Uncertainty.
    @variable(model, vec_lambda[1: m] >= 0)
    @variable(model, mat_mu[1: m, 1: num_x] >= 0)
    @variable(model, vec_z[1: num_x] >= 0)
    @constraint(model, - vec_z .<= vec_x)
    @constraint(model, vec_x .<= vec_z)
    for i = 1: m
        # sum() of variables without coefficients can be used directly
        @constraint(model, transpose(mat_a_bar[i, :]) * vec_x + vec_gammaCap[i] * vec_lambda[i] +
            sum(mat_mu[i, :]) <= vec_b_bar[i])
        @constraint(model, vec_lambda[i] .+ hcat(mat_mu[i, :]) .>=  hcat(mat_a_hat[i, :]) .* vec_z[:])
    end
    optimize!(model)
    obj_result = objective_value(model)
    vec_result_x = value(vec_x)
    # println("Result: vec_x = $vec_result_x.")
    vec_result_z = value(vec_z)
    mat_result_mu = value(mat_mu)
    vec_result_lambda = value(vec_lambda)
    println("---------------------------------------------------------\n",
            "   2/2. Nominal Ending\n",
            "---------------------------------------------------------\n")
    return (obj_result, vec_result_x, vec_result_z, mat_result_mu, vec_result_lambda)
end

"""
Fucntion for Two-Stage Stochastic MILP with Box Uncertainty
"""
function milpAdjustBox(;
    num_x,
    num_y,
    num_z,
    vec_min_y,
    vec_max_y,
    vec_c,
    vec_f,
    vec_g,
    vec_b,
    mat_a,
    mat_b,
    mat_d,
    mat_a_bar,
    mat_a_hat,
    mat_b_bar,
    mat_d_bar,
    vec_b_bar
)
    println("---------------------------------------------------------\n",
            "   1/2. Adjustable Robust LP with Box Uncertainty\n",
            "---------------------------------------------------------\n")
    model = Model(with_optimizer(GLPK.Optimizer))
    # 1. Standard LP
    @variable(model, vec_y[1: num_y], Int)
    @constraint(model, vec_y >= vec_min_y)
    @constraint(model, vec_max_y >=  vec_y)
    @variable(model, vec_x[1: num_x] >= 0)
    @variable(model, gamma)
    @variable(model, vec_alpha[1: num_z])
    @variable(model, vec_theta1[1: num_z])
    @objective(model, Min, (transpose(vec_f) * vec_y + transpose(vec_c) * vec_x)[1] + gamma)
    @constraint(model, mat_b * vec_y + mat_a * vec_x + mat_d * (vec_alpha + vec_theta1) .>= vec_b)
    # 2. Get rid of the uncertainty in objective function
    @variable(model, vec_beta[1: num_z])
    @constraint(model, gamma >= (transpose(vec_g) * (vec_alpha + vec_theta1))[1])
    @constraint(model, vec_theta1 .<= vec_beta)
    @constraint(model, - vec_beta .<= vec_theta1)
    @constraint(model, vec_alpha .>= - vec_theta1)  # z(zeta) must be greater than 0
    # 3. Transformation of Box Uncertainty.
    @variable(model, vec_theta2[1: num_z])
    vec_theta2_rep = (ones(num_y, num_z) .* vec_theta2')[:]  # Repeat every theta2 num_y times and then concatenate.
    @constraint(model, mat_b_bar * vec_y + mat_a_bar * vec_x + mat_a_hat * vec_theta2_rep +
        mat_d_bar * (vec_alpha + vec_theta1) .<= vec_b_bar)
    @constraint(model, vec_theta2_rep .<= vec_x)
    @constraint(model, - vec_x .<= vec_theta2_rep)
    # @constraint(model, [m = 1: num_y], sum([mat_a_bar[m, (p-1) * num_y + m] * vec_x[(p-1) * num_y + m] + vec_theta2[p]]
    #     for p = 1: 10)[1] + mat_b_bar[m, m] * vec_y[m] <= vec_b_bar[m])
    # @constraint(model, [m = 1: num_y, p = 1: num_z], - mat_a_hat[m, (p-1) * num_y + m] * vec_x[(p-1) * num_y + m]
    #     <= vec_theta2[p])
    # @constraint(model, [m = 1: num_y, p = 1: num_z], vec_theta2[p] <=
    #     mat_a_hat[m, (p-1) * num_y + m] * vec_x[(p-1) * num_y + m])
    # Solve the model
    optimize!(model)
    obj_result = objective_value(model)
    vec_result_y = value(vec_y)
    println("Result: vec_y = $vec_result_y.")
    vec_result_x = value(vec_x)
    println("Result: vec_x = $vec_result_x.")
    vec_result_theta1 = value(vec_theta1)
    vec_result_theta2 = value(vec_theta2)
    vec_result_gamma = value(gamma)
    println("---------------------------------------------------------\n",
            "   2/2. Nominal Ending\n",
            "---------------------------------------------------------\n")
    return (obj_result, vec_result_y, vec_result_x, vec_result_gamma, vec_result_theta1, vec_result_theta2)
end

"""
Fucntion for Two-Stage Stochastic LP with Box Uncertainty
"""
function lpAdjustBox(; num_x, m, num_z, vec_c, vec_g, vec_b,
    mat_a, mat_d, mat_a_bar, mat_a_hat, mat_d_bar, vec_b_bar)
    println("---------------------------------------------------------\n",
            "   1/2. Adjustable Robust LP with Box Uncertainty\n",
            "---------------------------------------------------------\n")
    model = Model(with_optimizer(GLPK.Optimizer))
    # 1. Standard LP
    @variable(model, vec_x[1: num_x] >= 0)
    @variable(model, gamma)
    @variable(model, vec_alpha[1: num_z])
    @variable(model, vec_theta1[1: num_z])
    @objective(model, Min, (transpose(vec_c) * vec_x)[1] + gamma)
    @constraint(model, mat_a * vec_x + mat_d * (vec_alpha + vec_theta1) .>= vec_b)
    # 2. Get rid of the uncertainty in objective function
    @variable(model, vec_beta[1: num_z])
    @constraint(model, gamma >= (transpose(vec_g) * (vec_alpha + vec_theta1))[1])
    @constraint(model, vec_theta1 .<= vec_beta)
    @constraint(model, - vec_beta .<= vec_theta1)
    @constraint(model, vec_alpha .>= - vec_theta1)  # z(zeta) must be greater than 0
    # 3. Transformation of Box Uncertainty.
    @variable(model, vec_theta2[1: num_z])
    vec_theta2_rep = (ones(m, num_z) .* vec_theta2')[:]  # Repeat every theta2 num_y times and then concatenate.
    @constraint(model, mat_a_bar * vec_x + mat_a_hat * vec_theta2_rep +
        mat_d_bar * (vec_alpha + vec_theta1) .<= vec_b_bar)
    @constraint(model, vec_theta2_rep .<= vec_x)
    @constraint(model, - vec_x .<= vec_theta2_rep)
    # Solve the model
    optimize!(model)
    obj_result = objective_value(model)
    # println("Result: obj = $(obj).")
    vec_result_x = value(vec_x)
    # println("Result: vec_x = $vec_result_x.")
    vec_result_theta1 = value(vec_theta1)
    vec_result_theta2 = value(vec_theta2)
    vec_result_gamma = value(gamma)
    println("---------------------------------------------------------\n",
            "   2/2. Nominal Ending\n",
            "---------------------------------------------------------\n")
    return (obj_result, vec_result_x, vec_result_gamma, vec_result_theta1, vec_result_theta2)
end

"""
Fucntion for LP with Box Uncertainty
"""
function lpBox(; num_x, m, num_z, vec_c, vec_b, mat_a, mat_a_bar, mat_a_hat, vec_b_bar)
    println("---------------------------------------------------------\n",
            "   1/2. Adjustable Robust LP with Box Uncertainty\n",
            "---------------------------------------------------------\n")
    model = Model(with_optimizer(GLPK.Optimizer))
    # 1. Standard LP
    @variable(model, vec_x[1: num_x] >= 0)
    @objective(model, Min, (transpose(vec_c) * vec_x)[1])
    @constraint(model, mat_a * vec_x .>= vec_b)
    # 2. Transformation of Box Uncertainty.
    @variable(model, vec_theta2[1: m])
    vec_theta2_rep = (ones(m, num_z) .* vec_theta2')[:]  # Repeat every theta2 num_y times and then concatenate.
    @constraint(model, mat_a_bar * vec_x + mat_a_hat * vec_theta2_rep .<= vec_b_bar)
    @constraint(model, vec_theta2_rep .<= vec_x)
    @constraint(model, - vec_x .<= vec_theta2_rep)
    # Solve the model
    optimize!(model)
    obj_result = objective_value(model)
    # println("Result: obj = $(obj).")
    vec_result_x = value(vec_x)
    # println("Result: vec_x = $vec_result_x.")
    vec_result_theta2 = value(vec_theta2)
    println("---------------------------------------------------------\n",
            "   2/2. Nominal Ending\n",
            "---------------------------------------------------------\n")
    return (obj_result, vec_result_x, vec_result_theta2)
end

ext: jl


text:
function save_mesh(::Val{:PLOT3D}, io::IO, m::AirfoilMesh; cfl3d::Bool=false)
    f = FortranFile(io)
    write(f, Int32(1))
    NI, NJ = size(m.grid)
    if cfl3d
        write(f, Int32[2, NI, NJ])
        blk = Array{Float64, 4}(undef, 2, NI, NJ, 3)
        for k in 1:NJ
            for j in 1:NI
                blk[:,j,k,1] .= m.grid[j,k][1] # x
                blk[:,j,k,3] .= m.grid[j,k][2] # z
            end
        end
        blk[1,:,:,2] .= 1 # y
        blk[2,:,:,2] .= 0
        write(f, blk)
    else
        write(f, Int32[NI, NJ])
        blk = Array{Float64, 3}(undef, NI, NJ, 2)
        for j in 1:NJ
            for i in 1:NI
                blk[i,j,1] = m.grid[i,j][1]
                blk[i,j,2] = m.grid[i,j][2]
            end
        end
        write(f, blk)
    end
    return nothing
end


text:
"The perfect sidekick to your scientific inquiries"
module DrWatson
import Pkg, LibGit2

const PATH_SEPARATOR = joinpath("_", "_")[2]

# Misc functions for kw-macros
convert_to_kw(ex::Expr) = Expr(:kw,ex.args...)
convert_to_kw(ex) = error("invalid keyword argument syntax \"$ex\"")

# Pure Julia implementation
include("project_setup.jl")
include("naming.jl")
include("saving_tools.jl")

using UnPack
export @pack!, @unpack

# Functionality that saves/loads
using FileIO
export save, load
export wsave, wload

_wsave(filename, data...; kwargs...) = FileIO.save(filename, data...; kwargs...)

"""
    wsave(filename, data...; kwargs...)

Save `data` at `filename` by first creating the appropriate paths.
Default fallback is `FileIO.save`. Extend this for your types
by extending `DrWatson._wsave(filename, data::YourType...; kwargs...)`.
"""
function wsave(filename, data...; kwargs...)
    mkpath(dirname(filename))
    return _wsave(filename, data...; kwargs...)
end

"Currently equivalent with `FileIO.load`."
wload(data...; kwargs...) = FileIO.load(data...; kwargs...)

include("saving_files.jl")
include("dict_list.jl")

# Functionality that requires Dataframes and other heavy dependencies:
using Requires
function __init__()
    @require DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" begin
        include("result_collection.jl")
    end
end

# Update messages
const display_update = true
const update_version = "2.0.0"
const update_name = "update_v$update_version"
if display_update
if !isfile(joinpath(@__DIR__, update_name))
printstyled(stdout,
"""
\nUpdate message: DrWatson v$update_version

In this new major release, the following breaking changes have occured:
1. DrWatson now uses, and suggests using, JLD2.jl instead of BSON.jl
   for saving files.
2. The behavior of `savename` with respect to rounding has changed,
   see its docstring for more. (no longer is `1.0` output as `1` in `savename`)
\n
"""; color = :light_magenta)
touch(joinpath(@__DIR__, update_name))
end
end

end


text:
using LinearAlgebra
using GFlops
using Memoize
using MKL
using BenchmarkTools
BLAS.get_config()

(typeof ∘ typeof ∘ typeof ∘ typeof ∘ typeof ∘ typeof)(1)

@gflops sum(rand(10, 10))

abstract type AbstractArray{T,N} end

function eltype(::Type{<:AbstractArray{T}}) where {T}
    T
end

AbstractArray{Float64,4}

Type{Array}

eltype(AbstractArray{Real,4})

typeof(Type{Real})

Type{Real}

using HCubature
using Plots
ENV["GKSwstype"] = "100"
plot(x -> x + 87, 5, 10)|>display
popdisplay()
scatter(rand(20000), rand(20000)) |> display

isa(1.0,Type{Float64})
begin
    return 1
    return 0
end

@memoize function fact(n)
           n < 0 && error("n must be non-negative")
           n <2  && return 1
           abs(fact(n-1)-fact(n-2))-n
end

function fact2(n)
        n < 0 && error("n must be non-negative")
        n <2  && return 2
        abs(fact2(n-1)-fact2(n-2))-n
end

@time fact2(20)

@time fact(50)

testv=rand(0:BigInt(10),1000)

@benchmark fact.($testv)

@benchmark fact2.($testv)

memoize_cache(fact)


text:
mutable struct ParetoCandidate{U, V, W}
    point::AbstractArray{U}
    parameter::AbstractArray{U}
    iter::V
    threshold::W
end

point(x::ParetoCandidate) = x.point
parameter(x::ParetoCandidate) = x.parameter
iter(x::ParetoCandidate) = x.iter
threshold(x::ParetoCandidate) = x.threshold

function update!(p::ParetoCandidate, point, parameter, iter, threshold)
    p.point .= point
    p.parameter .= parameter
    p.iter = iter
    p.threshold = threshold
    return
end

function update!(p::ParetoCandidate, y::ParetoCandidate)
    p.point .= point(y)
    p.parameter .= parameter(y)
    p.iter = iter(y)
    p.threshold = threshold(y)
    return
end

struct WeightedSum <: AbstractScalarizationMethod
    w::Union{UniformScaling, AbstractArray}
    f::Function
end

"""
    WeightedSum()
    WeightedSum(weights, function)

Scalarize the multi-objective optimization via a weighted sum such that the
objective becomes `sum(w_i ~f_i(x))`.
"""
WeightedSum() = WeightedSum(I, x->identity(x))

(w::WeightedSum)(x::ParetoCandidate) = sum(w.w*w.f(point(x)))

struct GoalProgramming <: AbstractScalarizationMethod
    n::Function
    f::Function
end

"""
    GoalProgramming()
    GoalProgramming(norm, function)

Scalarize the multi-objective optimization via a goal programming such that the
objective becomes `norm(f(x), p)`.
"""
GoalProgramming() = GoalProgramming(x->norm(x, 2), x->identity(x))

(g::GoalProgramming)(x::ParetoCandidate) = g.n(g.f(point(x)))

"""
    WeightedExponentialSum()
    WeightedExponentialSum(weights, function, p)

Scalarize the multi-objective optimization via a goal programming such that the
objective becomes `sum(w_i*f(x)_i)^p`.
"""
struct WeightedExponentialSum <: AbstractScalarizationMethod
    w::Union{UniformScaling, AbstractArray}
    f::Function
    p::Real
end

WeightedExponentialSum() = WeightedExponentialSum(I, identity, 2)

(w::WeightedExponentialSum)(x::ParetoCandidate) = sum(w.w*w.f(point(x)).^w.p)

mutable struct ParetoFront{S}
    candidates::AbstractArray{ParetoCandidate}
    scalarization::S
end

candidates(x::ParetoFront) = x.candidates

function ParetoFront(n::Int64; scalarization::AbstractScalarizationMethod = WeightedSum())
    candidates = Array{ParetoCandidate}(undef, n)
    return ParetoFront(candidates, scalarization)
end

(x::ParetoFront)(y::ParetoCandidate) = x.scalarization(y)

function assert_dominance(x::ParetoFront, y::ParetoFront)
    [x(cx) < y(cy) for (cx, cy) in zip(candidates(x), candidates(y))]
end

function conditional_add!(x::ParetoFront, y::ParetoFront)
    for (i, c) in enumerate(assert_dominance(y, x))
        c ? update!(x.candidates[i], y.candidates[i]) : nothing
    end
end

function set_candidate!(x::ParetoFront, idx, point, parameter, iter, threshold)
    if isdefined(candidates(x), idx)
        update!(x.candidates[idx], point, parameter, iter, threshold)
    else
        x.candidates[idx] = ParetoCandidate(point, parameter, iter, threshold)
    end
    return
end

Base.getindex(x::ParetoFront, idx) = getindex(x.candidates, idx)
iter(x::ParetoFront) =  maximum(iter.(candidates(x)))
threshold(x::ParetoFront) = minimum(threshold.(candidates(x)))


text:
"
    Abstract environment for an infinite state
    distinct from finite, because we have to recalculate everything when the state changes
"
abstract type AbstractInfEnv <: Cache end;

function leftenv(pars::AbstractInfEnv,pos::Int,state)
    check_recalculate!(pars,state)
    pars.lw[pos,:]
end

function rightenv(pars::AbstractInfEnv,pos::Int,state)
    check_recalculate!(pars,state)
    pars.rw[pos,:]
end

function leftenv(pars::AbstractInfEnv,row::Int,col::Int,state)
    check_recalculate!(pars,state)
    pars.lw[row,col]
end

function rightenv(pars::AbstractInfEnv,row::Int,col::Int,state)
    check_recalculate!(pars,state)
    pars.rw[row,col]
end

function poison!(pars::AbstractInfEnv)
    pars.dependency = similar(pars.dependency);
end

function check_recalculate!(pars::AbstractInfEnv,state)
    if !(pars.dependency === state)
        lock(pars.lock);

        #we have acquired the lock; maybe state has already been updated?
        !(pars.dependency === state) && recalculate!(pars,state)

        unlock(pars.lock);
    end
end


text:
using Test
using PyCall
using DFTK

ase = pyimport_e("ase")
if !ispynull(ase)
    py"""
    import ase
    """

    @testset "ase.Atoms -> DFTK lattice / atoms -> ase.Atoms" begin
        py"""
        import ase
        cell = [[3.21, 0.0, 0.0], [-1.605, 2.7799415461480477, 0.0], [0.0, 0.0, 5.21304]]
        positions = [[0.0, 0.0, 0.0], [0.0, 1.85329, 2.60652]]
        atoms = ase.Atoms(symbols='Mg2', pbc=True, cell=cell, positions=positions)
        """

        lattice = load_lattice(py"atoms")
        @test lattice[:, 1] ≈ DFTK.units.Ǎ * [3.21, 0.0, 0.0]
        @test lattice[:, 2] ≈ DFTK.units.Ǎ * [-1.605, 2.7799415461480477, 0.0]
        @test lattice[:, 3] ≈ DFTK.units.Ǎ * [0.0, 0.0, 5.21304]

        @test load_lattice(py"atoms") == load_lattice(py"atoms.cell")

        atoms = load_atoms(py"atoms")
        @test length(atoms) == 1
        @test all(at isa ElementCoulomb for (at, positions) in atoms)
        @test atoms[1][1].symbol == :Mg
        @test atoms[1][2][1] ≈ [0, 0, 0]
        @test atoms[1][2][2] ≈ [1/3, 2/3, 1/2] atol=1e-5

        newatoms = ase_atoms(lattice, atoms)
        @test all(py"$newatoms.cell == atoms.cell")
        @test all(py"$newatoms.symbols == atoms.symbols")
        @test py"$newatoms.get_positions()" ≈  py"atoms.get_positions()"
    end
end


text:
"""
    abstract type AbstractGrid{T,N} <: AbstractArray{T,N}

Grids are arrays of points.
"""
abstract type AbstractGrid{T,N} <: AbstractArray{T,N}
end

const AbstractGrid1d{T <: Real} = AbstractGrid{T,1}
dimension(grid::AbstractGrid{T}) where T = DomainSets.dimension_type(T)
prectype(::Type{G}) where {G<:AbstractGrid} = float(eltype(eltype(G)))
prectype(g::AbstractGrid) = prectype(typeof(g))

@propagate_inbounds function getindex(grid::AbstractGrid{T,1}, i::Int) where {T}
    checkbounds(grid, i)
    unsafe_getindex(grid, i)
end

export resize
"""
    resize(grid::AbstractGrid, dims...)

Create a grid of same structure (such as product structure)
but with different points in the different dimensions.
"""
resize(grid::AbstractGrid{T}, dims...) where {T} = similargrid(grid, T, dims...)

"""
    hasextension(grid::AbstractGrid)

Is it possible to use the `resize` function.
See also [`resize`](@ref)
"""
hasextension(::AbstractGrid) = false

"""
    iscomposite(grid::AbstractGrid)

Does the grid consist of multiple grids such as e.g. ProductGrid?
Used for pretty printing.
"""
iscomposite(::AbstractGrid) = false


text:
module CmdStan

using Compat, NamedArrays, Pkg, Documenter, DelimitedFiles, Unicode, Statistics

"""
The directory which contains the cmdstan executables such as `bin/stanc` and
`bin/stansummary`. Inferred from the environment variable `JULIA_CMDSTAN_HOME` or `ENV["JULIA_CMDSTAN_HOME"]`
when available.

If these are not available, use `set_cmdstan_home!` to set the value of CMDSTAN_HOME.

Example: `set_cmdstan_home!(homedir() * "/Projects/Stan/cmdstan/")`

Executing `versioninfo()` will display the value of `JULIA_CMDSTAN_HOME` if defined.
"""
CMDSTAN_HOME=""

function __init__()
  global CMDSTAN_HOME = if isdefined(Main, :JULIA_CMDSTAN_HOME)
    Main.JULIA_CMDSTAN_HOME
  elseif haskey(ENV, "JULIA_CMDSTAN_HOME")
    ENV["JULIA_CMDSTAN_HOME"]
  elseif haskey(ENV, "CMDSTAN_HOME")
    ENV["CMDSTAN_HOME"]
  else
    @warn("Environment variable CMDSTAN_HOME not set. Use set_cmdstan_home!.")
    ""
  end
end

"""Set the path for the `CMDSTAN_HOME` environment variable.

Example: `set_cmdstan_home!(homedir() * "/Projects/Stan/cmdstan/")`
"""
set_cmdstan_home!(path) = global CMDSTAN_HOME=path

include("main/stanmodel.jl")
include("main/stancode.jl")

# preprocessing

include("utilities/update_model_file.jl")
include("utilities/create_r_files.jl")
include("utilities/create_cmd_line.jl")

# run cmdstan

include("utilities/parallel.jl")

# used in postprocessing

include("utilities/read_samples.jl")
include("utilities/read_variational.jl")
include("utilities/read_diagnose.jl")
include("utilities/read_optimize.jl")
include("utilities/convert_a3d.jl")

# type definitions

include("types/sampletype.jl")
include("types/optimizetype.jl")
include("types/diagnosetype.jl")
include("types/variationaltype.jl")

export

# from this file
set_cmdstan_home!,
CMDSTAN_HOME,

# From stanmodel.jl
Stanmodel,

# From stancode.jl
stan,

# From sampletype.jl
Sample,

# From optimizetype.jl
Optimize,

# From diagnosetype.jl
Diagnose,

# From variationaltype.jl
Variational

end # module


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Functions to convert anomalies related to the orbit.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export M_to_E, M_to_f
export E_to_f, E_to_M
export f_to_E, f_to_M

################################################################################
#                              From Mean Anomaly
################################################################################

"""
    M_to_E(e::Number, M::Number, tol::Number = 1e-10)

Compute the eccentric anomaly (0,2π) \\[rad] given the eccentricity `e` and the
mean anomaly `M` [rad]. This function uses the Newton-Raphson algorithm and the
tolerance to accept the solution is `tol`.

"""
@inline function M_to_E(e::Number, M::Number, tol::Number = 1e-10)
    # Compute the eccentric anomaly using the Newton-Raphson method.
    # ==============================================================

    # Make sure that M is in the interval [0,2π].
    M = mod(M,2π)

    # Initial guess.
    #
    # See [1, p. 75].
    E = (M > π) ? M - e : M + e

    sin_E, cos_E = sincos(E)

    # Newton-Raphson iterations.
    while ( abs(E - e*sin_E - M) > tol )
        E = E - (E - e*sin_E - M)/(1-e*cos_E)

        sin_E, cos_E = sincos(E)
    end

    # Return the eccentric anomaly in the interval [0, 2π].
    mod(E, 2π)
end

"""
    M_to_f(e::Number, M::Number, tol::Number = 1e-10)

Compute the true anomaly (0,2π) \\[rad] given the eccentricity `e` and the mean
anomaly `M` [rad]. This function uses the Newton-Raphson algorithm and the
tolerance to accept the solution is `tol`.

"""
@inline function M_to_f(e::Number, M::Number, tol::Number = 1e-10)
    # Compute the eccentric anomaly.
    E = M_to_E(e, M, tol)

    # Compute the true anomaly in the interval [0,2π].
    E_to_f(e,E)
end

################################################################################
#                            From Eccentric Anomaly
################################################################################

"""
    E_to_f(e::Number, E::Number)

Compute the true anomaly (0,2π) \\[rad] given the eccentricity `e` and the
eccentric anomaly `E` [rad].

"""
@inline function E_to_f(e::Number, E::Number)
    sin_Eo2, cos_Eo2 = sincos(E/2)

    # Compute the true anomaly in the interval [0, 2*π].
    mod( 2atan(sqrt(1+e)*sin_Eo2, sqrt(1-e)*cos_Eo2) , 2π )
end

"""
    E_to_M(e::Number, E::Number)

Compute the mean anomaly (0,2π) \\[rad] given the eccentricity `e` and the
eccentric anomaly `E` [rad].

"""
E_to_M(e::Number, E::Number) = mod(E - e*sin(E), 2π)

################################################################################
#                              From True Anomaly
################################################################################

"""
    f_to_E(e::Number,f::Number)

Compute the eccentric anomaly (0,2π) \\[rad] given the eccentricity `e` and
the true anomaly `f` [rad].

"""
@inline function f_to_E(e::Number, f::Number)
    sin_fo2, cos_fo2 = sincos(f/2)

    mod( 2atan(sqrt(1-e)*sin_fo2, sqrt(1+e)*cos_fo2), 2π )
end

"""
    f_to_E(orb::Orbit)

Compute the eccentric anomaly (0,2π) \\[rad] given the orbit `orb` (see
`Orbit`).

"""
f_to_E(orb::Orbit) = f_to_E(orb.e, orb.f)

"""
    f_to_M(e::Number, f::Number)

Compute the mean anomaly (0,2π) \\[rad] given the eccentricity `e` and the
true anomaly `f` [rad].

"""
@inline function f_to_M(e::Number, f::Number)
    # Compute the eccentric anomaly.
    E = f_to_E(e, f)

    # Compute the true anomaly in the interval [0, 2π].
    E_to_M(e,E)
end

"""
    f_to_M(orb::Orbit)

Compute the mean anomaly (0,2π) \\[rad] given the orbit `orb` (see `Orbit`).

"""
f_to_M(orb::Orbit) = f_to_M(orb.e, orb.f)


text:
abstract type InputMapping end
struct InputAllToAll <: InputMapping
end
struct InputAllToFirstN <: InputMapping
  n::Int
end
struct InputDiag <: InputMapping
end
struct InputMatrices{T<:AbstractFloat} <: InputMapping
  matrices::Dict{Symbol, Matrix{T}}
end

function generate_sens_matrices(im::InputMapping, s_in, n_total; T=Float32)
  _w = zeros(T, s_in, n_total)
  _E = zeros(T, s_in, n_total)
  for (src,dst) in _input_idxs(im, s_in, n_total)
    add_synapse!(src, dst, _w, _E)
  end
  _w, _E
end

function generate_sens_matrices(im::InputMatrices, s_in, n_total; T=Float32)
  w = im.matrices[:w_sens]
  E = im.matrices[:E_sens]
  w, E
end

_input_idxs(m::InputAllToAll, s_in, n_total)    = [(src, dst) for dst in 1:n_total, src in 1:s_in]
_input_idxs(m::InputAllToFirstN, s_in, n_total) = [(src, dst) for dst in 1:m.n, src in 1:s_in]
_input_idxs(m::InputDiag, s_in, n_total)        = [(src, src) for src in 1:s_in]

get_n_in(m::InputMapping, n_total) = n_total
get_n_in(m::InputAllToFirstN, n_total) = m.n

abstract type OutputMapping end

struct OutputAll <: OutputMapping
end
struct OutputLastN <: OutputMapping
  n::Int
end
struct OutputIdxs <: OutputMapping
  idxs::Vector{Int}
end
_output_idxs(m::OutputAll, n_total)   = 1:n_total
_output_idxs(m::OutputLastN, n_total) = n_total-m.n+1:n_total
_output_idxs(m::OutputIdxs, n_total)  = m.idxs

get_n_out(m::OutputAll, n_total) = n_total
get_n_out(m::OutputLastN, n_total) = m.n
get_n_out(m::OutputIdxs, n_total) = length(m.idxs)

abstract type SynsMapping end

struct SynsFullyConnected <: SynsMapping
end

struct SynsNCP <: SynsMapping
  n_sensory::Int
  n_inter::Int
  n_command::Int
  n_motor::Int

  sensory_inter::Int
  inter_command::Int
  command_command::Int
  command_motor::Int
end

struct SynsMatrices{T<:AbstractFloat} <: SynsMapping
  matrices::Dict{Symbol, Matrix{T}}
end

SynsNCP(; n_sensory=2, n_inter=5, n_command=5, n_motor=1,
  sensory_inter=2, inter_command=3, command_command=2, command_motor=2
) = SynsNCP(n_sensory, n_inter, n_command, n_motor,
    sensory_inter, inter_command, command_command, command_motor)

function generate_syns_matrices(sm::SynsMapping, n_total::Integer; T=Float32)
  _w = zeros(T, n_total, n_total)
  _E = zeros(T, n_total, n_total)
  for (src,dst) in _syns_idxs(sm, n_total)
    add_synapse!(src, dst, _w, _E)
  end
  _w, _E
end

function generate_syns_matrices(sm::SynsMatrices, n_total::Integer; T=Float32)
  w = sm.matrices[:w_syns]
  E = sm.matrices[:E_syns]
  w, E
end

_syns_idxs(m::SynsFullyConnected, n_total) = [(src, dst) for dst in 1:n_total, src in 1:n_total]

function _syns_idxs(m::SynsNCP, n_total)
  idxs = Vector{Tuple{Int,Int}}()

  n_sensory       = m.n_sensory
  n_inter         = m.n_inter
  n_command       = m.n_command
  n_motor         = m.n_motor
  sensory_inter   = m.sensory_inter
  inter_command   = m.inter_command
  command_command = m.command_command
  command_motor   = m.command_motor

  #n_neurons = n_sensory + n_inter + n_command + n_motor
  n_neurons = n_total

  s = 1
  r_sens = s:s+n_sensory-1
  s += n_sensory
  r_inter = s:s+n_inter-1
  s += n_inter
  r_command = s:s+n_command-1
  s += n_command
  r_motor = s:s+n_motor-1

  w_syns = zeros(Float32, n_neurons, n_neurons)
  E_syns = zeros(Float32, n_neurons, n_neurons)

  connect_layers!(r_sens,    r_inter,   w_syns, E_syns; conns=sensory_inter)
  connect_layers!(r_inter,   r_command, w_syns, E_syns; conns=inter_command)
  connect_layers!(r_command, r_command, w_syns, E_syns; conns=command_command, handle_unconnected=false)
  connect_layers!(r_command, r_motor,   w_syns, E_syns; conns=command_motor)

  for dst in 1:size(w_syns,2)
    for src in 1:size(w_syns,1)
      w_syns[src,dst] == 0 && continue
      push!(idxs, (src,dst))
    end
  end

  idxs
end

mutable struct WiringConfig{T<:AbstractFloat, IM<:InputMapping, SM<:SynsMapping, OM<:OutputMapping}
  s_in::Int
  n_total::Int
  matrices::Dict{Symbol, Matrix{T}}
  input_mapping::IM
  syns_mapping::SM
  output_mapping::OM
end

WiringConfig(s_in::Integer, im::InputMapping, sm::SynsNCP, om::OutputMapping) =
  WiringConfig(s_in, sm.n_sensory+sm.n_inter+sm.n_command+sm.n_motor, im, sm, om)

function WiringConfig(s_in::Integer, n_total::Integer, im::InputMapping, sm::SynsMapping, om::OutputMapping; T=Float32)
  w_sens, E_sens = generate_sens_matrices(im, s_in, n_total; T)
  w_syns, E_syns = generate_syns_matrices(sm, n_total; T)

  matrices = Dict(
    :w_sens => w_sens,
    :E_sens => E_sens,
    :w_syns => w_syns,
    :E_syns => E_syns,
  )

  WiringConfig(s_in, n_total, matrices, im, sm, om)
end

function Base.getproperty(m::WiringConfig{T, IM, SM, OM}, s::Symbol) where {T<:AbstractFloat, IM<:InputMapping, SM<:SynsMapping, OM<:OutputMapping}
  if s === :s_in
    return getfield(m, :s_in)
  elseif s === :n_total
    return getfield(m, :n_total)
  elseif s === :n_out
    return get_n_out(m.output_mapping, m.n_total)
  elseif s === :input_mapping
    return getfield(m, :input_mapping)
  elseif s === :output_mapping
    return getfield(m, :output_mapping)
  elseif s === :syns_mapping
    return getfield(m, :syns_mapping)
  elseif s === :matrices
    return getfield(m, :matrices)
  else
    return getfield(m, s)
  end
end

default_weight() = 1
random_polarity(p=[-1,1]) = p[rand(1:length(p))]

function add_synapse!(src, dst, w, E)
  w[src,dst] = default_weight()
  #E[src,dst] = random_polarity([-50,30])
  E[src,dst] = random_polarity()
end

function connect_layers!(srcs, dsts, w, E; conns=Inf, handle_unconnected=true)
  for src in srcs
    n = 0
    for dst in dsts
      n > conns && break
      add_synapse!(src, dst, w, E)
      n += 1
    end
  end

  handle_unconnected == false && return

  mean_out = Int(ceil(mean([sum(w[i,dsts]) for i in srcs])))
  no_in = Int[]
  for dst in dsts
    sum(w[srcs,dst]) == 0 && push!(no_in, dst)
  end
  for dst in no_in
    for _ in 1:mean_out
      src = rand(srcs)
      add_synapse!(src, dst, w, E)
    end
  end
end

# struct Wiring{T<:AbstractFloat}
#   s_in::Int
#   n_total::Int
#   n_in::Int
#   n_out::Int

#   n_sensory::Int
#   n_inter::Int
#   n_command::Int
#   n_motor::Int

#   matrices::Dict{Symbol, Matrix{T}}
# end

# Wiring(s_in, n_total, n_in, n_out, matrices; n_sensory=0, n_inter=0, n_command=0, n_motor=0) =
#   Wiring(s_in, n_total, n_in, n_out, n_sensory, n_inter, n_command, n_motor, matrices)

# FullyConnected(s_in::Int, n_total::Int; T::DataType=Float32) = FullyConnected(s_in, n_total, n_total, n_total; T)
# function FullyConnected(s_in, n_total, n_in, n_out; T::DataType=Float32)
#   @assert n_out <= n_total
#   @assert n_in <= n_total

#   n_neurons = n_total

#   w_sens = zeros(T, s_in, n_neurons)
#   E_sens = zeros(T, s_in, n_neurons)
#   w_syns = zeros(T, n_neurons, n_neurons)
#   E_syns = zeros(T, n_neurons, n_neurons)

#   for s in 1:s_in
#     for d in 1:n_in
#       add_synapse!(s, d, w_sens, E_sens)
#     end
#   end
#   for s in 1:n_neurons
#     for d in 1:n_neurons
#       add_synapse!(s, d, w_syns, E_syns)
#     end
#   end

#   matrices = Dict(
#     :w_sens => w_sens,
#     :E_sens => E_sens,
#     :w_syns => w_syns,
#     :E_syns => E_syns,
#   )
#   Wiring(s_in, n_total, n_in, n_out, matrices)
# end

# function FullSensNCPWiring(s_in::Integer; T::DataType=Float32,
#                 n_sensory=2, n_inter=3, n_command=5, n_motor=1,
#                 sensory_inter=2, inter_command=2, command_command=2, command_motor=3)

#   FullSensNCPWiring(s_in, n_sensory, n_motor; T, n_sensory, n_inter, n_command, n_motor,
#                 sensory_inter, inter_command, command_command, command_motor)
# end
# function FullSensNCPWiring(s_in::Integer, n_in::Integer, n_out::Integer; T::DataType=Float32,
#                 n_sensory=2, n_inter=3, n_command=5, n_motor=1,
#                 sensory_inter=2, inter_command=2, command_command=2, command_motor=3)

#   n_neurons = n_sensory + n_inter + n_command + n_motor
#   n_total = n_neurons

#   r_in = 1:s_in

#   s = 1
#   r_sens = s:s+n_sensory-1
#   s += n_sensory
#   r_inter = s:s+n_inter-1
#   s += n_inter
#   r_command = s:s+n_command-1
#   s += n_command
#   r_motor = s:s+n_motor-1

#   w_sens = zeros(T, s_in, n_neurons)
#   E_sens = zeros(T, s_in, n_neurons)
#   w_syns = zeros(T, n_neurons, n_neurons)
#   E_syns = zeros(T, n_neurons, n_neurons)

#   connect_layers!(r_in, r_sens, w_sens, E_sens; handle_unconnected=false)

#   connect_layers!(r_sens,    r_inter,   w_syns, E_syns; conns=sensory_inter)
#   connect_layers!(r_inter,   r_command, w_syns, E_syns; conns=inter_command)
#   connect_layers!(r_command, r_command, w_syns, E_syns; conns=command_command, handle_unconnected=false)
#   connect_layers!(r_command, r_motor,   w_syns, E_syns; conns=command_motor)

#   matrices = Dict(
#     :w_sens => w_sens,
#     :E_sens => E_sens,
#     :w_syns => w_syns,
#     :E_syns => E_syns,
#   )

#   Wiring(s_in, n_total, n_in, n_out, matrices; n_sensory, n_inter, n_command, n_motor)
# end

# #w = FullSensNCPWiring(10,3)
# #heatmap(w.matrices[:w_sens])


text:
function value(obj::UnivariateFunction, arg0::jdouble)
    return jcall(obj, "value", jdouble, (jdouble,), arg0)
end


text:
using GenomicTiles
using BioFeatures
using DataFrames

function get_sub_regions(start, stop)
    subregions = DataFrame( start = [start, start + 50], stop = [start+60, start + 100])
    return Features( subregions )
end

#
# GET SOME COORDINATES TO SERVE AS A BASE THAT WILL BE SPLIT INTO TILES
#
region_lengths = int( rand(100)*1000 )
starts = 1:1000:100000
stops = starts + region_lengths

# Read these into a BioFeature - FeatureRegion
regions = DataFrame(start = starts, stop= stops)
FeaturesToBeSplitIntoTiles = Features( regions )

tiles = nothing
tiles_matrix = nothing
for features = values( FeaturesToBeSplitIntoTiles.features )
    for feature in values(features)
         subregions = get_sub_regions( feature[:start][1],feature[:stop][1] )
         tiles = make_tiles(feature, num_tiles=30, default_value = 1)
         #println(tiles)
         for tile in tiles
            tile = calculate_tile_value( tile, subregions, 0 )
         end
         tiles_matrix = add_tiles_to_tiles_matrix(tiles_matrix, tiles)
    end
end

(tiles_matrix_clustered,cluster_assignments) = cluster_tiles(tiles_matrix)
path = "/home/nic/Desktop/test.png"
plot_matrix2(tiles_matrix_clustered,path)
path = "/home/nic/Desktop/test-bar.png"
plot_cluster_bar_chart(tiles_matrix,cluster_assignments,path)

cluster_assignments
tiles_matrix'[1,:]
tiles_matrix'[100,:]
cluster_assignments
#clustered_tiles = cluster_tiles_matrix(tiles_matrix; method="kmeans",num_clusters=5)
#clustered_tiles is a tuple that includes the categories/assignments
#plot_cluster_tiles(clustered_tiles, path)
#plot_cluster_barplot(clustered_tiles, path)

#matrix = nothing
#tiles
# This first example considers the case where you havese several regions
# that overlap tiles. If the region overlaps the tiles then that tile should
# be counted as an overlap.

# The regions might have categories - but cannot overlap each other.
#  e.g. -----111---222-11-21-11

# for region in FeaturesToBeSplitIntoTiles
#   * generate some sub regions that will be tiled
#       subregions = get_sub_region(region.start,region.stop)
#   * get an array of tiles for the region to be tiled
#      tiles = make_tile(region)
#   for tile in tiles
#    *see which tiles the sub regions overlap and count, value_assignment is a vector for the values that should be assigned,
#    * if this is nothing - then just a count is made should an overlap be registered
#    * this sets a property in the tile
#       set_tile(tile,subregions,value_assignment)
# end

## Then put tiles into a matrix
#   matrix = nothing
#   foreach region ..
#      matrix_row = map(x -> x.value, tiles )
#      if matrix = nothing
#         matrix = matrix_row
#      else
#         hcat(matrix,matrix_row)
#      end
#

### Then
#
#   cluster_tiles = kmeans_tile_cluster(matrix,num_clusters)
#   plot_cluster_tile(cluster)
#   plot_cluster_barplot(cluster)
#
#   or
#
#   cluster_tiles = aff_prop_cluster(matrix, num_clusters)
#   ....

#=
starts_a = zeros(Int64,100)
stops_a  = zeros(Int64,100)
for i = 1:50
   starts_a[i] = starts[i]
   stops_a[i]  = starts[i] + 100
end

for i = 51:100
    starts_a[i] = stops[i] - 50
    stops_a[i] =  stops[i]
end

regions_a = DataFrame( start = starts_a, stop = stops_a)
FeaturesToBeCountedOverTiles = FeatureRegions(regions_a)

=#


text:
module OCToolbox

include("optimisers.jl")
include("tools.jl")
include("cost_functions.jl")
include("problems.jl")
include("integrators.jl")

export GRAPE, GRAPE_no_reshape
export export_pulse
export C1, C2, C3, C4, C5, C6, C7
export ProblemInfo
export pw_integrate, pw_full_evolution, ode_integrate, me_integrate

end


text:
@testset "Construction" for coeff in [
    Int64[1, 1, 1, 1],
    Float32[1, -4, 2],
    ComplexF64[1 - 1im, 2 + 3im],
    [3 // 4, -2 // 1, 1 // 1]
                                      ]
    p = Bernstein(coeff)
    @test p.coeffs == coeff
    @test coeffs(p) == coeff
    @test degree(p) <= length(coeff) - 1
    @test p.var == :x
    @test length(p) == length(coeff)
    @test size(p) == size(coeff)
    @test size(p, 1) == size(coeff, 1)
    @test typeof(p).parameters[end] == eltype(coeff)
    @test eltype(p) == eltype(coeff)
end

@testset "Other Construction" begin
    # Leading 0s are  trimmed
    p = Bernstein([1, 2, 0, 0])
    @test p.coeffs == [1, 2]

    # Different type
    p = Bernstein{3, Float64}(ones(Int32, 4))
    @test p.coeffs == ones(Float64, 4)

    p = Bernstein{1, Int}(30)
    @test p.coeffs == [30, 30]

    p = zero(Bernstein{1, Int})
    @test p.coeffs == Int[]

    p = zero(Bernstein{5, Int})
    @test length(p.coeffs)  ==  0

    p = one(Bernstein{1, Int})
    @test p.coeffs == ones(Int,  2)

    @test iszero(Bernstein{2, Int}(Int[]))

    p0 = Bernstein([0])
    @test iszero(p0)
    @test degree(p0) == -1
end

@testset "Roots" begin
    r = [1/4, 1/2, 3/4]
    c = fromroots(Bernstein, r)
    @test roots(c) ≈ r

    r = [1im, -1im]
    c = fromroots(Bernstein, r)
    @test roots(c) ≈ r
end

@testset "Values" begin
    c1 = Bernstein([2.5, 2.0, 1.5]) # 1 + 2x + 3x^2
    x =  [1/4, 1/2, 3/4]
    expected = [2.25, 2, 1.75]
    @test c1.(x) ≈ expected
end

@testset "Conversions" begin
    c1 = Bernstein([2.5, 2.0, 1.5])
    @test c1 ≈ Polynomial([2.5,  -1.0])
    p = convert(Polynomial{Float64}, c1)
    @test p == Polynomial{Float64}([2.5, -1.0])
    @test convert(Bernstein{2, Float64}, p) == c1

end

@testset "Mapdomain" begin
    x = -30:20
    mx = mapdomain(Bernstein, x)
    @test extrema(mx) == (0, 1)

    x = 0.5:0.01:0.6
    mx = mapdomain(Bernstein, x)
    @test extrema(mx) == (0, 1)
end

@testset "Arithmetic" begin
    # multiplication
    c1 = Bernstein([1, 2, 3])
    c2 = Bernstein([3, 2, 1])
    @test convert(Polynomial, c1 * c2) ≈ convert(Polynomial, c1) * convert(Polynomial,c2)

    c1 = Bernstein([1, 2, 3])
    c2 = Bernstein([1, 2, 3, 4])
    p1, p2 = convert(Polynomial,c1), convert(Polynomial,c2)
    @test p1*p2 ≈ convert(Polynomial, c1 * c2)

    c1 = Bernstein([1, 2, 3])
    c2 = Bernstein([3, 2, 1])
    d, r = divrem(c1, c2)
    @test d.coeffs ≈ [-1.0]
    @test r.coeffs ≈ [4.0]
    @test d * c2 + r ≈ c1

    c2 = Bernstein([0, 1, 2, 3])
    d, r = divrem(c2, c1)

    @test d.coeffs ≈ [1.5]
    @test r.coeffs ≈ [-1.5]
    z = d *  c1  +  r - c2
    @test coeffs(z)  ==  zeros(length(z))

    # GCD
    c1 = Bernstein(Float64[1//1, 2, 3])
    c2 = Bernstein(Float64[3//1, 2, 1])
    @test gcd(c1, c2) ≈ Bernstein([4.0])

    for i in 0:5, j in 0:5
        # multiplication
        target = zeros(i + j + 1)
        target[end] += 1
        c1 = Bernstein(vcat(zeros(i), 1))
        c2 = Bernstein(vcat(zeros(j), 1))
        @test c1 * c2 ≈ Bernstein(target)

        # divrem
        target = c1 + c2
        quo, rem = divrem(target, c1)
        res = quo * c1 + rem
        @test res ≈ target
    end
end

@testset "integrals derivatives" begin
    c1 = one(Bernstein{1,Int})
    @test integrate(c1) == Bernstein([0, 0.5, 1])
    for k in [-3, 0, 2]
        @test integrate(c1, k) == Bernstein(k .+ [0, 0.5, 1])
    end

    for i in 0:4
        scl = i + 1
        p = Polynomial(vcat(zeros(i), 1))
        pint = integrate(p, i)
        cheb = convert(Bernstein, p)
        cint = integrate(cheb, i)
        res = convert(Polynomial, cint)
        @test res ≈ pint
        @test derivative(cint) ≈ cheb
    end

    for i in 1:10
        p = Bernstein{5, Float64}(rand(1:5, 6))
        pp = p - (integrate ∘ derivative)(p)
        qq = p - (derivative ∘ integrate)(p)
        @test degree(truncate(pp, atol=1e-13)) <= 0
        @test degree(truncate(qq, atol=1e-13)) <  0
    end

end


text:
module RocketLazyObservableTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "LazyObservable" begin

    println("Testing: lazy")

    @testset begin

        values = []
        stream = lazy(Int)
        subscription = subscribe!(stream, (x) -> push!(values, x))
        set!(stream, from(1:3))
        @test values == [ 1, 2, 3 ]
        subscription = subscribe!(stream, (x) -> push!(values, x))
        @test values == [ 1, 2, 3, 1, 2, 3 ]

    end

    @testset begin

        values = []
        stream = lazy(Int)
        subscription = subscribe!(stream, (x) -> push!(values, x))
        unsubscribe!(subscription)
        set!(stream, from(1:3))
        @test values == [ ]
        subscription = subscribe!(stream, (x) -> push!(values, x))
        @test values == [ 1, 2, 3 ]

    end

    @testset begin

        values1 = []
        values2 = []
        stream = lazy(Int)
        subscription1 = subscribe!(stream, (x) -> push!(values1, x))
        subscription2 = subscribe!(stream, (x) -> push!(values2, x))
        set!(stream, from(1:3))
        @test values1 == [ 1, 2, 3 ]
        @test values2 == [ 1, 2, 3 ]

    end

    @testset begin

        values1 = []
        values2 = []
        stream = lazy(Int)
        subscription1 = subscribe!(stream, (x) -> push!(values1, x))
        subscription2 = subscribe!(stream, (x) -> push!(values2, x))
        unsubscribe!(subscription1)
        set!(stream, from(1:3))
        @test values1 == [ ]
        @test values2 == [ 1, 2, 3 ]

    end

    @testset begin

        values1 = []
        values2 = []
        stream = lazy(Int)
        subscription1 = subscribe!(stream, (x) -> push!(values1, x))
        subscription2 = subscribe!(stream, (x) -> push!(values2, x))
        unsubscribe!(subscription2)
        set!(stream, from(1:3))
        @test values1 == [ 1, 2, 3]
        @test values2 == [ ]

    end

    @testset begin

        values = []
        stream = lazy(Int)
        subject = Subject(Int)
        next!(subject, 1)
        subscription = subscribe!(stream, (x) -> push!(values, x))
        set!(stream, subject)
        @test values == [ ]
        next!(subject, 2)
        @test values == [ 2 ]
        unsubscribe!(subscription)
        next!(subject, 3)
        @test values == [ 2 ]

    end

end

end


text:
function mirror_image(k, b::Barrier, grd::Grid)

    _, proj_js, _ = findnz(grd.Proj)

    x, y = grd.points[:,k]
    k_cart = proj_js[k]
    j = Int(ceil(k_cart / grd._Nx))
    i = k_cart - (j - 1) * grd._Nx
    dist = 0
    maxsteps = 0

    while b.orientation*(b.func(x,y) - b.val) < 0
        fs = hcat([[b.orientation*(b.func(x+grd.dx*ii, y+grd.dy*jj) - b.val) for ii=-1:1] for jj=-1:1]...)
        di, dj = Tuple(findmax(fs)[2]) .- 2
        x += di*grd.dx
        y += dj*grd.dy
        i += di
        j += dj
        dist += sqrt(di^2 + dj^2)
        maxsteps += di^2 + dj^2
    end

    steps = 0
    while (dist > 0) && (steps < maxsteps)
        fs = hcat([[b.orientation*(b.func(x+grd.dx*ii, y+grd.dy*jj) - b.val) for ii=-1:1] for jj=-1:1]...)
        di, dj = Tuple(findmax(fs)[2]) .- 2
        x += di*grd.dx
        y += dj*grd.dy
        i += di
        j += dj
        dist -= sqrt(di^2 + dj^2)
        steps += 1
    end

    k_cart_new = i + (j-1)*grd._Nx

    return k_cart_new
end

struct GhostConditions

    Proj :: SparseMatrixCSC
    Mirror :: SparseMatrixCSC
    swaps :: Matrix{Int64}

end

function GhostConditions(bars::Vector{Barrier}, grd::Grid)

    chnl = RemoteChannel(()->Channel{Bool}(), 1)
    p = Progress(grd.Nk, desc="Finding mirror images... ")
    @async while take!(chnl)
        next!(p)
    end

    M = DArray((grd.Nk,length(bars)+2), workers(), (length(workers()), 1)) do inds

        ks = Int64[inds[1]...]
        swaps = ones(Int64, (length(inds[1]), length(bars)))
        j_mir = Int64[]
        j_proj = Int64[]
        _, proj_js, _ = findnz(grd.Proj)

        for i=1:length(ks)
            x, y = grd.points[:,ks[i]]

            flags = Bool[b.orientation*(b.func(x,y)-b.val) < 0 for b in bars]
            trueflag = findfirst(x->x, flags)
            if isnothing(trueflag)
                j_mir = [j_mir; proj_js[ks[i]]]
                j_proj = [j_proj; ks[i]]
            else
                swaps[i,trueflag] = -1
                j_mir = [j_mir; mirror_image(ks[i], bars[trueflag], grd)]
                j_proj = [j_proj; 0]
            end
            put!(chnl, true)

        end

        hcat(j_mir, j_proj, swaps)
    end

    j_mir = Array(M[:,1])
    j_proj = Array(M[:,2])
    j_proj = j_proj[j_proj .!= 0]
    swaps = Array(M[:,3:end])

    Mirror = sparse([1:grd.Nk...], j_mir, ones(grd.Nk), grd.Nk, grd._Nx*grd._Ny) * transpose(grd.Proj)
    Proj = sparse([1:length(j_proj)...], j_proj, ones(length(j_proj)), length(j_proj), grd.Nk)

    put!(chnl, false)
    return GhostConditions(Proj, Mirror, swaps)
end

function swap_sign(gc::GhostConditions, inds...)

    Mirror = gc.Mirror

    for i in inds
        Mirror = Diagonal(gc.swaps[:,i]) * Mirror
    end

    return GhostConditions(gc.Proj, Mirror, gc.swaps)
end

function extrapolate_ghosts(x::Vector{Float64}, xb::Vector{Float64}, gc::GhostConditions)

    return gc.Mirror*transpose(gc.Proj)*x + (I - gc.Mirror)*xb
end

function require_boundary_conditions(A::SparseMatrixCSC, gc::GhostConditions)

    P = gc.Proj
    Pt = transpose(P)

    Anew = P*A*gc.Mirror*Pt
    Bterm = P*A*(I - gc.Mirror)

    return Anew, Bterm
end


text:
module Diagnostics

export
    HorizontalAverage, TimeSeries, FieldMaximum,
    CFL, AdvectiveCFL, DiffusiveCFL, NaNChecker,
    run_diagnostic

using Oceananigans,
      Oceananigans.Operators

using Oceananigans: AbstractDiagnostic

include("diagnostics_kernels.jl")
include("nan_checker.jl")
include("horizontal_average.jl")
include("time_series.jl")
include("field_maximum.jl")
include("cfl.jl")

end


text:
using Bonsai, FilePaths
using FilePathsBase: /
using HTTP

folder = Path(@__DIR__) / "data"
app = App()
app.hot_reload = true

app.get("/") do stream
    Bonsai.write(stream, folder / "index.html")
    HTTP.setheader(stream, "Content-Type" => "text/plain")
end

app.get("**") do stream, next
    try
        next(stream)
    catch e
        @error e
        Bonsai.write(stream, Body(error=repr(e)))
    end
end

start(app, port=10001)
wait(app)
stop(app)

max_stars_repo_name: onetonfoot/Tree

max_issues_repo_name: onetonfoot/Tree

max_forks_repo_name: onetonfoot/Tree


text:
using BasisFunctions, FrameFun.ExtensionFrames, FrameFun.ExtensionFramePlatforms,
    DomainSets, Test, StaticArrays, FrameFun.Platforms, FrameFun.FrameFunInterface

@testset "ExtensionFramePlatforms" begin
    d1 = 0.0..0.5
    d2 = (0.0..0.5)^2
    d3 = disk(.4,SVector(.5,.5))
    dict1 = extensionframe(Fourier(10),d1)
    dict2 = extensionframe(Fourier(10)^2,d2)
    dict3 = extensionframe(Fourier(10)^2,d3)
    plat1 = (platform(dict1),10)
    plat2 = (platform(dict2),(10,10))
    plat3 = (platform(dict3),(10,10))

    @test dictionary(plat1...) isa ExtensionFrame
    @test dictionary(plat2...) isa ExtensionFrameTensor
    @test dictionary(plat3...) isa ExtensionFrame

    @test oversampling_grid(plat1...;oversamplingfactor=2)≈subgrid(FourierGrid(38),d1)
    @test oversampling_grid(plat2...;oversamplingfactor=2)≈subgrid(FourierGrid(27)^2,d2)
    @test oversampling_grid(plat3...;oversamplingfactor=2)≈subgrid(FourierGrid(21)^2,d3)
    @test length(oversampling_grid(plat1...;oversamplingfactor=2)) == 20
    @test abs(length(oversampling_grid(plat2...;oversamplingfactor=2)) - 200)<20
    @test abs(length(oversampling_grid(plat3...;oversamplingfactor=2)) - 200 )< 20

    @test SamplingStyle(platform(dict1)) == OversamplingStyle()
    @test SamplingStyle(platform(dict2)) == ProductSamplingStyle(OversamplingStyle(),OversamplingStyle())
    @test SamplingStyle(platform(dict3)) == OversamplingStyle()

    @test SolverStyle(platform(dict1),OversamplingStyle()) == AZStyle()
    @test SolverStyle(platform(dict2),ProductSamplingStyle(OversamplingStyle(),OversamplingStyle())) == ProductSolverStyle(AZStyle(),AZStyle())
    @test SolverStyle(platform(dict3),OversamplingStyle()) == AZStyle()

    @test DictionaryStyle(platform(dict1)) == FrameStyle()
    @test DictionaryStyle(platform(dict2)) == FrameStyle()
    @test DictionaryStyle(platform(dict3)) == FrameStyle()

    ddict1 = dualdictionary(plat1...,measure(plat1[1]))
    @test ddict1 isa ExtensionFrame
    @test operator(basis(ddict1)) ≈ operator(dual(Fourier(10),FourierMeasure()))
    ddict2 = dualdictionary(plat2...,measure(plat2[1]))
    @test ddict2 isa ExtensionFrameTensor
    ddict3 = dualdictionary(plat3...,measure(plat3[1]))
    @test ddict3 isa ExtensionFrame
end


text:
using BenchmarkTools
using Meshing
using GeometryTypes

# Define a parent BenchmarkGroup to contain our suite
const suite = BenchmarkGroup()
suite["SDF Construction"] = BenchmarkGroup()
suite["SDF Mesh"] = BenchmarkGroup()
suite["Function Mesh"] = BenchmarkGroup()

println("Benchmarking Meshing.jl...")

#
# Algorithms to benchmark
#

algos_sdf = [MarchingCubes(), MarchingTetrahedra(), NaiveSurfaceNets()]
algos_fn = [MarchingCubes(), NaiveSurfaceNets()]

#
# Benchmark SDF constructon for SDF and Direct sample comparisons
#

suite["SDF Construction"]["Torus"] = @benchmarkable SignedDistanceField(HyperRectangle(Vec(-2,-2,-2.),Vec(4,4,4.)),0.05) do v
    (sqrt(v[1]^2+v[2]^2)-0.5)^2 + v[3]^2 - 0.25 # torus
end

#
# Torus Constructions
#

sdf_torus = SignedDistanceField(HyperRectangle(Vec(-2,-2,-2.),Vec(4,4,4.)),0.05) do v
    (sqrt(v[1]^2+v[2]^2)-0.5)^2 + v[3]^2 - 0.25 # torus
end

fn_torus(v) = (sqrt(v[1]^2+v[2]^2)-0.5)^2 + v[3]^2 - 0.25 # torus

#
# Benchmark algorithms
#

for algo in algos_sdf
    suite["SDF Mesh"][string(typeof(algo))] = @benchmarkable HomogenousMesh(sdf_torus, $algo)
end

for algo in algos_fn
    suite["Function Mesh"][string(typeof(algo))] = @benchmarkable HomogenousMesh(fn_torus, HyperRectangle(Vec(-2,-2,-2.),Vec(4,4,4.)), (81,81,81), $algo)
end

# If a cache of tuned parameters already exists, use it, otherwise, tune and cache
# the benchmark parameters. Reusing cached parameters is faster and more reliable
# than re-tuning `suite` every time the file is included.
paramspath = joinpath(dirname(@__FILE__), "params.json")

if isfile(paramspath)
    loadparams!(suite, BenchmarkTools.load(paramspath)[1], :evals);
else
    tune!(suite)
    BenchmarkTools.save(paramspath, params(suite));
end

#
# Perform benchmarks and print results
#

results = run(suite)

for trial in results
    ctx = IOContext(stdout, :verbose => true, :compact => false)
    println(ctx)
    println(ctx, trial.first)
    println(ctx, trial.second)
end


text:
@inline function initialize!(integrator,cache::BS3ConstantCache,f=integrator.f)
  integrator.kshortsize = 2
  integrator.k = eltype(integrator.sol.k)(integrator.kshortsize)
  integrator.fsalfirst = f(integrator.t,integrator.uprev) # Pre-start fsal
end

@inline function perform_step!(integrator,cache::BS3ConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  @unpack a21,a32,a41,a42,a43,c1,c2,b1,b2,b3,b4 = cache
  k1 = integrator.fsalfirst
  a1 = dt*a21
  k2 = f(@muladd(t+c1*dt),@muladd(uprev+a1*k1))
  a2 = dt*a32
  k3 = f(@muladd(t+c2*dt),@muladd(uprev+a2*k2))
  u = @muladd uprev+dt*(a41*k1+a42*k2+a43*k3)
  k4 = f(t+dt,u); integrator.fsallast = k4
  if integrator.opts.adaptive
    utilde = @muladd uprev + dt*(b1*k1 + b2*k2 + b3*k3 + b4*k4)
    integrator.EEst = integrator.opts.internalnorm( ((utilde-u)./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol)))
  end
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  @pack integrator = t,dt,u
end

@inline function initialize!(integrator,cache::BS3Cache,f=integrator.f)
  integrator.kshortsize = 2
  integrator.k = eltype(integrator.sol.k)(integrator.kshortsize)
  integrator.fsalfirst = cache.fsalfirst  # done by pointers, no copying
  integrator.fsallast = cache.k4
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  f(integrator.t,integrator.uprev,integrator.fsalfirst) # Pre-start fsal
end

@inline function perform_step!(integrator,cache::BS3Cache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  uidx = eachindex(integrator.uprev)
  @unpack k2,k3,k4,utilde,tmp,atmp = cache
  k1 = cache.fsalfirst
  @unpack a21,a32,a41,a42,a43,c1,c2,b1,b2,b3,b4 = cache.tab
  a1 = dt*a21
  for i in uidx
    tmp[i] = @muladd uprev[i]+a1*k1[i]
  end
  f(@muladd(t+c1*dt),tmp,k2)
  a2 = dt*a32
  for i in uidx
    tmp[i] = @muladd uprev[i]+a2*k2[i]
  end
  f(@muladd(t+c2*dt),tmp,k3)
  for i in uidx
    u[i] = @muladd uprev[i]+dt*(a41*k1[i]+a42*k2[i]+a43*k3[i])
  end
  f(t+dt,u,k4)
  if integrator.opts.adaptive
    for i in uidx
      utilde[i] = @muladd uprev[i] + dt*(b1*k1[i] + b2*k2[i] + b3*k3[i] + b4*k4[i])
      atmp[i] = ((utilde[i]-u[i])./@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol))
    end
    integrator.EEst = integrator.opts.internalnorm(atmp)
  end
  @pack integrator = t,dt,u,k
end

@inline function initialize!(integrator,cache::BS5ConstantCache,f=integrator.f)
  integrator.kshortsize = 8
  integrator.k = eltype(integrator.sol.k)(integrator.kshortsize)
  integrator.fsalfirst = f(integrator.t,integrator.uprev) # Pre-start fsal
end

@inline function perform_step!(integrator,cache::BS5ConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
        @unpack c1,c2,c3,c4,c5,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,a81,a83,a84,a85,a86,a87,bhat1,bhat3,bhat4,bhat5,bhat6,btilde1,btilde2,btilde3,btilde4,btilde5,btilde6,btilde7,btilde8 = cache
  k1 = integrator.fsalfirst
  a = dt*a21
  k2 = f(@muladd(t+c1*dt),@muladd(uprev+a*k1))
  k3 = f(@muladd(t+c2*dt),@muladd(uprev+dt*(a31*k1+a32*k2)))
  k4 = f(@muladd(t+c3*dt),@muladd(uprev+dt*(a41*k1+a42*k2+a43*k3)))
  k5 = f(@muladd(t+c4*dt),@muladd(uprev+dt*(a51*k1+a52*k2+a53*k3+a54*k4)))
  k6 = f(@muladd(t+c5*dt),@muladd(uprev+dt*(a61*k1+a62*k2+a63*k3+a64*k4+a65*k5)))
  k7 = f(t+dt,@muladd(uprev+dt*(a71*k1+a72*k2+a73*k3+a74*k4+a75*k5+a76*k6)))
  u = @muladd uprev+dt*(a81*k1+a83*k3+a84*k4+a85*k5+a86*k6+a87*k7)
  integrator.fsallast = f(t+dt,u); k8 = integrator.fsallast
  if integrator.opts.adaptive
    uhat   = dt*@muladd(bhat1*k1 + bhat3*k3 + bhat4*k4 + bhat5*k5 + bhat6*k6)
    utilde = @muladd uprev + dt*(btilde1*k1 + btilde2*k2 + btilde3*k3 + btilde4*k4 + btilde5*k5 + btilde6*k6 + btilde7*k7 + btilde8*k8)
    EEst1 = integrator.opts.internalnorm( sum(((uhat)./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol))))
    EEst2 = integrator.opts.internalnorm( sum(((utilde-u)./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol))))
    integrator.EEst = max(EEst1,EEst2)
  end
  integrator.k[1]=k1; integrator.k[2]=k2; integrator.k[3]=k3;integrator.k[4]=k4;integrator.k[5]=k5;integrator.k[6]=k6;integrator.k[7]=k7;integrator.k[8]=k8
  @pack integrator = t,dt,u,k
end

@inline function initialize!(integrator,cache::BS5Cache,f=integrator.f)
  integrator.kshortsize = 8
  integrator.k = eltype(integrator.sol.k)(integrator.kshortsize)
  integrator.k = eltype(integrator.sol.k)(integrator.kshortsize)
  integrator.k[1]=cache.k1; integrator.k[2]=cache.k2;
  integrator.k[3]=cache.k3; integrator.k[4]=cache.k4;
  integrator.k[5]=cache.k5; integrator.k[6]=cache.k6;
  integrator.k[7]=cache.k7; integrator.k[8]=cache.k8
  integrator.fsalfirst = cache.k1; integrator.fsallast = cache.k8  # setup pointers
  f(integrator.t,integrator.uprev,integrator.fsalfirst) # Pre-start fsal
end

@inline function perform_step!(integrator,cache::BS5Cache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  uidx = eachindex(integrator.uprev)
  @unpack k1,k2,k3,k4,k5,k6,k7,k8,utilde,uhat,tmp,atmp,atmptilde = cache
  @unpack c1,c2,c3,c4,c5,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,a81,a83,a84,a85,a86,a87,bhat1,bhat3,bhat4,bhat5,bhat6,btilde1,btilde2,btilde3,btilde4,btilde5,btilde6,btilde7,btilde8 = cache.tab
  a = dt*a21
  for i in uidx
    tmp[i] = @muladd uprev[i]+a*k1[i]
  end
  f(@muladd(t+c1*dt),tmp,k2)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a31*k1[i]+a32*k2[i])
  end
  f(@muladd(t+c2*dt),tmp,k3)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a41*k1[i]+a42*k2[i]+a43*k3[i])
  end
  f(@muladd(t+c3*dt),tmp,k4)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a51*k1[i]+a52*k2[i]+a53*k3[i]+a54*k4[i])
  end
  f(@muladd(t+c4*dt),tmp,k5)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a61*k1[i]+a62*k2[i]+a63*k3[i]+a64*k4[i]+a65*k5[i])
  end
  f(@muladd(t+c5*dt),tmp,k6)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a71*k1[i]+a72*k2[i]+a73*k3[i]+a74*k4[i]+a75*k5[i]+a76*k6[i])
  end
  f(t+dt,tmp,k7)
  for i in uidx
    u[i] = @muladd uprev[i]+dt*(a81*k1[i]+a83*k3[i]+a84*k4[i]+a85*k5[i]+a86*k6[i]+a87*k7[i])
  end
  f(t+dt,u,k8)
  if integrator.opts.adaptive
    for i in uidx
      uhat[i]   = dt*@muladd(bhat1*k1[i] + bhat3*k3[i] + bhat4*k4[i] + bhat5*k5[i] + bhat6*k6[i])
      utilde[i] = @muladd uprev[i] + dt*(btilde1*k1[i] + btilde2*k2[i] + btilde3*k3[i] + btilde4*k4[i] + btilde5*k5[i] + btilde6*k6[i] + btilde7*k7[i] + btilde8*k8[i])
      atmp[i] = ((uhat[i])./@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol))
      atmptilde[i] = ((utilde[i]-u[i])./@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol))
    end
    EEst1 = integrator.opts.internalnorm(atmp)
    EEst2 = integrator.opts.internalnorm(atmptilde)
    integrator.EEst = max(EEst1,EEst2)
  end
  @pack integrator = t,dt,u,k
end

@inline function initialize!(integrator,cache::Tsit5ConstantCache,f=integrator.f)
  integrator.kshortsize = 7
  integrator.k = eltype(integrator.sol.k)(integrator.kshortsize)
  integrator.fsalfirst = f(integrator.t,integrator.uprev) # Pre-start fsal
end

@inline function perform_step!(integrator,cache::Tsit5ConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  @unpack c1,c2,c3,c4,c5,c6,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,b1,b2,b3,b4,b5,b6,b7 = cache
  k1 = integrator.fsalfirst
  a = dt*a21
  k2 = f(@muladd(t+c1*dt),@muladd(uprev+a*k1))
  k3 = f(@muladd(t+c2*dt),@muladd(uprev+dt*(a31*k1+a32*k2)))
  k4 = f(@muladd(t+c3*dt),@muladd(uprev+dt*(a41*k1+a42*k2+a43*k3)))
  k5 = f(@muladd(t+c4*dt),@muladd(uprev+dt*(a51*k1+a52*k2+a53*k3+a54*k4)))
  k6 = f(t+dt,@muladd(uprev+dt*(a61*k1+a62*k2+a63*k3+a64*k4+a65*k5)))
  u = @muladd uprev+dt*(a71*k1+a72*k2+a73*k3+a74*k4+a75*k5+a76*k6)
  integrator.fsallast = f(t+dt,u); k7 = integrator.fsallast
  if integrator.opts.adaptive
    utilde = @muladd uprev + dt*(b1*k1 + b2*k2 + b3*k3 + b4*k4 + b5*k5 + b6*k6 + b7*k7)
    integrator.EEst = integrator.opts.internalnorm(((utilde-u)./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol)))
  end
  integrator.k[1] = k1
  integrator.k[2] = k2
  integrator.k[3] = k3
  integrator.k[4] = k4
  integrator.k[5] = k5
  integrator.k[6] = k6
  integrator.k[7] = k7
  @pack integrator = t,dt,u,k
end

@inline function initialize!(integrator,cache::Tsit5Cache,f=integrator.f)
  integrator.kshortsize = 7
  integrator.fsalfirst = cache.k1; integrator.fsallast = cache.k7 # setup pointers
  integrator.k = eltype(integrator.sol.k)(integrator.kshortsize)
  # Setup k pointers
  integrator.k[1] = cache.k1
  integrator.k[2] = cache.k2
  integrator.k[3] = cache.k3
  integrator.k[4] = cache.k4
  integrator.k[5] = cache.k5
  integrator.k[6] = cache.k6
  integrator.k[7] = cache.k7
  f(integrator.t,integrator.uprev,integrator.fsalfirst) # Pre-start fsal
end

@inline function perform_step!(integrator,cache::Tsit5Cache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  uidx = eachindex(integrator.uprev)
  @unpack c1,c2,c3,c4,c5,c6,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,b1,b2,b3,b4,b5,b6,b7 = cache.tab
  @unpack k1,k2,k3,k4,k5,k6,k7,utilde,tmp,atmp = cache
  a = dt*a21
  for i in uidx
    tmp[i] = @muladd uprev[i]+a*k1[i]
  end
  f(@muladd(t+c1*dt),tmp,k2)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a31*k1[i]+a32*k2[i])
  end
  f(@muladd(t+c2*dt),tmp,k3)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a41*k1[i]+a42*k2[i]+a43*k3[i])
  end
  f(@muladd(t+c3*dt),tmp,k4)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a51*k1[i]+a52*k2[i]+a53*k3[i]+a54*k4[i])
  end
  f(@muladd(t+c4*dt),tmp,k5)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a61*k1[i]+a62*k2[i]+a63*k3[i]+a64*k4[i]+a65*k5[i])
  end
  f(t+dt,tmp,k6)
  for i in uidx
    u[i] = @muladd uprev[i]+dt*(a71*k1[i]+a72*k2[i]+a73*k3[i]+a74*k4[i]+a75*k5[i]+a76*k6[i])
  end
  f(t+dt,u,k7)
  if integrator.opts.adaptive
    for i in uidx
      utilde[i] = @muladd uprev[i] + dt*(b1*k1[i] + b2*k2[i] + b3*k3[i] + b4*k4[i] + b5*k5[i] + b6*k6[i] + b7*k7[i])
      atmp[i] = ((utilde[i]-u[i])./@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol))
    end
    integrator.EEst = integrator.opts.internalnorm(atmp)
  end
  @pack integrator = t,dt,u,k
end

@inline function initialize!(integrator,cache::DP5ConstantCache,f=integrator.f)
  integrator.kshortsize = 4
  integrator.k = eltype(integrator.sol.k)(integrator.kshortsize)
  integrator.fsalfirst = f(integrator.t,integrator.uprev) # Pre-start fsal
end

@inline function perform_step!(integrator,cache::DP5ConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  @unpack a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a73,a74,a75,a76,b1,b3,b4,b5,b6,b7,c1,c2,c3,c4,c5,c6 = cache
  @unpack d1,d3,d4,d5,d6,d7 = cache
  k1 = integrator.fsalfirst
  a = dt*a21
  k2 = f(@muladd(t+c1*dt),@muladd(uprev+a*k1))
  k3 = f(@muladd(t+c2*dt),@muladd(uprev+dt*(a31*k1+a32*k2)))
  k4 = f(@muladd(t+c3*dt),@muladd(uprev+dt*(a41*k1+a42*k2+a43*k3)))
  k5 = f(@muladd(t+c4*dt),@muladd(uprev+dt*(a51*k1+a52*k2+a53*k3+a54*k4)))
  k6 = f(t+dt,@muladd uprev+dt*(a61*k1+a62*k2+a63*k3+a64*k4+a65*k5))
  update = @muladd a71*k1+a73*k3+a74*k4+a75*k5+a76*k6
  u = @muladd uprev+dt*update
  integrator.fsallast = f(t+dt,u); k7 = integrator.fsallast
  if integrator.opts.adaptive
    utilde = uprev + dt*(b1*k1 + b3*k3 + b4*k4 + b5*k5 + b6*k6 + b7*k7)
    integrator.EEst = integrator.opts.internalnorm( ((utilde-u)./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol)))
  end
  integrator.k[1] = update
  bspl = k1 - update
  integrator.k[2] = bspl
  integrator.k[3] = update - k7 - bspl
  integrator.k[4] = @muladd (d1*k1+d3*k3+d4*k4+d5*k5+d6*k6+d7*k7)
  @pack integrator = t,dt,u,k
end

@inline function initialize!(integrator,cache::DP5Cache,f=integrator.f)
  integrator.kshortsize = 4
  integrator.k = [cache.update,cache.bspl,cache.dense_tmp3,cache.dense_tmp4]
  integrator.fsalfirst = cache.k1; integrator.fsallast = cache.k7
  f(integrator.t,integrator.uprev,integrator.fsalfirst) # Pre-start fsal
end

@inline function perform_step!(integrator,cache::DP5Cache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  uidx = eachindex(integrator.uprev)
  @unpack a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a73,a74,a75,a76,b1,b3,b4,b5,b6,b7,c1,c2,c3,c4,c5,c6 = cache.tab
  @unpack k1,k2,k3,k4,k5,k6,k7,dense_tmp3,dense_tmp4,update,bspl,utilde,tmp,atmp = cache
  @unpack d1,d3,d4,d5,d6,d7 = cache.tab
  a = dt*a21
  for i in uidx
    tmp[i] = @muladd uprev[i]+a*k1[i]
  end
  f(@muladd(t+c1*dt),tmp,k2)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a31*k1[i]+a32*k2[i])
  end
  f(@muladd(t+c2*dt),tmp,k3)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a41*k1[i]+a42*k2[i]+a43*k3[i])
  end
  f(@muladd(t+c3*dt),tmp,k4)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a51*k1[i]+a52*k2[i]+a53*k3[i]+a54*k4[i])
  end
  f(@muladd(t+c4*dt),tmp,k5)
  for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a61*k1[i]+a62*k2[i]+a63*k3[i]+a64*k4[i]+a65*k5[i])
  end
  f(t+dt,tmp,k6)
  for i in uidx
    update[i] = @muladd a71*k1[i]+a73*k3[i]+a74*k4[i]+a75*k5[i]+a76*k6[i]
    u[i] = @muladd uprev[i]+dt*update[i]
  end
  f(t+dt,u,k7);
  if integrator.opts.adaptive
    for i in uidx
      utilde[i] = @muladd uprev[i] + dt*(b1*k1[i] + b3*k3[i] + b4*k4[i] + b5*k5[i] + b6*k6[i] + b7*k7[i])
      atmp[i] = ((utilde[i]-u[i])./@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol))
    end
    integrator.EEst = integrator.opts.internalnorm(atmp)
  end
  for i in uidx
    bspl[i] = k1[i] - update[i]
    integrator.k[3][i] = update[i] - k7[i] - bspl[i]
    integrator.k[4][i] = @muladd(d1*k1[i]+d3*k3[i]+d4*k4[i]+d5*k5[i]+d6*k6[i]+d7*k7[i])
  end
  @pack integrator = t,dt,u,k
end


text:
using Test

include("test_geocoding.jl")
include("test_geonames.jl")