text:
using DrWatson
@quickactivate :FractalDimension # uses DynamicalSystems, PyPlot
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, BenchmarkTools

import ChaosTools: float_to_int, molteno_boxing, _molteno_boxing, estimate_r0_buenoorovio, correlation_boxing

function float_to_int(data::Dataset{D,T}) where {D, T}
    N = length(data)
    mins, maxs = minmaxima(data)
    sizes = maxs .- mins
    ε0 = maximum(sizes)
    I = T(typemax(UInt64))
    # Let f:[min,max] -> [0+eps(T),1-eps(T)]*typemax(UInt64), then f(x) = m*x + b
    m = (1 - 2eps(T)) * (I - 2eps(I)) / ε0
    b = eps(I) .- mins .* m

    res = Vector{SVector{D,UInt64}}()
    sizehint!(res, N)
    for x in data
        int_val = floor.(UInt64, m .* x + b)
        push!(res, int_val)
    end
    Dataset(res), ε0
end

return_T_eps(X::Dataset{D,T}) where {D,T} = eps(T(typemax(UInt64)))

function molteno_boxing(data::Dataset, k0 = 10)
    integers, ε0 = float_to_int(data)
    boxes = _molteno_boxing(integers, k0)
    εs = ε0 ./ 2 .^ (1:length(boxes))
    return boxes, εs
end

function estimate_r0_buenoorovio(data)
    r0 = zero(eltype(data))
    mini, maxi = minmaxima(data)
    N = length(data)
    m = DelayEmbeddings.dimension(data)
    R = maximum(maxi .- mini)
    # The possibility of a bad pick exists, if so, the calculation is repeated.
    ν = zero(eltype(data))
    # Sample N/10 datapoints out of data for rough estimate of effective size.
    data_sample2 = data[unique(rand(1:N, N÷10))] |> Dataset
    r_ℓ = R / 10
    η_ℓ = length(correlation_boxing(data_sample2, r_ℓ)[1])
    r0 = zero(eltype(data))
    while true
        # Sample √N datapoints for rough dimension estimate
        data_sample1 = data[unique(rand(1:N, ceil(Int, sqrt(N))))] |> Dataset
        # Define logarithmic series of radii.
        lower = log10(min_pairwise_distance(data_sample1)[2])
        εs = 10 .^ range(lower, stop = log10(R), length = 16)
        # Estimate ν from a sample using the Grassberger Procaccia algorithm.
        cm = correlationsum(data_sample1, εs)
        ν = linear_region(log.(εs), log.(cm), tol = 0.5)[2]
        # Estimate the effictive size of the chaotic attractor.
        ℓ = r_ℓ * η_ℓ^(1/ν)
        # Calculate the optimal number of filled boxes according to Bueno-Orovio
        η_opt = N^(2/3) * ((3^ν - 1/2) / (3^m - 1))^(1/2)
        # The optimal box size is the effictive size divided by the box number # to the power of the inverse dimension.
        r0 = ℓ / η_opt^(1/ν)
        !isnan(r0) && break
    end
    r0
end

@btime()

# %% TODO: Also compare an embedded set for different `d` and `τ`.

# %% sensitivity to data dimensionality
N, dt = 2000, 0.1
fig, axs = subplots(4,1; sharex = true, figsize = (10,16))
for (j, D) in enumerate(4:3:16)
    color = "C$(j-1)"
    ds = Systems.lorenz96(D, range(0; length = D, step = 0.1); F = 8.0)
    X = trajectory(ds, N*dt; dt = dt, Ttr = 100.0)
    r0 = estimate_r0_buenoorovio(X)

    εs = 10.0 .^ range(-4 + D/4, 2, length = 32)
    bεs = r0 .* 10 .^ range(log10(min_pairwise_distance(X)[2]/r0), stop = 0, length = 16)
    m, mεs = molteno_boxing(X)

    H2 = genentropy.(2, εs, Ref(X); base = 10.0)
    H1 = genentropy.(2, m, base = 10.0)
    C1 = correlationsum(X, εs; w = 5)
    C2 = boxed_correlationsum(X, bεs, r0)
    axs[1].plot(log10.(εs), -H2, label = "\$ D = $D\$", color = color)
    if D == 4 || D == 7
        i = findfirst(z -> z > 0, 10 .^ (-H1))
        x, y = log10.(mεs)[i:end], -H1[i:end]
        is, d = linear_region(x, y)
        axs[2].plot(x, y, zorder = 1, color = color, label = "\$ d = $(round(d, digits = 2)) \$")
        axs[2].plot([x[is[1]], x[is[end]]], [y[is[1]], y[is[end]]],
        zorder = 2, ms = 5, marker = "o", ls = "None", color = color)
    else
        axs[2].plot(log10.(mεs), -H1, color = color)
    end
    axs[4].plot(log10.(bεs), log10.(C2))
    axs[4].set_xlabel("\$\\log_{10}(\\epsilon)\$")
    axs[4].set_ylabel("\$\\log_{10}(C_2)\$")
    axs[3].set_ylabel("\$\\log_{10}(C_1)\$")
    axs[2].set_ylabel("\$-H_1\$")
    axs[1].set_ylabel("\$-H_2\$")
    for (number, (set, es)) in enumerate(zip([C1, C2], [εs, bεs]))
        i = findfirst(z -> z > 0, set)
        x, y = log10.(es)[i:end], log10.(set)[i:end]
        is, d = linear_region(x, y)
        axs[2+number].plot(x, y, zorder = 1, color = color, label = "\$ d = $(round(d, digits = 2)) \$")
        axs[2+number].plot([x[is[1]], x[is[end]]], [y[is[1]], y[is[end]]],
        zorder = 2, ms = 5, marker = "o", ls = "None", color = color)
    end
end
for ax in axs; ax.legend(loc = "upper left"); end
fig.tight_layout()
axs[1].tick_params(labelbottom=false)
axs[1].set_title("increasing data dimension (Lorenz-96, N=$(N))")
wsave("plots/comparing_datadim.png", fig)
# comment: N = 100_000 took so much time on my pc I had to shut it down.

# %% Sensititivy to initial condition

N, dt, D = 5000, 0.1, 4
fig, axs = subplots(4,1; sharex = true, figsize = (10,16))
for (j, __D) in enumerate(4:3:16)
    color = "C$(j-1)"
    ds = Systems.lorenz96(D; F = 8.0)
    X = trajectory(ds, N*dt; dt = dt, Ttr = 10.0)

    r0 = estimate_r0_buenoorovio(X)

    εs = 10.0 .^ range(-4 + D/4, 2, length = 32)
    bεs = r0 .* 10 .^ range(log10(min_pairwise_distance(X)[2]/r0), stop = 0, length = 16)
    m, mεs = molteno_boxing(X)

    H2 = genentropy.(2, εs, Ref(X); base = 10.0)
    H1 = genentropy.(2, m, base = 10.0)
    C1 = correlationsum(X, εs; w = 5)
    C2 = boxed_correlationsum(X, bεs, r0)
    axs[1].plot(log10.(εs), -H2, label = nothing, color = color)
    axs[2].plot(log10.(mεs), -H1, color = color)
    axs[4].plot(log10.(bεs), log10.(C2))
    axs[4].set_xlabel("\$\\log_{10}(\\epsilon)\$")
    axs[4].set_ylabel("\$\\log_{10}(C_2)\$")
    axs[3].set_ylabel("\$\\log_{10}(C_1)\$")
    axs[2].set_ylabel("\$-H_1\$")
    axs[1].set_ylabel("\$-H_2\$")
    for (number, (set, es)) in enumerate(zip([10 .^ -H1, C1, C2], [mεs, εs, bεs]))
        i = findfirst(z -> z > 0, set)
        x, y = log10.(es)[i:end], log10.(set)[i:end]
        is, d = linear_region(x, y)
        axs[1+number].plot(x, y, zorder = 1, color = color, label = "\$ d = $(round(d, digits = 2))\$")
        axs[1+number].plot([x[is[1]], x[is[end]]], [y[is[1]], y[is[end]]],
        zorder = 2, ms = 5, marker = "o", ls = "None", color = color)
    end
end
for ax in axs; ax.legend(loc = "upper left"); ax.grid(); end
axs[1].tick_params(labelbottom=false)
axs[1].set_title("different initial conditions (Lorenz-96, N=$(N))")
fig.savefig("plots/comparing_initcon.png")

# %% Sensititivy to trajectory length

dt = 0.1
fig, axs = subplots(4,1; sharex = true, figsize = (10, 16))
for (j, N) in enumerate([500, 1000, 5000, 10000, 20000, 50000])
    color = "C$(j-1)"
    ds = Systems.lorenz()
    X = trajectory(ds, N*dt; dt = dt, Ttr = 10.0)

    r0 = estimate_r0_buenoorovio(X)

    εs = 10.0 .^ range(-4 + 3/4, 2, length = 40)
    bεs = r0 .* 10 .^ range(log10(min_pairwise_distance(X)[2]/r0), stop = 0, length = 40)
    m, mεs = molteno_boxing(X)

    H2 = genentropy.(2, εs, Ref(X); base = 10.0)
    H1 = genentropy.(2, m, base = 10.0)
    C1 = [correlationsum(X, ε; w = 5) for ε in εs]
    C2 = boxed_correlationsum(X, bεs, r0)
    axs[1].plot(log10.(εs), -H2, label = "\$N = $N\$", color = color)
    #axs[2].plot(log10.(mεs), -H1, color = color)
    #axs[4].plot(log10.(bεs), log10.(C2))
    axs[4].set_xlabel("\$\\log_{10}(\\epsilon)\$")
    axs[4].set_ylabel("\$\\log_{10}(C_2)\$")
    axs[3].set_ylabel("\$\\log_{10}(C_1)\$")
    axs[2].set_ylabel("\$-H_1\$")
    axs[1].set_ylabel("\$-H_2\$")
    for (number, (set, es)) in enumerate(zip([10 .^ -H1, C1, C2], [mεs, εs, bεs]))
        i = findfirst(z -> z > 0, set)
        x, y = log10.(es)[i:end], log10.(set)[i:end]
        is, d = linear_region(x, y)
        axs[1+number].plot(x, y, zorder = 1, color = color, label = "\$ d = $(round(d,digits=2))\$")
        axs[1+number].plot([x[is[1]], x[is[end]]], [y[is[1]], y[is[end]]],
        zorder = 2, ms = 5, marker = "o", ls = "None", color = color)
    end
end
for ax in axs; ax.legend(loc = "upper left"); ax.grid(); end
axs[1].tick_params(labelbottom=false)
axs[1].set_title("different trajectory lengths (Lorenz63)")
fig.savefig("plots/comparing_length.png")

# %% sensitivity to underlying chaotic set
N, dt = 10000, 0.1
D = 7

lo67 = Systems.lorenz()
X1 = trajectory(lo67, N*dt; dt = dt, Ttr = 10.0)

lo96 = Systems.lorenz96(D, range(0; length = D, step = 0.1); F = 8.0)
X2 = trajectory(lo96, N*dt; dt = dt, Ttr = 10.0)

he = Systems.henon()
X3 = trajectory(he, N; Ttr = 10)

to = Systems.towel()
X4 = trajectory(to, N; Ttr = 10)

sm = Systems.coupledstandardmaps(4; ks = fill(0.5, 4), Γ = 0.5)
X5 = trajectory(sm, N; Ttr = 100)

sm = Systems.coupledstandardmaps(4; ks = fill(0.5, 4), Γ = 0.5)
X5 = trajectory(sm, N; Ttr = 100)

using Distributions # for multivariate gaussian
dist = MvNormal(6, 0.1)
X6 = Dataset(rand(dist, N)')

names = ["Lorenz63", "Lorenz96 (7)", "Hénon", "towel", "coupled SM (8)", "gaussian (6)"]
lss = ["-", "--", ":", "-."]
fig, axs = subplots(4,1; sharex = true, figsize = (10,16))

for (j, X) in enumerate((X1, X2, X3, X4, X5, X6))
    color = "C$(j-1)"
    ls = lss[mod1(j, 4)]
    name = names[j]

    r0 = estimate_r0_buenoorovio(X)

    εs = 10.0 .^ range(-4 + 3/4, 2, length = 40)
    bεs = r0 .* 10 .^ range(log10(min_pairwise_distance(X)[2]/r0), stop = 0, length = 40)
    m, mεs = molteno_boxing(X)

    H2 = genentropy.(2, εs, Ref(X); base = 10.0)
    H1 = genentropy.(2, m, base = 10.0)
    C1 = [correlationsum(X, ε; w = 5) for ε in εs]
    C2 = boxed_correlationsum(X, bεs, r0)
    axs[1].plot(log10.(εs), -H2, label = name, color = color, ls = ls)
    #axs[2].plot(log10.(mεs), -H1, color = color)
    #axs[4].plot(log10.(bεs), log10.(C2))
    axs[4].set_xlabel("\$\\log_{10}(\\epsilon)\$")
    axs[4].set_ylabel("\$\\log_{10}(C_2)\$")
    axs[3].set_ylabel("\$\\log_{10}(C_1)\$")
    axs[2].set_ylabel("\$-H_1\$")
    axs[1].set_ylabel("\$-H_2\$")
    for (number, (set, es)) in enumerate(zip([10 .^ -H1, C1, C2], [mεs, εs, bεs]))
        i = findfirst(z -> z > 0, set)
        x, y = log10.(es)[i:end], log10.(set)[i:end]
        is, d = linear_region(x, y)
        axs[1+number].plot(x, y, zorder = 1, color = color, label = "\$ d = $(round(d, digits = 2))\$", ls = ls)
        axs[1+number].plot([x[is[1]], x[is[end]]], [y[is[1]], y[is[end]]],
        zorder = 2, ms = 5, marker = "o", ls = "None", color = color)
    end
end
for ax in axs; ax.legend(loc = "upper left"); ax.grid(); end
axs[1].tick_params(labelbottom=false)
axs[1].set_title("different datasets (N=$(N))")
fig.savefig("plots/comparing_differentedata.png")

# %% computation times for different lengths

dt = 0.1
Ns = 500:500:10000
times_lo67 = [[] for i in 1:4]
times_he = [[] for i in 1:4]
for (j, N) in enumerate(Ns)
    X1 = trajectory(lo67, N*dt; dt = dt, Ttr = 10.0)
    X2 = trajectory(he, N; Ttr = 10)
    # standard generalized entropy algorithm
    t_lo67 = @benchmark begin
        ε0 = min_pairwise_distance(X1)[2]
        ε1 = minimum(maxima(X1) - minima(X1))
        εs = 10.0 .^ range(log10(ε0), log10(ε1), length = 12)
        H = genentropy.(2, εs, Ref(X1); base = 10.0)
    end
    t_he = @benchmark begin
        ε0 = min_pairwise_distance(X2)[2]
        ε1 = minimum(maxima(X2) - minima(X2))
        εs = 10.0 .^ range(log10(ε0), log10(ε1), length = 12)
        H = genentropy.(2, εs, Ref(X2); base = 10.0)
    end
    push!(times_lo67[1], minimum(t_lo67.times))
    push!(times_he[1], minimum(t_he.times))

    # molteno's generalized entropy algorithm
    t_lo67 = @benchmark begin
        boxes, εs = molteno_boxing(X1)
        H = genentropy.(2, boxes; base = 10.0)
    end
    t_he = @benchmark begin
        boxes, εs = molteno_boxing(X2)
        H = genentropy.(2, boxes; base = 10.0)
    end
    push!(times_lo67[2], minimum(t_lo67.times))
    push!(times_he[2], minimum(t_he.times))

    # classic correlation algorithm
    t_lo67 = @benchmark begin
        ε0 = min_pairwise_distance(X1)[2]
        ε1 = minimum(maxima(X1) - minima(X1))
        εs = 10 .^ range(log10(ε0), log10(ε1), length = 12)
        C = log10.(correlationsum(X1, εs; w = 5))
    end
    t_he = @benchmark begin
        ε0 = min_pairwise_distance(X2)[2]
        ε1 = minimum(maxima(X2) - minima(X2))
        εs = 10 .^ range(log10(ε0), log10(ε1), length = 12)
        C = log10.(correlationsum(X2, εs; w = 5))
    end
    push!(times_lo67[3], minimum(t_lo67.times))
    push!(times_he[3], minimum(t_he.times))

    # boxed correlation algorithm
    t_lo67 = @benchmark begin
        r0 = estimate_r0_buenoorovio(X1)
        ε0 = min_pairwise_distance(X1)[2]
        εs = 10.0 .^ range(log10(ε0), log10(r0), length = 12)
        C = log10.(boxed_correlationsum(X1, εs, r0))
    end
    t_he = @benchmark begin
        r0 = estimate_r0_buenoorovio(X2)
        ε0 = min_pairwise_distance(X2)[2]
        εs = 10.0 .^ range(log10(ε0), log10(r0), length = 12)
        C = log10.(boxed_correlationsum(X2, εs, r0))
    end
    push!(times_lo67[4], minimum(t_lo67.times))
    push!(times_he[4], minimum(t_he.times))
end
fig, axs = subplots(2,1; sharex = true, figsize = (10, 16))
names = ["Generalized dimension" "Molteno's method" "Correlation dimension" "Boxed correlation dimension"]
for i in 1:4
    color = "C$(i-1)"
    name = names[i]
    axs[1].plot(Ns, times_lo67[i], label = name, color = color)
    axs[2].plot(Ns, times_he[i], color = color)
end
axs[2].set_xlabel("\$N\$")
axs[2].set_ylabel("\$t_{Hénon}\$")
axs[1].set_ylabel("\$t_{Lorenz63}\$")

for ax in axs; ax.legend(loc = "upper left"); ax.grid(); end
axs[1].tick_params(labelbottom=false)
axs[1].set_title("minimum computation time")
fig.savefig("plots/comparing_btime.png")

# %% sensitivity to noise

max_stars_repo_name: Datseris/FractalDimension

max_issues_repo_name: Datseris/FractalDimension

max_forks_repo_name: Datseris/FractalDimension


text:
using Test

#@testset "AlgebraicRelations" begin
#  include("AlgebraicRelations.jl")
#end

@testset "DB" begin
  include("DB.jl")
end

@testset "Queries" begin
  include("Queries.jl")
end

@testset "Presentations" begin
  include("Presentations.jl")
end

ext: jl


text:
# We provide a non-`Variable` implementation of the `AbstractVariable` interface
# to test that only the interface is used (and not, e.g. direct field access).
solver = MOI.OptimizerWithAttributes(SCS.Optimizer, "verbose" => 0)
TOL = 1e-3
module DictVectors
using Convex

# To make sure `Convex` isn't using field access on `AbstractVariable`'s
# we'll use a global dictionary to store information about each instance
# our of mock variable type, `DictVector`.
const global_cache = Dict{UInt64,Any}()

mutable struct DictVector{T} <: Convex.AbstractVariable
    head::Symbol
    id_hash::UInt64
    size::Tuple{Int,Int}
    function DictVector{T}(d) where {T}
        this = new(:DictVector, 0, (d, 1))
        this.id_hash = objectid(this)
        global_cache[this.id_hash] = Dict(
            :value => nothing,
            :sign => T <: Complex ? ComplexSign() : NoSign(),
            :vartype => ContVar,
            :constraints => Constraint[],
            :vexity => Convex.AffineVexity(),
        )
        return this
    end
end

Convex.evaluate(x::DictVector) = global_cache[x.id_hash][:value]

function Convex.set_value!(x::DictVector, v::AbstractArray)
    return global_cache[x.id_hash][:value] = v
end
function Convex.set_value!(x::DictVector, v::Number)
    return global_cache[x.id_hash][:value] = v
end

Convex.vexity(x::DictVector) = global_cache[x.id_hash][:vexity]
function Convex.vexity!(x::DictVector, v::Convex.Vexity)
    return global_cache[x.id_hash][:vexity] = v
end

Convex.sign(x::DictVector) = global_cache[x.id_hash][:sign]
Convex.sign!(x::DictVector, s::Convex.Sign) = global_cache[x.id_hash][:sign] = s

Convex.vartype(x::DictVector) = global_cache[x.id_hash][:vartype]
function Convex.vartype!(x::DictVector, s::Convex.VarType)
    return global_cache[x.id_hash][:vartype] = s
end

Convex.constraints(x::DictVector) = global_cache[x.id_hash][:constraints]
function Convex.add_constraint!(x::DictVector, s::Convex.Constraint)
    return push!(global_cache[x.id_hash][:constraints], s)
end

end

import .DictVectors

@testset "AbstractVariable interface" begin
    # Let us solve a basic problem from `test_affine.jl`

    x = DictVectors.DictVector{BigFloat}(1)
    y = DictVectors.DictVector{BigFloat}(1)
    p = minimize(x + y, [x >= 3, y >= 2])
    @test vexity(p) == Convex.AffineVexity()
    solve!(p, solver)
    @test p.optval ≈ 5 atol = TOL
    @test evaluate(x + y) ≈ 5 atol = TOL

    add_constraint!(x, x >= 4)
    solve!(p, solver)
    @test p.optval ≈ 6 atol = TOL
    @test evaluate(x + y) ≈ 6 atol = TOL
    @test length(constraints(x)) == 1
end

# Let us do another example of custom variable types, but using field access for simplicity

module DensityMatricies
using Convex

mutable struct DensityMatrix{T} <: Convex.AbstractVariable
    head::Symbol
    id_hash::UInt64
    size::Tuple{Int,Int}
    value::Convex.ValueOrNothing
    vexity::Convex.Vexity
    function DensityMatrix(d)
        this = new{ComplexF64}(
            :DensityMatrix,
            0,
            (d, d),
            nothing,
            Convex.AffineVexity(),
        )
        this.id_hash = objectid(this)
        return this
    end
end
Convex.constraints(ρ::DensityMatrix) = [ρ ⪰ 0, tr(ρ) == 1]
Convex.sign(::DensityMatrix) = Convex.ComplexSign()
Convex.vartype(::DensityMatrix) = Convex.ContVar

end

import .DensityMatricies
import LinearAlgebra
@testset "DensityMatrix custom variable type" begin
    X = randn(4, 4) + im * rand(4, 4)
    X = X + X'
    # `X` is Hermitian and non-degenerate (with probability 1)
    # Let us calculate the projection onto the eigenspace associated
    # to the maximum eigenvalue
    e_vals, e_vecs = LinearAlgebra.eigen(LinearAlgebra.Hermitian(X))
    e_val, idx = findmax(e_vals)
    e_vec = e_vecs[:, idx]
    proj = e_vec * e_vec'

    # found it! Now let us do it again via an SDP
    ρ = DensityMatricies.DensityMatrix(4)

    prob = maximize(real(tr(ρ * X)))
    solve!(prob, solver)

    @test prob.optval ≈ e_val atol = TOL
    @test evaluate(ρ) ≈ proj atol = TOL
end

module ProbabilityVectors
using Convex
mutable struct ProbabilityVector <: Convex.AbstractVariable
    head::Symbol
    id_hash::UInt64
    size::Tuple{Int,Int}
    value::Convex.ValueOrNothing
    vexity::Convex.Vexity
    function ProbabilityVector(d)
        this =
            new(:ProbabilityVector, 0, (d, 1), nothing, Convex.AffineVexity())
        this.id_hash = objectid(this)
        return this
    end
end
Convex.constraints(p::ProbabilityVector) = [sum(p) == 1]
Convex.sign(::ProbabilityVector) = Convex.Positive()
Convex.vartype(::ProbabilityVector) = Convex.ContVar

(p::ProbabilityVector)(x) = dot(p, x)

end

using .ProbabilityVectors

@testset "ProbabilityVectors" begin
    p = ProbabilityVectors.ProbabilityVector(3)
    x = [1.0, 2.0, 3.0]

    @test p(x) isa AbstractExpr
    @test sign(p) == Positive()
    prob = minimize(p(x))
    solve!(prob, solver)
    @test prob.optval ≈ 1.0 atol = TOL
    @test evaluate(p(x)) ≈ 1.0 atol = TOL
    @test evaluate(p) ≈ [1.0, 0.0, 0.0] atol = TOL
end


text:
"""
    ProjectTo(x::T)

Returns a `ProjectTo{T,...}` functor able to project a differential `dx` onto the same tangent space as `x`.
This functor encloses over what ever is needed to be able to be able to do that projection.
For example, when projecting `dx=ZeroTangent()` on an array `P=Array{T, N}`, the size of `x`
is not available from `P`, so it is stored in the functor.

    (::ProjectTo{T})(dx)

Projects the differential `dx` on the onto the tangent space used to create the `ProjectTo`.
"""
struct ProjectTo{P,D<:NamedTuple}
    info::D
end
ProjectTo{P}(info::D) where {P,D<:NamedTuple} = ProjectTo{P,D}(info)
ProjectTo{P}(; kwargs...) where {P} = ProjectTo{P}(NamedTuple(kwargs))

backing(project::ProjectTo) = getfield(project, :info)
Base.getproperty(p::ProjectTo, name::Symbol) = getproperty(backing(p), name)
Base.propertynames(p::ProjectTo) = propertynames(backing(p))

function Base.show(io::IO, project::ProjectTo{T}) where {T}
    print(io, "ProjectTo{")
    show(io, T)
    print(io, "}")
    if isempty(backing(project))
        print(io, "()")
    else
        show(io, backing(project))
    end
end

# fallback (structs)
function ProjectTo(x::T) where {T}
    # Generic fallback for structs, recursively make `ProjectTo`s all their fields
    fields_nt::NamedTuple = backing(x)
    return ProjectTo{T}(map(ProjectTo, fields_nt))
end
function (project::ProjectTo{T})(dx::Tangent) where {T}
    sub_projects = backing(project)
    sub_dxs = backing(canonicalize(dx))
    _call(f, x) = f(x)
    return construct(T, map(_call, sub_projects, sub_dxs))
end

# should not work for Tuples and NamedTuples, as not valid tangent types
function ProjectTo(x::T) where {T<:Union{<:Tuple,NamedTuple}}
    return throw(
        ArgumentError("The `x` in `ProjectTo(x)` must be a valid differential, not $x")
    )
end

# Generic
(project::ProjectTo)(dx::AbstractThunk) = project(unthunk(dx))
(::ProjectTo{T})(dx::T) where {T} = dx  # not always true, but we can special case for when it isn't
(::ProjectTo{T})(dx::AbstractZero) where {T} = zero(T)

# Number
ProjectTo(::T) where {T<:Number} = ProjectTo{T}()
(::ProjectTo{T})(dx::Number) where {T<:Number} = convert(T, dx)
(::ProjectTo{T})(dx::Number) where {T<:Real} = convert(T, real(dx))

# Arrays
ProjectTo(xs::T) where {T<:Array} = ProjectTo{T}(; elements=map(ProjectTo, xs))
function (project::ProjectTo{T})(dx::Array) where {T<:Array}
    _call(f, x) = f(x)
    return T(map(_call, project.elements, dx))
end
function (project::ProjectTo{T})(dx::AbstractZero) where {T<:Array}
    return T(map(proj -> proj(dx), project.elements))
end
(project::ProjectTo{<:Array})(dx::AbstractArray) = project(collect(dx))

# Arrays{<:Number}: optimized case so we don't need a projector per element
function ProjectTo(x::T) where {E<:Number,T<:Array{E}}
    return ProjectTo{T}(; element=ProjectTo(zero(E)), size=size(x))
end
(project::ProjectTo{<:Array{T}})(dx::Array) where {T<:Number} = project.element.(dx)
function (project::ProjectTo{<:Array{T}})(dx::AbstractZero) where {T<:Number}
    return zeros(T, project.size)
end
function (project::ProjectTo{<:Array{T}})(dx::Tangent{<:SubArray}) where {T<:Number}
    return project(dx.parent)
end

# Diagonal
ProjectTo(x::T) where {T<:Diagonal} = ProjectTo{T}(; diag=ProjectTo(diag(x)))
(project::ProjectTo{T})(dx::AbstractMatrix) where {T<:Diagonal} = T(project.diag(diag(dx)))
(project::ProjectTo{T})(dx::AbstractZero) where {T<:Diagonal} = T(project.diag(dx))

# :data, :uplo fields
for SymHerm in (:Symmetric, :Hermitian)
    @eval begin
        function ProjectTo(x::T) where {T<:$SymHerm}
            return ProjectTo{T}(; uplo=Symbol(x.uplo), parent=ProjectTo(parent(x)))
        end
        function (project::ProjectTo{<:$SymHerm})(dx::AbstractMatrix)
            return $SymHerm(project.parent(dx), project.uplo)
        end
        function (project::ProjectTo{<:$SymHerm})(dx::AbstractZero)
            return $SymHerm(project.parent(dx), project.uplo)
        end
        function (project::ProjectTo{<:$SymHerm})(dx::Tangent)
            return $SymHerm(project.parent(dx.data), project.uplo)
        end
    end
end

# :data field
for UL in (:UpperTriangular, :LowerTriangular)
    @eval begin
        ProjectTo(x::T) where {T<:$UL} = ProjectTo{T}(; parent=ProjectTo(parent(x)))
        (project::ProjectTo{<:$UL})(dx::AbstractMatrix) = $UL(project.parent(dx))
        (project::ProjectTo{<:$UL})(dx::AbstractZero) = $UL(project.parent(dx))
        (project::ProjectTo{<:$UL})(dx::Tangent) = $UL(project.parent(dx.data))
    end
end

# Transpose
ProjectTo(x::T) where {T<:Transpose} = ProjectTo{T}(; parent=ProjectTo(parent(x)))
function (project::ProjectTo{<:Transpose})(dx::AbstractMatrix)
    return transpose(project.parent(transpose(dx)))
end
(project::ProjectTo{<:Transpose})(dx::AbstractZero) = transpose(project.parent(dx))

# Adjoint
ProjectTo(x::T) where {T<:Adjoint} = ProjectTo{T}(; parent=ProjectTo(parent(x)))
(project::ProjectTo{<:Adjoint})(dx::AbstractMatrix) = adjoint(project.parent(adjoint(dx)))
(project::ProjectTo{<:Adjoint})(dx::AbstractZero) = adjoint(project.parent(dx))

# PermutedDimsArray
ProjectTo(x::P) where {P<:PermutedDimsArray} = ProjectTo{P}(; parent=ProjectTo(parent(x)))
function (project::ProjectTo{<:PermutedDimsArray{T,N,perm,iperm,AA}})(
    dx::AbstractArray
) where {T,N,perm,iperm,AA}
    return PermutedDimsArray{T,N,perm,iperm,AA}(permutedims(project.parent(dx), perm))
end
function (project::ProjectTo{P})(dx::AbstractZero) where {P<:PermutedDimsArray}
    return P(project.parent(dx))
end

# SubArray
ProjectTo(x::T) where {T<:SubArray} = ProjectTo(copy(x))  # don't project on to a view, but onto matching copy


text:
@testset "Kruskal" begin
    g4 = CompleteGraph(4)

    distmx = [
        0  1  5  6
        1  0  4  10
        5  4  0  3
        6  10  3  0
    ]

    vec_mst = Vector{Edge}([Edge(1, 2), Edge(3, 4), Edge(2, 3)])
    max_vec_mst = Vector{Edge}([Edge(2, 4), Edge(1, 4), Edge(1, 3)])
    for g in testgraphs(g4)
    # Testing Kruskal's algorithm
        mst = @inferred(kruskal_mst(g, distmx))
        @test mst == vec_mst
	@test @inferred(kruskal_mst(g, distmx, minimize=false)) == max_vec_mst
    end
    #second test
    distmx_sec = [
        0       0       0.26    0       0.38    0       0.58    0.16
        0       0       0.36    0.29    0       0.32    0       0.19
        0.26    0.36    0       0.17    0       0       0.4     0.34
        0       0.29    0.17    0       0       0       0.52    0
        0.38    0       0       0       0       0.35    0.93    0.37
        0       0.32    0       0       0.35    0       0       0.28
        0.58    0       0.4     0.52    0.93    0       0       0
        0.16    0.19    0.34    0       0.37    0.28    0       0
    ]

    gx = SimpleGraph(distmx_sec)
    vec2 = Vector{Edge}([Edge(1, 8), Edge(3, 4), Edge(2, 8), Edge(1, 3), Edge(6, 8), Edge(5, 6), Edge(3, 7)])
    max_vec2 = Vector{Edge}([Edge(5, 7), Edge(1, 7), Edge(4, 7), Edge(3, 7), Edge(5, 8), Edge(2, 3), Edge(5, 6)])
    for g in testgraphs(gx)
        mst2 = @inferred(kruskal_mst(g, distmx_sec))
        @test mst2 == vec2
	@test @inferred(kruskal_mst(g, distmx_sec, minimize=false)) == max_vec2
    end
end


text:
######
##    TableRegressionModel methods to forward
######
@delegate TableRegressionModel.model [
    leverage,
    weights,
    workingweights,
    dispersion,
    scale,
    tauscale,
    fitted,
    isfitted,
    islinear,
    Estimator,
    projectionmatrix,
]

fit!(p::TableRegressionModel, args...; kwargs...) = (fit!(p.model, args...; kwargs...); p)
refit!(p::TableRegressionModel, args...; kwargs...) =
    (refit!(p.model, args...; kwargs...); p)

"""
    RobustLinearModel

Robust linear model representation

## Fields

* `resp`: the [`RobustLinResp`](@ref) structure.
* `pred`: the predictor structure, of type [`DensePredChol`](@ref), [`SparsePredChol`](@ref), [`DensePredCG`](@ref), [`SparsePredCG`](@ref) or [`RidgePred`](@ref).
* `fitdispersion`: if true, the dispersion is estimated otherwise it is kept fixed
* `fitted`: if true, the model was already fitted
"""
mutable struct RobustLinearModel{T<:AbstractFloat,R<:RobustResp{T},L<:LinPred} <:
               AbstractRobustModel{T}
    resp::R
    pred::L
    fitdispersion::Bool
    fitted::Bool
end

function Base.getproperty(
    r::TableRegressionModel{M},
    s::Symbol,
) where {M<:RobustLinearModel}
    if s ∈ (:resp, :pred, :fitdispersion, :fitted)
        getproperty(r.model, s)
    else
        getfield(r, s)
    end
end

######
##    AbstractRobustModel methods
######

dof(m::AbstractRobustModel) = length(coef(m))

dof_residual(m::AbstractRobustModel) = nobs(m) - dof(m)

function coeftable(m::AbstractRobustModel; level::Real=0.95)
    cc = coef(m)
    se = stderror(m)
    tt = cc ./ se
    ci = se * quantile(TDist(dof_residual(m)), (1 - level) / 2)
    p = ccdf.(Ref(FDist(1, dof_residual(m))), abs2.(tt))
    levstr = isinteger(level * 100) ? string(Integer(level * 100)) : string(level * 100)
    CoefTable(
        hcat(cc, se, tt, p, cc + ci, cc - ci),
        ["Coef.", "Std. Error", "t", "Pr(>|t|)", "Lower $(levstr)%", "Upper $(levstr)%"],
        ["x$i" for i in 1:length(cc)],
        4,
        3,
    )
end

function confint(m::AbstractRobustModel; level::Real=0.95)
    alpha = quantile(TDist(dof_residual(m)), (1 - level) / 2)
    hcat(coef(m), coef(m)) + stderror(m) * alpha * hcat(1.0, -1.0)
end
confint(m::AbstractRobustModel, level::Real) = confint(m; level=level)

## TODO: specialize to make it faster
leverage(p::AbstractRobustModel) = diag(projectionmatrix(p))

######
##    RobustLinearModel methods
######

function show(io::IO, obj::RobustLinearModel)
    println(io, "Robust regression with $(obj.resp.est)\n\nCoefficients:\n", coeftable(obj))
end

function show(io::IO, obj::TableRegressionModel{M,T}) where {T,M<:RobustLinearModel}
    println(
        io,
        "Robust regression with $(obj.model.resp.est)\n\n$(obj.mf.f)\n\nCoefficients:\n",
        coeftable(obj),
    )
end

islinear(m::RobustLinearModel) = true

"""
    deviance(m::RobustLinearModel)

The sum of twice the loss/objective applied to the scaled residuals.

It is consistent with the definition of the deviance for OLS.
"""
deviance(m::RobustLinearModel) = deviance(m.resp)

nulldeviance(m::RobustLinearModel) = nulldeviance(m.resp)

"""
    dispersion(m::RobustLinearModel, sqr::Bool=false)

The dispersion is the (weighted) sum of robust residuals. If `sqr` is true, return the squared dispersion.
"""
dispersion(m::RobustLinearModel, sqr::Bool=false) = dispersion(m.resp, dof_residual(m), sqr)

nobs(m::RobustLinearModel)::Int = nobs(m.resp)

coef(m::RobustLinearModel) = coef(m.pred)

"""
    Estimator(m::RobustLinearModel)

The robust estimator object used to fit the model.
"""
Estimator(m::RobustLinearModel) = Estimator(m.resp)

stderror(m::RobustLinearModel) =
    location_variance(m.resp, dof_residual(m), false) .* sqrt.(diag(vcov(m)))

loglikelihood(m::RobustLinearModel) = loglikelihood(m.resp)

nullloglikelihood(m::RobustLinearModel) = nullloglikelihood(m.resp)

weights(m::RobustLinearModel) = weights(m.resp)

"""
    workingweights(m::RobustLinearModel)

The robust weights computed by the model.

This can be used to detect outliers, as outliers weights are lower than the
weights of valid data points.
"""
workingweights(m::RobustLinearModel) = workingweights(m.resp)

response(m::RobustLinearModel) = response(m.resp)

isfitted(m::RobustLinearModel) = m.fitted

fitted(m::RobustLinearModel) = fitted(m.resp)

residuals(m::RobustLinearModel) = residuals(m.resp)

"""
    scale(m::RobustLinearModel, sqr::Bool=false)

The robust scale estimate used for the robust estimation.

If `sqr` is `true`, the square of the scale is returned.
"""
scale(m::RobustLinearModel, sqr::Bool=false) = scale(m.resp, sqr)

"""
    tauscale(m::RobustLinearModel, sqr::Bool=false; kwargs...)

The robust τ-scale that is minimized in τ-estimation.

If `sqr` is `true`, the square of the τ-scale is returned.
"""
tauscale(m::RobustLinearModel, args...; kwargs...) = tauscale(m.resp, args...; kwargs...)

function predict(
    m::RobustLinearModel,
    newX::AbstractMatrix;
    offset::FPVector=eltype(newX)[],
)
    mu = newX * coef(m)
    if !isempty(m.resp.offset)
        if !(length(offset) == size(newX, 1))
            mess =
                "fit with offset, so `offset` keyword argument" *
                " must be an offset of length `size(newX, 1)`"
            throw(ArgumentError(mess))
        end
        broadcast!(+, mu, mu, offset)
    else
        if length(offset) > 0
            mess = "fit without offset, so value of `offset` kwarg does not make sense"
            throw(ArgumentError(mess))
        end
    end
    mu
end
predict(m::RobustLinearModel) = fitted(m)

"""
    rlm(X, y, args...; kwargs...)

An alias for `fit(RobustLinearModel, X, y, est; kwargs...)`.

The arguments `X` and `y` can be a `Matrix` and a `Vector` or a `Formula` and a `DataFrame`.
"""
rlm(X, y, args...; kwargs...) = fit(RobustLinearModel, X, y, args...; kwargs...)

"""
    fit(::Type{M},
        X::Union{AbstractMatrix{T},SparseMatrixCSC{T}},
        y::AbstractVector{T},
        est::Estimator;
        method::Symbol       = :chol, # :cg
        dofit::Bool          = true,
        wts::FPVector        = similar(y, 0),
        offset::FPVector     = similar(y, 0),
        fitdispersion::Bool  = false,
        ridgeλ::Real         = 0,
        ridgeG::Union{UniformScaling, AbstractArray} = I,
        βprior::AbstractVector = [],
        quantile::Union{Nothing, AbstractFloat} = nothing,
        initial_scale::Union{Symbol, Real}=:mad,
        σ0::Union{Nothing, Symbol, Real}=initial_scale,
        initial_coef::AbstractVector=[],
        β0::AbstractVector=initial_coef,
        correct_leverage::Bool=false
        fitargs...) where {M<:RobustLinearModel, T<:AbstractFloat}

Create a robust model with the model matrix (or formula) X and response vector (or dataframe) y,
using a robust estimator.

# Arguments

- `X`: the model matrix (it can be dense or sparse) or a formula
- `y`: the response vector or a dataframe.
- `est`: a robust estimator

# Keywords

- `method::Symbol = :chol`: the method to use for solving the weighted linear system, `chol` (default) or `cg`;
- `dofit::Bool = true`: if false, return the model object without fitting;
- `wts::Vector = []`: a weight vector, should be empty if no weights are used;
- `offset::Vector = []`: an offset vector, should be empty if no offset is used;
- `fitdispersion::Bool = false`: reevaluate the dispersion;
- `ridgeλ::Real = 0`: if positive, perform a robust ridge regression with shrinkage parameter `ridgeλ`. [`RidgePred`](@ref) object will be used;
- `ridgeG::Union{UniformScaling, AbstractArray} = I`: define a custom regularization matrix. Default to unity matrix (with 0 for the intercept);
- `βprior::AbstractVector = []`: define a custom prior for the coefficients for ridge regression. Default to `zeros(p)`;
- `quantile::Union{Nothing, AbstractFloat} = nothing`: only for [`GeneralizedQuantileEstimator`](@ref), define the quantile to estimate;
- `initial_scale::Union{Symbol, Real}=:mad`: the initial scale estimate, for non-convex estimator it helps to find the global minimum. Automatic computation using `:mad`, `L1` or `extrema` (non-robust).
- `σ0::Union{Nothing, Symbol, Real}=initial_scale`: alias of `initial_scale`;
- `initial_coef::AbstractVector=[]`: the initial coefficients estimate, for non-convex estimator it helps to find the global minimum.
- `β0::AbstractVector=initial_coef`: alias of `initial_coef`;
- `correct_leverage::Bool=false`: apply the leverage correction weights with [`leverage_weights`](@ref).
- `fitargs...`: other keyword arguments used to control the convergence of the IRLS algorithm (see [`pirls!`](@ref)).

# Output

the RobustLinearModel object.

"""
function fit(
    ::Type{M},
    X::Union{AbstractMatrix{T},SparseMatrixCSC{T}},
    y::AbstractVector{T},
    est::V;
    method::Symbol=:chol, # :cg
    dofit::Bool=true,
    wts::FPVector=similar(y, 0),
    offset::FPVector=similar(y, 0),
    fitdispersion::Bool=false,
    ridgeλ::Real=0,
    ridgeG::Union{UniformScaling,AbstractArray}=I,
    βprior::AbstractVector=[],
    quantile::Union{Nothing,AbstractFloat}=nothing,
    fitargs...,
) where {M<:RobustLinearModel,V<:AbstractEstimator,T<:AbstractFloat}

    # Check that X and y have the same number of observations
    n, p = size(X)
    if n != size(y, 1)
        throw(DimensionMismatch("number of rows in X and y must match"))
    end

    # Change quantile
    if !isnothing(quantile)
        if !isa(est, AbstractQuantileEstimator)
            throw(
                TypeError(
                    :fit,
                    "arguments, quantile can be changed only if isa(est, AbstractQuantileEstimator)",
                    AbstractQuantileEstimator,
                    est,
                ),
            )
        end
        est.quantile = quantile
    end

    # Response object
    rr = RobustLinResp(est, y, offset, wts)

    # Predictor object
    pp = if ridgeλ > 0
        # With ridge regularization
        G = if isa(ridgeG, UniformScaling)
            # Has an intersect
            if all(X[:, 1] .== 1)
                spdiagm(0 => [float(i != 1) for i in 1:p])
            else
                I(p)
            end
        else
            ridgeG
        end
        if method == :cg
            cgpred(X, float(ridgeλ), G, βprior)
        else
            cholpred(X, float(ridgeλ), G, βprior)
        end
    else
        # No regularization
        if method == :cg
            cgpred(X)
        else
            cholpred(X)
        end
    end

    m = RobustLinearModel(rr, pp, fitdispersion, false)
    return dofit ? fit!(m; fitargs...) : m
end

function fit(
    ::Type{M},
    X::Union{AbstractMatrix,SparseMatrixCSC},
    y::AbstractVector,
    est::AbstractEstimator;
    kwargs...,
) where {M<:AbstractRobustModel}
    fit(M, float(X), float(y), est; kwargs...)
end

"""
    refit!(m::RobustLinearModel, [y::FPVector];
                                 wts::Union{Nothing, FPVector} = nothing,
                                 offset::Union{Nothing, FPVector} = nothing,
                                 quantile::Union{Nothing, AbstractFloat} = nothing,
                                 ridgeλ::Union{Nothing, Real} = nothing,
                                 kwargs...)

Refit the [`RobustLinearModel`](@ref).

This function assumes that `m` was correctly initialized and the model is refitted with
the new values for the response, weights, offset, quantile and ridge shrinkage.

Defining a new `quantile` is only possible for [`GeneralizedQuantileEstimator`](@ref).

Defining a new `ridgeλ` is only possible for [`RidgePred`](@ref) objects.
"""
function refit!(m::RobustLinearModel, y::FPVector; kwargs...)
    r = m.resp
    # Check that old and new y have the same number of observations
    if size(r.y, 1) != size(y, 1)
        mess = "the new response vector should have the same dimension: "*
               "$(size(r.y, 1)) != $(size(y, 1))"
        throw(DimensionMismatch(mess))
    end
    # Update y
    copyto!(r.y, y)

    refit!(m; kwargs...)
end

function refit!(
    m::RobustLinearModel{T};
    wts::Union{Nothing,FPVector}=nothing,
    offset::Union{Nothing,FPVector}=nothing,
    quantile::Union{Nothing,AbstractFloat}=nothing,
    ridgeλ::Union{Nothing,Real}=nothing,
    kwargs...,
) where {T}

    if haskey(kwargs, :method)
        @warn("the method cannot be changed when refitting,"*
              " ignore the method argument $(kwargs[:method])."
        )
        delete!(kwargs, :method)
    end

    r = m.resp

    n = length(r.y)
    if !isa(wts, Nothing) && (length(wts) in (0, n))
        copy!(r.wts, wts)
    end
    if !isa(offset, Nothing) && (length(offset) in (0, n))
        copy!(r.offset, offset)
    end

    # Update quantile, if the estimator is AbstractQuantileEstimator
    if !isnothing(quantile)
        isa(r.est, AbstractQuantileEstimator) || throw(
            TypeError(
                :refit!,
                "arguments, quantile can be changed only if isa(r.est, AbstractQuantileEstimator)",
                AbstractQuantileEstimator,
                r.est,
            ),
        )
        r.est.quantile = quantile
    end

    # Update ridge shrinkage parameter
    if !isnothing(ridgeλ)
        isa(m.pred, RidgePred) || throw(
            TypeError(
                :refit!,
                "arguments, ridgeλ can be changed only if the predictor is a RidgePred",
                RidgePred,
                r.pred,
            ),
        )
        m.pred.λ = float(ridgeλ)
        # reset beta0 because it needs to be zero for the first estimation
        resetβ0!(m.pred)
    end

    # Reinitialize the coefficients and the response
    fill!(coef(m), zero(T))
    initresp!(r)

    m.fitted = false
    fit!(m; kwargs...)
end

"""
    fit!(m::RobustLinearModel; initial_scale::Union{Symbol, Real}=:mad,
              σ0::Union{Nothing, Symbol, Real}=initial_scale,
              initial_coef::AbstractVector=[],
              β0::AbstractVector=initial_coef,
              correct_leverage::Bool=false, kwargs...)

Optimize the objective of a `RobustLinearModel`.  When `verbose` is `true` the values of the
objective and the parameters are printed on stdout at each iteration.

This function assumes that `m` was correctly initialized.

This function returns early if the model was already fitted, instead call `refit!`.
"""
function fit!(
    m::RobustLinearModel;
    initial_scale::Union{Symbol,Real}=:mad,
    σ0::Union{Nothing,Symbol,Real}=initial_scale,
    initial_coef::AbstractVector=[],
    β0::AbstractVector=initial_coef,
    correct_leverage::Bool=false,
    kwargs...,
)

    # Return early if model has the fit flag set
    m.fitted && return m

    # Compute the initial values
    σ0, β0 = process_σ0β0(m, σ0, β0)

    if correct_leverage
        wts = m.resp.wts
        copy!(wts, leverage_weights(m))
        ## TODO: maybe multiply by the old wts?
    end

    # Get type
    V = typeof(m.resp.est)

    _fit!(m, V; σ0=σ0, β0=β0, kwargs...)

    m.fitted = true
    m
end

## Error message
function _fit!(m::RobustLinearModel, ::Type{E}; kwargs...) where {E<:AbstractEstimator}
    allowed_estimators =
        (MEstimator, SEstimator, MMEstimator, TauEstimator, GeneralizedQuantileEstimator)
    mess = "only types $(allowed_estimators) are allowed, "*
           "you must define the `_fit!` method for the type: $(E)"
    error(mess)
end

# Fit M-estimator
function _fit!(
    m::RobustLinearModel,
    ::Type{E};
    σ0::AbstractFloat=1.0,
    β0::AbstractVector=[],
    verbose::Bool=false,
    kwargs...,
) where {E<:MEstimator}

    verbose && println("\nFit with M-estimator: $(Estimator(m))")
    ## Minimize the objective
    pirls!(m; sigma0=σ0, beta0=β0, verbose=verbose, kwargs...)

    ## TODO: update scale is fitdispersion is true
    m
end

# Fit Generalized M-Quantile estimator
function _fit!(
    m::RobustLinearModel,
    ::Type{E};
    σ0::AbstractFloat=1.0,
    β0::AbstractVector=[],
    verbose::Bool=false,
    kwargs...,
) where {E<:GeneralizedQuantileEstimator}

    verbose && println("\nFit with M-Quantile-estimator: $(Estimator(m))")
    ## Minimize the objective
    pirls!(m; sigma0=σ0, beta0=β0, verbose=verbose, kwargs...)

    ## TODO: update scale if fitdispersion is true
    m
end

# Fit S-estimator
function _fit!(
    m::RobustLinearModel,
    ::Type{E};
    σ0::AbstractFloat=1.0,
    β0::AbstractVector=[],
    verbose::Bool=false,
    resample::Bool=false,
    resampling_options::Dict{Symbol,F}=Dict{Symbol,Any}(:verbose => verbose),
    kwargs...,
) where {F,E<:SEstimator}

    ## Resampling algorithm to find a starting point close to the global minimum
    if resample
        σ0, β0 = resampling_best_estimate(m, E; resampling_options...)
    end

    verbose && println("\nFit with S-estimator: $(Estimator(m))")
    ## Minimize the objective
    pirls_Sestimate!(m; sigma0=σ0, beta0=β0, verbose=verbose, kwargs...)

    # Set the `fitdispersion` flag to true, because σ was estimated
    m.fitdispersion = true

    m
end

# Fit MM-estimator
function _fit!(
    m::RobustLinearModel,
    ::Type{E};
    σ0::AbstractFloat=1.0,
    β0::AbstractVector=[],
    verbose::Bool=false,
    resample::Bool=false,
    resampling_options::Dict{Symbol,F}=Dict{Symbol,Any}(:verbose => verbose),
    kwargs...,
) where {F,E<:MMEstimator}

    ## Set the S-Estimator for robust estimation of σ and β0
    set_SEstimator(Estimator(m.resp))

    ## Resampling algorithm to find a starting point close to the global minimum
    if resample
        σ0, β0 = resampling_best_estimate(m, E; resampling_options...)
    end

    verbose && println("\nFit with MM-estimator - 1. S-estimator: $(Estimator(m.resp))")
    ## Minimize the objective
    pirls_Sestimate!(m; sigma0=σ0, beta0=β0, verbose=verbose, kwargs...)

    ## Use an M-estimate to estimate coefficients
    β0 = coef(m)
    σ0 = scale(m)

    ## Set the M-Estimator for efficient estimation of β
    set_MEstimator(Estimator(m.resp))

    verbose && println("\nFit with MM-estimator - 2. M-estimator: $(Estimator(m.resp))")
    ## Minimize the objective
    pirls!(m; sigma0=σ0, beta0=β0, verbose=verbose, kwargs...)

    # Set the `fitdispersion` flag to true, because σ was estimated
    m.fitdispersion = true

    m
end

# Fit τ-estimator
function _fit!(
    m::RobustLinearModel,
    ::Type{E};
    σ0::AbstractFloat=1.0,
    β0::AbstractVector=[],
    verbose::Bool=false,
    resample::Bool=false,
    resampling_options::Dict{Symbol,F}=Dict{Symbol,Any}(:verbose => verbose),
    kwargs...,
) where {F,E<:TauEstimator}

    ## Resampling algorithm to find a starting point close to the global minimum
    if resample
        σ0, β0 = resampling_best_estimate(m, E; resampling_options...)
    end

    verbose && println("\nFit with τ-estimator: $(Estimator(m))")
    ## Minimize the objective
    pirls_τestimate!(m; sigma0=σ0, beta0=β0, verbose=verbose, kwargs...)

    # Set the `fitdispersion` flag to true, because σ was estimated
    m.fitdispersion = true

    m
end

function process_σ0β0(
    m::RobustLinearModel,
    σ0::Union{Real,Symbol}=scale(m),
    β0::AbstractVector=[],
)
    # Process scale σ0
    σ0 = if isa(σ0, Real)
        float(σ0)
    elseif isa(σ0, Symbol)
        initialscale(m, σ0)
    end

    # Process coefficients β0
    β0 = if isempty(β0) || size(β0, 1) != size(coef(m), 1)
        []
    else
        float(β0)
    end

    return σ0, β0
end

function initialscale(m::RobustLinearModel, method::Symbol=:mad; factor::AbstractFloat=1.0)
    factor > 0 || error("factor should be positive")

    y = response(m)
    wts = Vector(weights(m))

    allowed_methods = (:mad, :extrema, :L1)
    if method == :mad
        σ = if length(wts) == length(y)
            factor * mad(wts .* y; normalize=true)
        else
            factor * mad(y; normalize=true)
        end
    elseif method == :L1
        X = modelmatrix(m)
        σ = dispersion(quantreg(X, y; wts=wts))
    elseif method == :extrema
        # this is not robust
        σ = -(-(extrema(y)...)) / 2
    else
        error("only $(join(allowed_methods, ", ", " and ")) methods are allowed")
    end
    return σ
end

function setβ0!(m::RobustLinearModel{T}, β0::AbstractVector=[]) where {T<:AbstractFloat}
    r = m.resp
    p = m.pred

    initresp!(r)
    if isempty(β0)
        # Compute beta0 from solving the least square with the response value r.y
        # initresp!(r)
        delbeta!(p, r.wrkres, r.wrkwt)
        installbeta!(p)
    else
        copyto!(p.beta0, float(β0))
        fill!(p.delbeta, 0)
    end

    m
end

"""
    setinitη!(m)
Compute the predictor using the initial value of β0 and compute the residuals
"""
function setinitη!(m::RobustLinearModel{T}) where {T}
    r = m.resp
    p = m.pred

    ## Initially, β0 is defined but not ∇β, so use f=0
    linpred!(r.η, p, 0)
    updateres!(r; updatescale=false)

    m
end

"""
    setinitσ!(m)
Compute the predictor scale using the MAD of the residuals
Use only for rough estimate, like in the resampling phase.
"""
function setinitσ!(m::RobustLinearModel; kwargs...)
    m.resp.σ = madresidualscale(m.resp; kwargs...)
    m
end

"""
    setη!(m, f=1.0; updatescale=false, kwargs...)
Compute the ∇β using the current residuals and working weights (only if f=1,
which corresponds to the first iteration of linesearch), then compute
the predictor using the ∇β value and compute the new residuals and deviance.
The scaletype argument defines if the location or scale loss function should be used
If updatescale is true, the scale is also updated along with the residuals.
"""
function setη!(
    m::RobustLinearModel{T},
    f::T=1.0;
    updatescale::Bool=false,
    kwargs...,
) where {T}
    r = m.resp
    p = m.pred

    # First update of linesearch algorithm, compute ∇β
    if f == 1
        delbeta!(p, r.wrkres, r.wrkwt)
    end
    # Compute and set the predictor η from β0 and ∇β
    linpred!(r.η, p, f)

    # Update the residuals and weights (and scale if updatescale=true)
    updateres!(r; updatescale=updatescale, kwargs...)
    m
end

"""
    pirls!(m::RobustLinearModel{T}; verbose::Bool=false, maxiter::Integer=30,
           minstepfac::Real=1e-3, atol::Real=1e-6, rtol::Real=1e-6,
           beta0::AbstractVector=[], sigma0::Union{Nothing, T}=nothing)

(Penalized) Iteratively Reweighted Least Square procedure for M-estimation.
The Penalized aspect is not implemented (yet).
"""
function pirls!(
    m::RobustLinearModel{T};
    verbose::Bool=false,
    maxiter::Integer=30,
    minstepfac::Real=1e-3,
    atol::Real=1e-6,
    rtol::Real=1e-5,
    beta0::AbstractVector=[],
    sigma0::Union{Nothing,T}=nothing,
) where {T<:AbstractFloat}

    # Check arguments
    maxiter >= 1 || throw(ArgumentError("maxiter must be positive"))
    0 < minstepfac < 1 || throw(ArgumentError("minstepfac must be in (0, 1)"))

    # Extract fields and set convergence flag
    cvg, p, r = false, m.pred, m.resp

    ## Initialize σ, default to do not change
    if !isnothing(sigma0)
        r.σ = sigma0
    end

    ## Initialize β or set it to the provided values
    setβ0!(m, beta0)

    # Initialize μ and compute residuals
    setinitη!(m)

    # If σ==0, iterations will fail, so return here
    if iszero(r.σ)
        verbose && println("Initial scale is 0.0, no iterations performed.")
        return m
    end

    # Compute initial deviance
    devold = deviance(m)
    absdev = abs(devold)
    dev = devold
    Δdev = 0

    verbose && println("initial deviance: $(@sprintf("%.4g", devold))")
    for i in 1:maxiter
        f = 1.0 # line search factor
        # local dev
        absdev = abs(devold)

        # Compute the change to β, update μ and compute deviance
        dev = try
            deviance(setη!(m; updatescale=false))
        catch e
            isa(e, DomainError) ? Inf : rethrow(e)
        end

        # Assert the deviance is positive (up to rounding error)
        @assert dev > -atol

        verbose && println(
            "deviance at step $i: $(@sprintf("%.4g", dev)), crit=$((devold - dev)/absdev)",
        )

        # Line search
        ## If the deviance isn't declining then half the step size
        ## The rtol*abs(devold) term is to avoid failure when deviance
        ## is unchanged except for rounding errors.
        while dev > devold + rtol * absdev
            f /= 2
            f > minstepfac ||
                error("linesearch failed at iteration $(i) with beta0 = $(p.beta0)")

            dev = try
                # Update μ and compute deviance with new f. Do not recompute ∇β
                deviance(setη!(m, f))
            catch e
                isa(e, DomainError) ? Inf : rethrow(e)
            end
        end
        installbeta!(p, f)

        # Test for convergence
        Δdev = (devold - dev)
        verbose && println("Iteration: $i, deviance: $dev, Δdev: $(Δdev)")
        tol = max(rtol * absdev, atol)
        if -tol < Δdev < tol || dev < atol
            cvg = true
            break
        end
        @assert isfinite(dev)
        devold = dev
    end
    cvg || throw(ConvergenceException(maxiter))
    m
end

"""
    pirls_Sestimate!(m::RobustLinearModel{T}; verbose::Bool=false, maxiter::Integer=30,
           minstepfac::Real=1e-3, atol::Real=1e-6, rtol::Real=1e-6,
           beta0::AbstractVector=T[], sigma0::Union{Nothing, T}=nothing)

(Penalized) Iteratively Reweighted Least Square procedure for S-estimation.
The Penalized aspect is not implemented (yet).
"""
function pirls_Sestimate!(
    m::RobustLinearModel{T};
    verbose::Bool=false,
    maxiter::Integer=30,
    minstepfac::Real=1e-3,
    atol::Real=1e-6,
    rtol::Real=1e-5,
    miniter::Int=2,
    beta0::AbstractVector=[],
    sigma0::Union{Nothing,T}=nothing,
) where {T<:AbstractFloat}

    # Check arguments
    maxiter >= 1 || throw(ArgumentError("maxiter must be positive"))
    0 < minstepfac < 1 || throw(ArgumentError("minstepfac must be in (0, 1)"))

    # Extract fields and set convergence flag
    cvg, p, r = false, m.pred, m.resp

    ## Initialize σ, default to largest value
    maxσ = -(-(extrema(r.y)...)) / 2
    if !isnothing(sigma0)
        r.σ = sigma0
    else
        r.σ = maxσ
    end

    ## Initialize β or set it to the provided values
    setβ0!(m, beta0)

    # Initialize μ and compute residuals
    setinitη!(m)

    # Compute initial scale
    sigold =
        scale(setη!(m; updatescale=true, verbose=verbose, sigma0=sigma0, fallback=maxσ))
    installbeta!(p, 1)
    r.σ = sigold

    verbose && println("initial scale: $(@sprintf("%.4g", sigold))")
    for i in 1:maxiter
        f = 1.0 # line search factor
        local sig

        # Compute the change to β, update μ and compute deviance
        sig =
            scale(setη!(m; updatescale=true, verbose=verbose, sigma0=sigold, fallback=maxσ))

        # Assert the deviance is positive (up to rounding error)
        @assert sig > -atol

        verbose && println(
            "scale at step $i: $(@sprintf("%.4g", sig)), crit=$((sigold - sig)/sigold)",
        )

        # Line search
        ## If the scale isn't declining then half the step size
        ## The rtol*abs(sigold) term is to avoid failure when scale
        ## is unchanged except for rounding errors.
        while sig > sigold * (1 + rtol)
            f /= 2
            if f <= minstepfac
                if i <= miniter
                    sigold = maxσ
                    r.σ = sigold
                    verbose && println(
                        "linesearch failed at early iteration $(i), set scale to maximum value: $(sigold)",
                    )
                else
                    error("linesearch failed at iteration $(i) with beta0 = $(p.beta0)")
                end
            end
            # Update μ and compute deviance with new f. Do not recompute ∇β
            sig = scale(
                setη!(m; updatescale=true, verbose=verbose, sigma0=sigold, fallback=maxσ),
            )
        end
        installbeta!(p, f)
        r.σ = sig

        # Test for convergence
        Δsig = (sigold - sig)
        verbose && println("Iteration: $i, scale: $sig, Δsig: $(Δsig)")
        tol = max(rtol * sigold, atol)
        if -tol < Δsig < tol || sig < atol
            cvg = true
            break
        end
        @assert isfinite(sig) && !iszero(sig)
        sigold = sig
    end
    cvg || throw(ConvergenceException(maxiter))
    m
end

"""
    pirls_τestimate!(m::RobustLinearModel{T}; verbose::Bool=false, maxiter::Integer=30,
           minstepfac::Real=1e-3, atol::Real=1e-6, rtol::Real=1e-6,
           beta0::AbstractVector=T[], sigma0::Union{Nothing, T}=nothing)

(Penalized) Iteratively Reweighted Least Square procedure for τ-estimation.
The Penalized aspect is not implemented (yet).
"""
function pirls_τestimate!(
    m::RobustLinearModel{T};
    verbose::Bool=false,
    maxiter::Integer=30,
    minstepfac::Real=1e-3,
    atol::Real=1e-6,
    rtol::Real=1e-5,
    miniter::Int=2,
    beta0::AbstractVector=[],
    sigma0::Union{Nothing,T}=nothing,
) where {T<:AbstractFloat}

    # Check arguments
    maxiter >= 1 || throw(ArgumentError("maxiter must be positive"))
    0 < minstepfac < 1 || throw(ArgumentError("minstepfac must be in (0, 1)"))

    # Extract fields and set convergence flag
    cvg, p, r = false, m.pred, m.resp

    ## Initialize σ, default to largest value
    maxσ = -(-(extrema(r.y)...)) / 2
    if !isnothing(sigma0)
        r.σ = sigma0
    else
        r.σ = maxσ
    end

    ## Initialize β or set it to the provided values
    setβ0!(m, beta0)

    # Initialize μ and compute residuals
    setinitη!(m)

    # Compute initial τ-scale
    tauold = tauscale(setη!(m; updatescale=true); verbose=verbose)
    installbeta!(p, 1)

    verbose && println("initial τ-scale: $(@sprintf("%.4g", tauold))")
    for i in 1:maxiter
        f = 1.0 # line search factor
        local tau

        # Compute the change to β, update μ and compute deviance
        tau = tauscale(
            setη!(m; updatescale=true, verbose=verbose, fallback=maxσ);
            verbose=verbose,
        )

        # Assert the deviance is positive (up to rounding error)
        @assert tau > -atol

        verbose && println(
            "scale at step $i: $(@sprintf("%.4g", tau)), crit=$((tauold - tau)/tauold)",
        )

        # Line search
        ## If the scale isn't declining then half the step size
        ## The rtol*abs(sigold) term is to avoid failure when scale
        ## is unchanged except for rounding errors.
        while tau > tauold + rtol * tau
            f /= 2
            if f <= minstepfac
                if i <= miniter
                    tauold = maxσ
                    r.σ = tauold
                    verbose && println(
                        "linesearch failed at early iteration $(i), set scale to maximum value: $(tauold)",
                    )
                else
                    error("linesearch failed at iteration $(i) with beta0 = $(p.beta0)")
                end
            end

            # Update μ and compute deviance with new f. Do not recompute ∇β
            tau = tauscale(setη!(m; updatescale=true))
        end
        installbeta!(p, f)

        # Test for convergence
        Δtau = (tauold - tau)
        verbose && println("Iteration: $i, scale: $tau, Δsig: $(Δtau)")
        tol = max(rtol * tauold, atol)
        if -tol < Δtau < tol || tau < atol
            cvg = true
            break
        end
        @assert isfinite(tau) && !iszero(tau)
        tauold = tau
    end
    cvg || throw(ConvergenceException(maxiter))
    m
end

##########
###   Resampling
##########

"""
For S- and τ-Estimators, compute the minimum number of subsamples to draw
to ensure that with probability 1-α, at least one of the subsample
is free of outlier, given that the ratio of outlier/clean data is ε.
The number of data point per subsample is p, that should be at least
equal to the degree of freedom.
"""
function resampling_minN(p::Int, α::Real=0.05, ε::Real=0.5)
    ceil(Int, abs(log(α) / log(1 - (1 - ε)^p)))
end

function resampling_initialcoef(m, inds)
    # Get the subsampled model matrix, response and weights
    Xi = modelmatrix(m)[inds, :]
    yi = response(m)[inds]
    wi = Vector(weights(m)[inds])

    # Fit with OLS
    coef(lm(Xi, yi; wts=wi))
end

"""
    best_from_resampling(m::RobustLinearModel, ::Type{E}; kwargs...) where {E<:Union{SEstimator, MMEstimator, TauEstimator}}

Return the best scale σ0 and coefficients β0 from resampling of the S- or τ-Estimate.
"""
function resampling_best_estimate(
    m::RobustLinearModel,
    ::Type{E};
    propoutliers::Real=0.5,
    Nsamples::Union{Nothing,Int}=nothing,
    Nsubsamples::Int=10,
    Npoints::Union{Nothing,Int}=nothing,
    Nsteps_β::Int=2,
    Nsteps_σ::Int=1,
    verbose::Bool=false,
    rng::AbstractRNG=GLOBAL_RNG,
) where {E<:Union{SEstimator,MMEstimator,TauEstimator}}

    ## TODO: implement something similar to DetS (not sure it could apply)
    ## Hubert2015 - The DetS and DetMM estimators for multivariate location and scatter
    ## (https://www.sciencedirect.com/science/article/abs/pii/S0167947314002175)

    if isnothing(Nsamples)
        Nsamples = resampling_minN(dof(m), 0.05, propoutliers)
    end
    if isnothing(Npoints)
        Npoints = dof(m)
    end
    Nsubsamples = min(Nsubsamples, Nsamples)

    verbose && println("Start $(Nsamples) subsamples...")
    σis = zeros(Nsamples)
    βis = zeros(dof(m), Nsamples)
    for i in 1:Nsamples
        # TODO: to parallelize, make a deepcopy of m
        inds = sample(rng, 1:nobs(m), Npoints; replace=false, ordered=false)
        # Find OLS fit of the subsample
        βi = resampling_initialcoef(m, inds)
        verbose && println("Sample $(i)/$(Nsamples): β0 = $(βi)")

        ## Initialize β or set it to the provided values
        setβ0!(m, βi)
        # Initialize μ and compute residuals
        setinitη!(m)
        # Initialize σ as mad(residuals)
        setinitσ!(m)

        σi = 0
        for k in 1:Nsteps_β
            setη!(
                m;
                updatescale=true,
                verbose=verbose,
                sigma0=:mad,
                nmax=Nsteps_σ,
                approx=true,
            )

            σi = if E <: TauEstimator
                tauscale(m)
            else # if E <: Union{SEstimator, MMEstimator}
                scale(m)
            end
            installbeta!(m.pred, 1)
        end
        σis[i] = σi
        βis[:, i] .= coef(m)
        verbose && println("Sample $(i)/$(Nsamples): β1=$(βis[:, i])\tσ1=$(σi)")
    end

    verbose && println("Sorted scales: $(sort(σis))")
    inds = sortperm(σis)[1:Nsubsamples]
    σls = σis[inds]
    βls = βis[:, inds]

    verbose && println("Keep best $(Nsubsamples) subsamples: $(inds)")
    for l in 1:Nsubsamples
        σl, βl = σls[l], βls[:, l]
        # TODO: to parallelize, make a deepcopy of m

        if E <: Union{SEstimator,MMEstimator}
            try
                pirls_Sestimate!(m; verbose=verbose, beta0=βl, sigma0=σl, miniter=3)
                σls[l] = scale(m)
            catch e
                # Didn't converge, set to infinite scale
                σls[l] = Inf
            end
        elseif E <: TauEstimator
            try
                pirls_τestimate!(m; verbose=verbose, beta0=βl, sigma0=σl)
                σls[l] = tauscale(m)
            catch e
                # Didn't converge, set to infinite scale
                σls[l] = Inf
            end
        else
            error("estimator $E not supported.")
        end
        # Update coefficients
        βls[:, l] .= coef(m)

        verbose && println("Subsample $(l)/$(Nsubsamples): β2=$(βls[:, l])\tσ2=$(σls[l])")
    end
    N = argmin(σls)
    ## TODO: for τ-Estimate, the returned scale is τ not σ
    verbose && println("Best subsample: β=$(βls[:, N])\tσ=$(σls[N])")
    return σls[N], βls[:, N]
end


text:
# Copyright Henry S. Grasshorn Gebhardt
# This is a module to calculate a Spline.

# TODO:
# - Derivatives and integrals are not correct in extrapolation region.

module Splines

export Spline1D, Spline1Dtrans, Spline1Dloglog
export evaluate
export derivative
export integrate

#using TimerOutputs

@enum ExtrapolationType zero=0 boundary=1 linear=2 powerlaw=3 lastpoly=4

mutable struct Spline1D{k,T}
    x::Array{T,1}
    y::Array{T,1}
    ypp::Array{T,1}
    ilast::Int  # for caching
    betainv::Array{T,1} # workspace
    extrapolation::ExtrapolationType
end

# To make fused loops work
if VERSION < v"0.7.0"
    import Base.size, Base.getindex
    size(s::Spline1D) = (1,)
    getindex(s::Spline1D, i) = s
else
    import Base.length, Base.iterate
    length(s::Spline1D) = 1
    iterate(s::Spline1D) = s, nothing
    iterate(s::Spline1D, x) = nothing
end

################# All splines
function Spline1D(x, y; k=3, extrapolation=boundary)
    # checks
    if k < 0 || k == 2 || k > 3
        error("k=$k is not supported. :(")
    end
    if k == 0 && length(x) != length(y) + 1
        error("k=0 requires length(x) == length(y) + 1 (no OnCell() support)")
    elseif k != 0 && length(x) != length(y)
        error("Length of 'x' ($(length(x))) and 'y' ($(length(y))) must match.")
    end
    length(x) >= k+1 || error("Need at least k+1 points (length(x)=$(length(x))).")

    # do it, man
    T = promote_type(eltype(x), eltype(y))
    xin = convert(Array{T,1}, deepcopy(x))
    yin = convert(Array{T,1}, deepcopy(y))
    spline = Spline1D{k,T}(xin, yin, [], 1, [], extrapolation)
    #@show spline.extrapolation
    return Spline1D(xin, yin, spline)
end

Spline1D(; k=3, T=Float64, extrapolation=boundary) = Spline1D{k,T}([], [], [], 1, [], extrapolation)

function evaluate(spline::Spline1D, x)
    N = length(spline.x)
    i = findlargestsmaller(spline.x, x, spline)
    (i < 1 || i >= N) && return evaluate_past_boundary(spline, x)
    return eval_spline_piece_unsafe(spline, i, x)
end

(s::Spline1D)(x) = evaluate(s, x)

function evaluate_past_boundary(spl::Spline1D{k,T}, x::Number) where {k,T}
    idx = (x < spl.x[1]) ? 1 : length(spl.y)
    #@show x,k,spl.x[1],spl.x[end],length(spl.x),length(spl.y),idx
    if spl.extrapolation == zero
        return T(0)
    elseif spl.extrapolation == boundary
        return spl.y[idx]
    elseif spl.extrapolation == linear
        x0 = spl.x[idx]
        y0 = spl.y[idx]
        yp0 = derivative(spl, x0)
        return yp0 * (x - x0) + y0
    elseif spl.extrapolation == powerlaw
        x0 = spl.x[idx]
        y0 = spl.y[idx]
        yp0 = derivative(spl, x0)
        p = x0 * yp0 / y0
        #@show "powerlaw",x,x0,yp0,p
        return y0 * (x / x0)^p
    elseif spl.extrapolation == lastpoly
        if idx == length(spl.y)
            idx -= 1
        end
        return eval_spline_piece_unsafe(spl, idx, x)
    else
        error("Unkown ExtrapolationType '$(spl.extrapolation)'")
    end
end

#################### 0th-order spline ##################
function Spline1D(x, y, spl::Spline1D{0})
	spl.x = x
	spl.y = y
	spl.ilast = 1
        return spl
end

function eval_spline_piece_unsafe(spline::Spline1D{0}, i::Integer, x::Number)
	@inbounds return spline.y[i]
end

function derivative(spline::Spline1D{0,T}, x::Number) where {T}
        xx = spline.x
        yy = spline.y
	N = length(xx)
	i = findlargestsmaller(xx, x, spline)
	if i <= 1 || i >= N
            return T(0)
	end
        if x == xx[i]
            if yy[i-1] == yy[i]
                return T(0)
            end
            return (yy[i-1] < yy[i]) ? T(Inf) : T(-Inf)
        elseif x == xx[i+1]
            if yy[i] == yy[i+1]
                return T(0)
            end
            return (yy[i] < yy[i+1]) ? T(Inf) : T(-Inf)
        end
        return T(0)
end

#################### 1st-order spline ##################
function Spline1D(x, y, spl::Spline1D{1})
    spl.x = x
    spl.y = y
    spl.ilast = 1
    return spl
end

function eval_spline_piece_unsafe(spline::Spline1D{1}, i::Integer, x::Number)
        xx = spline.x
        yy = spline.y
        @inbounds return yy[i] + (x - xx[i]) / (xx[i+1] - xx[i]) * (yy[i+1] - yy[i])
end

function derivative(spline::Spline1D{1,T}, x::Number) where {T}
        xx = spline.x
        yy = spline.y
	N = length(xx)
	i = findlargestsmaller(xx, x, spline)
	if i < 1 || i >= N
            return T(0)
	end
        @inbounds return (yy[i+1] - yy[i]) / (xx[i+1] - xx[i])
end

#################### 3rd-order spline ##################
Spline1D(x, y, spl::Spline1D{3}) = Spline1Dk3(x, y, spl.betainv, spl.ypp, spl)

function Ai(x, y, i)
	6(y[i+1] - y[i]) / (x[i+1] - x[i]) - 6(y[i] - y[i-1]) / (x[i] - x[i-1])
end

function yderiv_k2(x, y, i, j, k)
	((y[k]-y[i])*(x[j]-x[i])^2 - (y[j]-y[i])*(x[k]-x[i])^2) / (
		(x[k]-x[i]) * (x[j]-x[k]) * (x[j] - x[i]))
end

function yderiv_k3(x, y, i, j, k, l)
	# Note: I don't trust this mathematics, yet
	x1 = x[i]
	x2 = x[j]
	x3 = x[k]
	x4 = x[l]
	y1 = y[i]
	y2 = y[j]
	y3 = y[k]
	y4 = y[l]
	hinv = 1 / (x4 - x1)
	t2 = (x2 - x1) * hinv
	t3 = (x3 - x1) * hinv
	yp1 = hinv * t3 / (t2*t3*(1-t2^2) - t2^2*(1-t2*(1+t3))) * (
		y2 - y1 - (y4-y1)*t2^3
		- t2^2*(1-t2) / (t3^2*(1-t3)) * (y3 - y1 - (y4-y1)*t3^3))
	return yp1
end

# for k = 3
function Spline1Dk3(x, y, betainv, ypp, spline::Spline1D{3})
	N = length(x)

	resize!(ypp, N)
	if length(betainv) < N-1
		resize!(betainv, N-1)
	end

	u = ypp  # alias for clarity

	##ypp[1] = 0.0
	#ypp[1] = Ai(x, y, 2) / (3 * (x[3] - x[1]))
	#A2adj = - (x[2] - x[1]) * ypp[1]
	#a22adj = 0.0

	##ypp[N] = 0.0
	#ypp[N] = Ai(x, y, N-1) / (3  * (x[N] - x[N-2]))
	#ANadj = - (x[N] - x[N-1]) * ypp[N]
	#aNNadj = 0.0

	yp1 = yderiv_k2(x, y, 1, 2, 3)
	#yp1 = yderiv_k3(x, y, 1, 2, 3, 4)
	A2adj = 3 * yp1 - 3 * (y[2] - y[1]) / (x[2] - x[1])
	a22adj = - (x[2] - x[1]) / 2

	ypN = yderiv_k2(x, y, N, N-1, N-2)
	#ypN = yderiv_k3(x, y, N, N-1, N-2, N-3)
	ANadj = 3 * (y[N] - y[N-1]) / (x[N] - x[N-1]) - 3 * ypN
	aNNadj = - (x[N] - x[N-1]) / 2

	# forward substitution
	betainv[2] = 1 / (2 * (x[3] - x[1]) + a22adj)
	u[2] = Ai(x, y, 2) + A2adj
	for i=3:N-2
		xixim1betainvim1 = (x[i] - x[i-1]) * betainv[i-1]
		betainv[i] = 1 / (2(x[i+1] - x[i-1]) - (x[i] - x[i-1]) * xixim1betainvim1)
		u[i] = Ai(x, y, i) - xixim1betainvim1 * u[i-1]
	end
	i = N-1
	xixim1betainvim1 = (x[i] - x[i-1]) * betainv[i-1]
	betainv[i] = 1 / (2(x[i+1] - x[i-1]) + aNNadj - (x[i] - x[i-1]) * xixim1betainvim1)
	u[i] = Ai(x, y, i) + ANadj - xixim1betainvim1 * u[i-1]

	# backward substitution
	ypp[N-1] = betainv[N-1] * u[N-1]
	for i=N-2:-1:2
		ypp[i] = betainv[i] * (u[i] - (x[i+1] - x[i]) * ypp[i+1])
	end

	ypp[1] = (a22adj * ypp[2] - A2adj) / (x[2] - x[1])
	ypp[N] = (aNNadj * ypp[N-1] - ANadj) / (x[N] - x[N-1])

        spline.x = x
	spline.y = y
	spline.ypp = ypp
	spline.ilast = 1
	spline.betainv = betainv
	return spline
end

function findlargestsmaller_bisect(a, x::Number, nmin, nmax)::Int64
	nmid = div(nmin + nmax, 2)
	while nmid != nmin
		@inbounds if x >= a[nmid]
			nmin = nmid
		else
			nmax = nmid
		end
		nmid = div(nmin + nmax, 2)
	end
	return nmin
end

#global ihist = zeros(Int, 2001)

function findlargestsmaller(a, x, ilast::Int)
	N = length(a)

	@inbounds if x < a[1]
		return 0
	elseif x > a[end]
		return N + 1
	end

        nmin = (1 <= ilast <= N) ? ilast : 1
	@inbounds while x < a[nmin]
		nmin = div(nmin, 2)
	end

        # handle most common cases
        @inbounds if nmin < N && x < a[nmin + 1]  # i == ilast == nmin
            #global ihist[0 + 1000 + 100] += 1
            return nmin
        elseif nmin < N - 1 && x < a[nmin + 2]  # i == ilast + 1
            #global ihist[1 + 1000 + 100] += 1
            return nmin + 1
        end

        nmax = nmin  # nmin <= N
        @inbounds while x > a[nmax]  # x > a[N] is already excluded
		nmax = div(nmax + N, 2) + 1
	end

        i = findlargestsmaller_bisect(a, x, nmin, nmax)
        #global ihist[i - ilast + 1000] += 1
        return i
end

function findlargestsmaller(a, x, spline)
    i = findlargestsmaller(a, x, spline.ilast)
    spline.ilast = i
    return i
end

function get_abcdh(spline, i, x)
	@inbounds x0 = spline.x[i]
	@inbounds x1 = spline.x[i+1]
	h = x1 - x0
	h2 = h^2
	@inbounds ypp0 = spline.ypp[i]
	@inbounds ypp1 = spline.ypp[i+1]

	@inbounds d = spline.y[i]
	b = h2 * ypp0 / 2
	a = h2 * (ypp1 - ypp0) / 6
	@inbounds c = spline.y[i+1] - a - b - d
	t = (x - x0) / h
	return a, b, c, d, t, h
end

function eval_spline_piece_unsafe(spline::Spline1D{3}, i::Integer, x::Number)
	a, b, c, d, t, = get_abcdh(spline, i, x)
	@evalpoly(t, d, c, b, a)
end

function derivative(spline::Spline1D{3,T}, x::Number) where {T}
	N = length(spline.x)
	i = findlargestsmaller(spline.x, x, spline)
	if i < 1
            return T(0)
	elseif i > N
            return T(0)
	end
	a, b, c, d, t, h = get_abcdh(spline, i, x)
	@evalpoly(t, c, 2b, 3a) / h
end

# integrate():
#   Integrate the spline 'spline' over the interval [xlo, xhi].
function integrate(spline::Spline1D{3,T}, xlo::Number, xhi::Number) where {T}
	a = findlargestsmaller(spline.x, xlo, spline)
	b = findlargestsmaller(spline.x, xhi, spline)
	N = length(spline.x)

        I = T(0)
        if spline.extrapolation == zero && a > N-1
            return T(0)
        elseif spline.extrapolation == boundary
            if a < 1
                    I += (spline.x[1] - xlo) * spline.y[1]
                    xlo = spline.x[1]
            elseif a > N-1
                    return (xhi - xlo) * spline.y[N]
            end
            if b > N-1
                    I += (xhi - spline.x[N]) * spline.y[N]
                    xhi = spline.x[N]
            end
        end
        (a < 1) && (a = 1)
        (b > N-1) && (b = N-1)

	x0 = spline.x[a]
	y0 = spline.y[a]
	ypp0 = spline.ypp[a]
	x1 = spline.x[a+1]
	y1 = spline.y[a+1]
	ypp1 = spline.ypp[a+1]
	h = x1 - x0
	h2 = h^2
	bi = h2 * ypp0 / 2
	ai = h2 * (ypp1 - ypp0) / 6
	ci = y1 - y0 - ai - bi
	t = (xlo - x0) / h
	I += - h * t * @evalpoly(t, y0, ci/2, bi/3, ai/4)
        #dI = h * (ai/4 + bi/3 + ci/2 + y0)
        I1 = T(0)
        I2 = T(0)
        dI1 = h * (y0 + y1)
        dI2 = h^3 * (ypp0 + ypp1)
	@inbounds for i=a+2:b+1
                #I += dI
                #I += dI1/2 - dI2/24
                I1 += dI1
                I2 += dI2

		# prepare the next:
		x0 = x1
		y0 = y1
		ypp0 = ypp1
		x1 = spline.x[i]
		y1 = spline.y[i]
		ypp1 = spline.ypp[i]

		h = x1 - x0
		#h2 = h^2
		#bi = h2 * ypp0 / 2
		#ai = h2 * (ypp1 - ypp0) / 6
		#ci = y1 - y0 - ai - bi
                #dI = h * (ai/4 + bi/3 + ci/2 + y0)

                #dI = h * ((y0 + y1) * (1/2) - h2 * (ypp0 + ypp1) * (1/24))
                dI1 = h * (y0 + y1)
                dI2 = h^3 * (ypp0 + ypp1)
	end
        I += I1/2 - I2/24
        h = x1 - x0
        h2 = h^2
        bi = h2 * ypp0 / 2
        ai = h2 * (ypp1 - ypp0) / 6
        ci = y1 - y0 - ai - bi
	t = (xhi - x0) / h
	I += h * t * @evalpoly(t, y0, ci/2, bi/3, ai/4)

	return I
end

######################## Spline1Dloglog
# This is a specialization of Spline1Dtrans. Reason: be faster.
struct Spline1Dloglog{Spline1DType}
	spline::Spline1DType
end

size(s::Spline1Dloglog) = (1,)
getindex(s::Spline1Dloglog, i) = s

# empty initialization
Spline1Dloglog(; k=3, T=Float64, extrapolation=boundary) = Spline1Dloglog(Spline1D(k=k, T=T, extrapolation=extrapolation))

# Spline1D(): Same as the other constructors, except that the inputs are
# transformed by the 'x_in_trans()' and 'y_in_trans()' functions.
# 'y_out_trans()' will typically be the inverse of 'y_in_trans()', as it is
# used to undo the transformation.
function Spline1Dloglog(x, y; k=3, extrapolation=boundary)
	xt = log.(x)
	yt = log.(y)
	spl = Spline1D(xt, yt; k=k, extrapolation=extrapolation)
        Spline1Dloglog(spl)
end

function Spline1D(x, y, spl::Spline1Dloglog)
    resize!(spl.spline.x, Int64(length(x)))
    resize!(spl.spline.y, Int64(length(y)))
    @. spl.spline.x = log(x)
    @. spl.spline.y = log(y)
    Spline1D(spl.spline.x, spl.spline.y, spl.spline)
    return spl
end

evaluate(s::Spline1Dloglog, x) = begin
    lnx = log(x)                  # ~20% of time
    #@show s.spline.extrapolation
    lnv = evaluate(s.spline, lnx) # ~30% of time
    v = exp(lnv)                  # ~50% of time
    return v
    #exp(evaluate(s.spline, log(x)))
end

(s::Spline1Dloglog)(x) = evaluate(s, x)

function derivative(s::Spline1Dloglog, x)
	xt = log(x)
	sval = evaluate(s.spline, xt)
	dyout = exp(sval)
	ds = derivative(s.spline, xt)
	dxt = inv(x)
	return dyout * ds * dxt
end

######################## Spline1Dtrans
struct Spline1Dtrans{Spline1DType}
	spline::Spline1DType
	x_in_trans::Function
	y_in_trans::Function
	y_out_trans::Function
	x_in_trans_deriv::Function
	y_out_trans_deriv::Function
end

size(s::Spline1Dtrans) = (1,)
getindex(s::Spline1Dtrans, i) = s

# Spline1D(): Same as the other constructors, except that the inputs are
# transformed by the 'x_in_trans()' and 'y_in_trans()' functions.
# 'y_out_trans()' will typically be the inverse of 'y_in_trans()', as it is
# used to undo the transformation.
function Spline1Dtrans(x, y, x_in_trans::Function, y_in_trans::Function, y_out_trans::Function,
                       x_in_trans_deriv::Function, y_out_trans_deriv::Function;
                       k=3, extrapolation=boundary)
	xt = x_in_trans.(x)
	yt = y_in_trans.(y)
	spl = Spline1D(xt, yt; k=k, extrapolation=extrapolation)
	Spline1Dtrans(spl, x_in_trans, y_in_trans, y_out_trans, x_in_trans_deriv, y_out_trans_deriv)
end

function find_inverse(fn)
	if fn == identity
		return identity
	elseif fn == log
		return exp
	elseif fn == exp
		return log
	else
		error("Inverse of function '$fn' not implemented here.")
		return exp
	end
end

function find_derivative(fn)
	if fn == identity
		return one
	elseif fn == log
		return inv
	elseif fn == exp
		return exp
	else
		error("Derivative of function '$fn' not implemented here.")
		return exp
	end
end

function Spline1Dtrans(x, y, x_in_trans::Function, y_in_trans::Function;
                       k=3, extrapolation=boundary)
	x_in_trans_deriv = find_derivative(x_in_trans)
	y_out_trans = find_inverse(y_in_trans)
	y_out_trans_deriv = find_derivative(y_out_trans)
	Spline1Dtrans(x, y, x_in_trans, y_in_trans, y_out_trans,
		      x_in_trans_deriv, y_out_trans_deriv; k=k, extrapolation=extrapolation)
end

function Spline1D(x, y, spl::Spline1Dtrans)
    resize!(spl.spline.x, length(x))
    resize!(spl.spline.y, length(y))
    @. spl.spline.x = spl.x_in_trans(x)
    @. spl.spline.y = spl.y_in_trans(y)
    Spline1D(spl.spline.x, spl.spline.y, spl.spline)
    return spl
end

evaluate(s::Spline1Dtrans, x) = s.y_out_trans(evaluate(s.spline, s.x_in_trans(x)))

(s::Spline1Dtrans)(x) = evaluate(s, x)

function derivative(s::Spline1Dtrans, x)
	xt = s.x_in_trans(x)
	sval = evaluate(s.spline, xt)
	dyout = s.y_out_trans(sval)
	ds = derivative(s.spline, xt)
	dxt = s.x_in_trans_deriv(x)
	return dyout * ds * dxt
end

end

# vim: set sw=4 et sts=4 :


text:
abstract type Solution{dType, tType, N} end

nsave(sol::Solution) = error("nsave() not implemented for ", typeof(sol))
ntime(sol::Solution) = error("ntime() not implemented for ", typeof(sol))
nsamples(sol::Solution) = error("nsamples() not implemented for ", typeof(sol))

counter(sol::Solution) = error("counter() not implemented for ", typeof(sol))
offset(sol::Solution) = error("offset() not implemented for ", typeof(sol))
lastentry(sol::Solution) = error("lastentry() not implemented for ", typeof(sol))
timesteps(sol::Solution) = error("timesteps() not implemented for ", typeof(sol))
eachtimestep(sol::Solution) = error("eachtimestep() not implemented for ", typeof(sol))

eachsample(sol::Solution) = 1:nsamples(sol)

Base.write(::IOStream, sol::Solution) = error("Base.write(::IOStream, ::Solution) not implemented for ", typeof(sol))

function test_interface(sol::Solution)

    @test_nowarn nsave(sol)
    @test_nowarn ntime(sol)
    @test_nowarn nsamples(sol)

    @test_nowarn counter(sol)
    @test_nowarn offset(sol)
    @test_nowarn lastentry(sol)
    @test_nowarn timesteps(sol)

    @test_nowarn eachtimestep(sol)
    @test_nowarn eachsample(sol)

end


text:
using Soss

μdist = @model ν begin
    s ~ Gamma(ν , ν)
    z ~ Normal()
    return sqrt(s)*z
end

σdist = @model begin
    x ~ Normal()
    return abs(x)
end

m = @model begin
    μ ~ μdist(ν=1.0)
    σ ~ σdist()
    x ~ Normal(μ,σ) |> iid(10)
    return x
end

simulate(μdist(ν=1.0))

simulate(σdist())

simulate(m())

rand(m())

x = randn(2)

xform(m() | (x=x, μ = (z = 1.0,)))

dynamicHMC(m() | (x=x, μ = (z = 1.0,)))

_data = (x=x, μ = (z = 1.0,))

# function show(io::IO, t::NamedTuple)
#     n = nfields(t)
#     for i = 1:n
#         # if field types aren't concrete, show full type
#         if typeof(getfield(t, i)) !== fieldtype(typeof(t), i)
#             show(io, typeof(t))
#             print(io, "(")
#             show(io, Tuple(t))
#             print(io, ")")
#             return
#         end
#     end
#     if n == 0
#         print(io, "NamedTuple()")
#     else
#         typeinfo = get(io, :typeinfo, Any)
#         print(io, "(")
#         for i = 1:n
#             print(io, fieldname(typeof(t),i), " = ")
#             show(IOContext(io, :typeinfo =>
#                            t isa typeinfo <: NamedTuple ? fieldtype(typeinfo, i) : Any),
#                  getfield(t, i))
#             if n == 1
#                 print(io, ",")
#             elseif i < n
#                 print(io, ", ")
#             end
#         end
#         print(io, ")")
#     end
# end


text:
module DnD

Pkg.activate(".")
using CSV
using DataFrames

creatures = CSV.read("data_d20datascience\\src\\creatures.csv")
fighter   = CSV.read("data_d20datascience\\src\\fighter.csv")
wizard    = CSV.read("data_d20datascience\\src\\wizard.csv")

teamA = creatures[Integer.(round.(rand(1)*300 .+ (1:4))),:]
teamB = creatures[Integer.(round.(rand(1)*300 .+ (1:4))),:]

# function combatsim(teamA, teamB)

  teamA.team = 1
  teamB.team = 2
  teams = vcat(teamA, teamB)

  println("Team 1\n  ", join(teamA.name, "\n  "))
  println("Team 2\n  ", join(teamB.name, "\n  "))

  teams.currentHP = copy(teams.hp)
  teams.order = rand(size(teams,1)) .* 20 .+ 1 + teams.init
  sort!(teams, :order, rev = true)

  team1count = size(teamA, 1)
  team2count = size(teamB, 1)

  currentround = 0

  while ((team1count > 0) & (team2count > 0))
    global currentround
    println(currentround)
    currentround += 1
    # Loop through each creature
    for i in 1:size(teams, 1)
        # If current creature HP is 0 skip
      (teams.currentHP[i] <= 0) && continue

      currentteam =  teams.team[i]
      otherteam = teams.team .!== currentteam

      for j in 1:teams.att_num[i]
        global team1count, team2count

        potential = findall(otherteam .& (teams.currentHP .> 0) .== 1)

        target = potential[1]

        attack_roll  = ceil.(rand(1)*20)[1] |> Integer
        attack_value = attack_roll + teams.hit_mod[i]

        hit = (teams.ac[target] <= attack_value)

        (attack_roll == 20) && (hit == true)
        (attack_roll == 1) && (hit == false)

        hit && (teams.currentHP[target] -= teams.dam_avg[i])

        team1count = sum(teams.currentHP[teams.team .== 1, :] .> 0)
        team2count = sum(teams.currentHP[teams.team .== 2, :] .> 0)

        ((team1count==0) | (team2count==0)) && break

      end

      ((team1count==0) | (team2count==0)) && break
    end

    (team1count==0) && (println("Team 2 wins!"); break)
    (team2count==0) && (println("Team 1 wins!"); break)

    (currentround > 10) && break

  end

# end

end # module


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   ModiaMath.Frames (ModiaMath/Frames/_module.jl)
#
import ..StaticArrays: SVector, @SMatrix

"""
    const ModiaMath.Vector3D = SVector{3,Float64}

Type of a vector in 3D space (e.g. position vector of the origin of a frame)
"""
const Vector3D = SVector{3,Float64}

"""
    const ModiaMath.ZeroVector3D = Vector3D(0.0, 0.0, 0.0)

Constant of a Vector3D where all elements are zero
"""
const ZeroVector3D = Vector3D(0.0, 0.0, 0.0)

"""
    M = ModiaMath.skew(e::AbstractVector)

Return the skew symmetric matrix `M::SMatrix{3,3,Float64,9}` of vector `e` (`length(e) = 3`)
"""
skew(e::AbstractVector) = @SMatrix([  0.0   -e[3]    e[2];
                                      e[3]   0.0    -e[1];
                                     -e[2]   e[1]    0.0 ])


text:
"""
    KLocalOperatorTensor

A KLocalOperator representing the sum of several KLocalOperator-s. Internally,
the sum is stored as a vector of local operators acting on some sites.
"""
struct KLocalOperatorTensor{H,T,O1,O2} <: AbsLinearSuperOperator{T}
    hilb::H
    sites::T

    # list of sites in this sum
    op_l::O1
    op_r::O2
end

"""
    Builds the tensor product of two operators
"""
function KLocalOperatorTensor(op_l, op_r)
    if op_l isa KLocalIdentity
        st = (Int[],       sites(op_r))
        hilb = basis(op_r)
    elseif op_r isa KLocalIdentity
        st = (sites(op_l), Int[])
        hilb = basis(op_l)
    else
        @assert basis(op_l) == basis(op_r)
        hilb = basis(op_l)
        st = (sites(op_l), sites(op_r))
    end

    KLocalOperatorTensor(SuperOpSpace(hilb), st, op_l, op_r)
end

function KLocalOperatorTensor(op_l::KLocalOperatorSum, op_r::KLocalIdentity)
    ops = [KLocalOperatorTensor(op, op_r) for op=operators(op_l)]
    return sum(ops)
end

function KLocalOperatorTensor(op_l::KLocalIdentity, op_r::KLocalOperatorSum)
    ops = [KLocalOperatorTensor(op_l, op) for op=operators(op_r)]
    return sum(ops)
end

KLocalOperatorSum(op::KLocalOperatorTensor) = KLocalOperatorSum(basis(op), [sites(op)], [op])

function KLocalOperatorTensor(op_l::KLocalOperatorSum, op_r::KLocalOperatorSum)
    throw("error not impl")
end

QuantumOpticsBase.basis(op::KLocalOperatorTensor) = op.hilb
sites(op::KLocalOperatorTensor) = op.sites

## Accessors
operators(op::KLocalOperatorTensor) = (op,)

conn_type(op::KLocalOperatorTensor{H,T,O1,O2}) where {H,T,O1,O2}      = conn_type(op)
#conn_type(op::KLocalOperatorTensor{H,T,O1,Nothing}) where {H,T,O1}    = conn_type(op.op_l)
#conn_type(op::KLocalOperatorTensor{H,T,Nothing,O2}) where {H,T,O2}    = conn_type(op.op_r)

conn_type(::Type{KLocalOperatorTensor{H,T,O1,O2}}) where {H,T,O1,O2}  =
    OpConnectionTensor{conn_type(O1), conn_type(O2)}
conn_type(::Type{KLocalOperatorTensor{H,T,O1,Nothing}}) where{H,T,O1} = #conn_type(O1)
    OpConnectionTensor{conn_type(O1), eye_conn_type(O1)}
conn_type(::Type{KLocalOperatorTensor{H,T,Nothing,O2}}) where{H,T,O2} = #conn_type(O2)
    OpConnectionTensor{eye_conn_type(O2), conn_type(O2)}
conn_type(::Type{KLocalOperatorTensor{H,T,O1,<:KLocalIdentity}}) where{H,T,O1} = #conn_type(O1)
    OpConnectionTensor{conn_type(O1), eye_conn_type(O1)}
conn_type(::Type{KLocalOperatorTensor{H,T,<:KLocalIdentity,O2}}) where{H,T,O2} = #conn_type(O2)
    OpConnectionTensor{eye_conn_type(O2), conn_type(O2)}

duplicate(::Nothing) = nothing

function duplicate(op::KLocalOperatorTensor)
    KLocalOperatorTensor(duplicate(op.op_l), duplicate(op.op_r))
end

function _row_valdiff!(conn::AbsOpConnection, op::KLocalOperatorTensor, v::ADoubleState)
    op_r = op.op_r
    op_l = op.op_l
    if op_r isa KLocalIdentity
        r_r = local_index(row(v), basis(op_l), sites(op_l))
        append!(conn, (op_l.op_conns[r_r], nothing))
    elseif op_l isa KLocalIdentity
        r_c = local_index(col(v), basis(op_r), sites(op_r))
        append!(conn, (nothing, op_r.op_conns[r_c]))
    else
        r_r = local_index(row(v), basis(op_l), sites(op_l))
        r_c = local_index(col(v), basis(op_r), sites(op_r))
        append!(conn, (op_l.op_conns[r_r], op_r.op_conns[r_c]))
    end
    return conn
end

function map_connections(fun::Function, op::KLocalOperatorTensor, v::ADoubleState)
    op_r = op.op_r
    op_l = op.op_l
    hilb_ph = physical(basis(op))
    if op_r isa KLocalIdentity
        r = local_index(row(v), hilb_ph, sites(op_l))

        for (mel,changes)=op_l.op_conns[r]
            #fun(mel, 1.0, changes, nothing, v)
            fun(mel, (changes, nothing), v)
        end
    elseif op_l isa KLocalIdentity
        r = local_index(col(v), hilb_ph, sites(op_r))

        for (mel,changes)=op_r.op_conns[r]
            #fun(1.0, mel, nothing, changes, v)
            fun(mel, (nothing, changes), v)
        end
    else
        r_r = local_index(row(v), hilb_ph, sites(op_l))
        r_c = local_index(col(v), hilb_ph, sites(op_r))

        for (mel_r, changes_r)=op_l.op_conns[r_r]
            for (mel_c, changes_c)=op_r.op_conns[r_c]
                #fun(mel_r, mel_c, changes_r, changes_c, v)
                fun(mel_r*mel_c, (changes_r, changes_c), v)
            end
        end
    end
    return nothing
end

function accumulate_connections!(acc::AbstractAccumulator, op::KLocalOperatorTensor, v::ADoubleState)
    op_l = op.op_l; op_r = op.op_r

    if op.op_r isa KLocalIdentity
        r = local_index(row(v), basis(op_l), sites(op_l))

        for (mel,changes)=op_l.op_conns[r]
            #fun(mel, 1.0, changes, nothing, v)
            acc(mel, changes, nothing, v)
        end
    elseif op.op_l isa KLocalIdentity
        r = local_index(col(v), basis(op_r), sites(op_r))

        for (mel,changes)=op_r.op_conns[r]
            #fun(1.0, mel, nothing, changes, v)
            acc(mel, nothing, changes, v)
        end
    else
        r_r = local_index(row(v), basis(op_l), sites(op_l))
        r_c = local_index(col(v), basis(op_r), sites(op_r))

        for (mel_r, changes_r)=op_l.op_conns[r_r]
            for (mel_c, changes_c)=op_r.op_conns[r_c]
                acc(mel_r*mel_c, changes_r, changes_c, v)
            end
        end
    end
    return acc
end

_add_samesite(op_l::KLocalOperatorTensor, op_r::KLocalOperatorTensor) = _add_samesite!(duplicate(op_l), op_r)

_add_samesite!(::Nothing, ::Nothing) = nothing
function _add_samesite!(op_l::KLocalOperatorTensor, op_r::KLocalOperatorTensor)
    @assert op_l.sites == op_r.sites

    _add_samesite!(op_l.op_l, op_r.op_l)
    _add_samesite!(op_l.op_r, op_r.op_r)
    return op_l
end

function Base.conj!(op::KLocalOperatorTensor)
    !isnothing(op.op_l) && conj!(op.op_l)
    !isnothing(op.op_r) && conj!(op.op_r)
    return op
end

function Base.transpose(op::KLocalOperatorTensor)
    KLocalOperatorTensor(op.op_r, op.op_l)

end

Base.conj(op::KLocalOperatorTensor) = conj!(duplicate(op))

Base.adjoint(op::KLocalOperatorTensor) = conj(transpose(op))

Base.:+(op_l::KLocalOperatorTensor, op_r::KLocalOperatorTensor) = begin
    sites(op_l) == sites(op_r) && return _add_samesite(op_l, op_r)
    return KLocalOperatorSum(op_l) + op_r
end

Base.:*(a::Number, b::KLocalOperatorTensor) =
    _op_alpha_prod(b,a)

function _op_alpha_prod(op::KLocalOperatorTensor, a::Number)
    if op.op_l isa KLocalIdentity
        return KLocalOperatorTensor(op.op_l, a*op.op_r)
    elseif op.op_r isa KLocalIdentity
        return KLocalOperatorTensor(a*op.op_l, op.op_r)
    else
        a = sqrt(a)
        return KLocalOperatorTensor(a*op.op_l, a*op.op_r)
    end
end

Base.isapprox(l::KLocalOperatorTensor, r::KLocalOperatorTensor; kwargs...) =
    isapprox(l.op_l, r.op_l; kwargs...) && isapprox(l.op_r, r.op_r; kwargs...)

Base.:(==)(l::KLocalOperatorTensor, r::KLocalOperatorTensor) =
    l.op_l == r.op_l && l.op_r == r.op_r

Base.eltype(::T) where {T<:KLocalOperatorTensor} = eltype(T)
Base.eltype(T::Type{<:KLocalOperatorTensor{A,B,C}}) where {A,B,C<:KLocalOperator} =
    eltype(C)
Base.eltype(T::Type{<:KLocalOperatorTensor{A,B,C,D}}) where {A,B,C,D<:KLocalOperator} =
    eltype(D)
Base.eltype(T::Type{<:KLocalOperatorTensor{A,B,C,D}}) where {A,B,C<:KLocalOperator,D<:KLocalOperator} =
    eltype(C)


text:
using Test
using Dates: now
using ExtensibleScheduler
using ExtensibleScheduler: MemoryJobStore, Job, get_job_id

@testset "JobStores" begin

    function return_args(id)
        "Returned from return_args $id"
    end

    @testset "MemoryJobStore" begin
        jobstore = MemoryJobStore()
        @test isempty(jobstore)

        action1 = Action(return_args, 1)
        @test run(action1) == "Returned from return_args 1"
        trigger1 = Trigger(DateTime(2010, 1, 1))
        id1 = get_job_id(jobstore)
        name1 = "Textual description of job1"
        priority1 = 0
        dt_created1 = now()
        dt_updated1 = dt_created1
        dt_next_fire1 = DateTime(0)
        n_triggered1 = 0
        config1 = JobConfig()
        job1 = Job(id1, action1, trigger1, name1, priority1, dt_created1, dt_updated1, dt_next_fire1, n_triggered1, config1)
        push!(jobstore, job1)
        @test length(jobstore) == 1

        action2 = Action(return_args, 2)
        @test run(action2) == "Returned from return_args 2"
        trigger2 = Trigger(DateTime(2010, 1, 1))
        id2 = get_job_id(jobstore)
        @test id2 != id1
        name2 = "Textual description of job2"
        priority2 = 0
        dt_created2 = now()
        dt_updated2 = dt_created2
        dt_next_fire2 = DateTime(0)
        n_triggered2 = 0
        config2 = JobConfig()
        job2 = Job(id2, action2, trigger2, name2, priority2, dt_created2, dt_updated2, dt_next_fire2, n_triggered2, config2)
        push!(jobstore, job2)
        @test length(jobstore) == 2
    end

end


text:
@with_kw type DynaMo2D <: DynaMo
  name::String
  bodies::Vector{Body2D} = Body2D[]
  springs::Vector{Spring2D} = Spring2D[]
  phases::Vector{Phase} = Phase[]
  contacts::Vector{Body2D} = Body2D[]
  J::Matrix{Sym} = Matrix{Sym}[]
  MM::Matrix{Sym} = Matrix{Sym}[]
  SpringForces::Vector{Sym} = Vector{Sym}[]
  PositionStates::Vector{Sym} = Vector{Sym}[]
  VelocityStates::Vector{Sym} = Vector{Sym}[]
  States::Vector{Sym} = Vector{Sym}[]
  dim::Int = 2
  angdim::Int = 1
  fulldim::Int = 3
end

function addBody(m::DynaMo2D,name::String;parent::Body=ground2D,
  connection::String=free)
  body = Body2D(name,parent,connection)
  push!(m.bodies,body)
  return body
end

function MaximalMassMatrix(m::DynaMo)
  N = length(m.bodies)*m.fulldim
  M = zeros(SymFloat,N,N)
  for (i,body) in enumerate(m.bodies)
    dex = m.dim*(i-1)+ 1
    M[dex,dex] = body.mass
    M[dex+1,dex+1] = body.mass
    M[dex+2,dex+2] = body.inertia
  end
  return M
end

function poseJacobian(m::DynaMo)
  J = zeros(SymFloat,3*length(m.bodies),m.DOF)
  for (i,body) in enumerate(m.bodies)
    J[m.fulldim*(i-1)+1:m.fulldim*i,:] = poseJacobian(body,m)
  end
  return J
end

function poseJacobian(b::Body,m::DynaMo)
  return vcat(translationVelocityJacobian(b,m),angularVelocityJacobian(b,m))
end

function poseJacobianDot(b::Body,m::DynaMo)
  poseJacobianDot = timederiv(poseJacobian(b,m))
end

function angularVelocityJacobian(m::DynaMo)
  J = zeros(SymFloat,m.angdim*length(m.bodies),m.DOF)

  for (i,body) in enumerate(m.bodies)
    J[m.angdim*(i-1)+1:m.angdim*i,:] = angularVelocityJacobian(body,m)
  end
  return J
end

function angularVelocityJacobian(b::Body,m::DynaMo)
  return linEq(b.angularvelocity,b.connection.u)
end

function translationalVelocityJacobian(m::DynaMo)
  J = zeros(SymFloat,m.dim*length(m.bodies),m.DOF)
  for (i,body) in enumerate(m.bodies)
    J[m.dim*(i-1)+1:m.dim*i,:] = translationalVelocityJacobian(body,m)
  end
  return J
end

function translationalVelocityJacobian(b::Body,m::DynaMo)
  return linEq(b.velocity,b.connection.u)
end

function J = linEq(M::Vector{Sym},x::Vector{Sym})
  return [coeff(i,j) for i in M, j in x]
end

function DOF(m::DynaMo)
  val = 0
  for body in m.bodies
    val += DOF(body)
  end
  return val
end

function DOF(b::Body)
  return length(b.q)
end


text:
"""
Continuation method implementations extend  `ContinuationMethod`
and the corresponding `show` and `step`.
Responsible for changing the project itself
"""
abstract ContinuationMethod

# generic show function. usually displays a GUI
Base.show(::ContinuationMethod) = error("implement!")

# step function. advances activeSolution by one step.
# atm called only by the cont gui.
# exposed so it can be called eg from a bifurcation view.
Base.step(::ContinuationMethod) = error("implement!")


text:
using Pkg
Pkg.update()
Pkg.develop(path = ".")
Pkg.add("DashBase")
Pkg.add("DashHtmlComponents")
Pkg.add("DashCoreComponents")
Pkg.add("DashTable")
Pkg.add("HTTP")
Pkg.build("Dash")
Pkg.build("DashHtmlComponents")
Pkg.build("DashCoreComponents")
Pkg.build("HTTP")
Pkg.test("Dash", coverage=true)
function precompile_pkgs()
     for pkg in collect(keys(Pkg.installed()))
        if !isdefined(Symbol(pkg), :Symbol) && pkg != "Compat.jl"
            @info("Importing $(pkg)...")
            try (@eval import $(Symbol(pkg))) catch end
        end
    end
end
precompile_pkgs()


text:
"""
    @! expr

Convert all supported mutating calls to double bang equivalent.

# Examples
```jldoctest
julia> using BangBang

julia> @! push!(empty!((0, 1)), 2, 3)
(2, 3)

julia> y = [1, 2];

julia> @! y .= 2 .* y
       y
2-element Array{Int64,1}:
 2
 4

julia> y = (1, 2);

julia> @! y .= 2 .* y
       y
(2, 4)
```
"""
macro !(expr)
    foldexpr(macroexpand(__module__, expr)) do x
        if Meta.isexpr(x, :call)
            isdotoperator(x.args[1]) && return x
            return Expr(:call, Expr(:call, _maybb, x.args[1]), x.args[2:end]...)
        elseif Meta.isexpr(x, :.=) && x.args[1] isa Symbol
            @assert length(x.args) == 2
            lhs, rhs = x.args
            return :($lhs = $materialize!!(
                $Base.@isdefined($lhs) ? $lhs : $(Undefined()),
                $air.($rhs),
            ))
        end
        return x
    end |> esc
end

foldexpr(f, x) = x
foldexpr(f, ex::Expr) = f(Expr(ex.head, foldexpr.(f, ex.args)...))

isdotoperator(x::Symbol) = undotoperator(x) !== nothing

function undotoperator(x::Symbol)
    startswith(string(x), ".") || return nothing
    op = Symbol(string(x)[2:end])
    Base.isoperator(op) || return nothing
    return op
end

function air end
struct Aired{T}
    value::T
end
@inline Broadcast.broadcasted(::typeof(air), x) = Aired(x)
@inline Broadcast.materialize(x::Aired) = x.value