text:
using BinaryBuilder

name = "MUMPS_seq_MKL"
version = v"5.3.5"

sources = [
  ArchiveSource("http://mumps.enseeiht.fr/MUMPS_$version.tar.gz",
                "e5d665fdb7043043f0799ae3dbe3b37e5b200d1ab7a6f7b2a4e463fd89507fa4"),
  DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
mkdir -p ${libdir}
cd $WORKSPACE/srcdir/MUMPS*
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/mumps_int32.patch

makefile="Makefile.INTEL.SEQ"
cp Make.inc/${makefile} Makefile.inc

libmkl=(-lmkl_intel_lp64 -lmkl_core -lmkl_intel_thread -liomp5 -lpthread -lm -ldl)
if [[ ${target} == *mingw* ]]; then
    libmkl=(${libdir}/mkl_rt.dll ${libdir}/mkl_intel_thread.dll ${libdir}/libiomp5md.dll ${libdir}/libwinpthread-1.dll)
fi

optf=(-O -DGEMMT_AVAILABLE)
optl=(-O)
optc=(-O)

make_args+=(OPTF=-O3
            CDEFS=-DAdd_
            LMETISDIR=${libdir}
            IMETIS=-I${prefix}/include
            LMETIS='-L$(LMETISDIR) -lmetis'
            ORDERINGSF="-Dpord -Dmetis"
            CC="$CC -fPIC ${CFLAGS[@]}"
            FC="gfortran -fPIC ${FFLAGS[@]}"
            FL="gfortran -fPIC"
            OPTF="${optf[@]}"
            OPTL="${optl[@]}"
            OPTC="${optc[@]}"
            LIBBLAS="${libmkl[@]}")

if [[ "${target}" == *-apple* ]]; then
  make_args+=(RANLIB=echo)
fi

# NB: parallel build fails
make alllib "${make_args[@]}"

# build shared libs
all_load="--whole-archive"
noall_load="--no-whole-archive"
extra=""
if [[ "${target}" == *-apple-* ]]; then
    all_load="-all_load"
    noall_load="-noall_load"
    extra="-Wl,-undefined -Wl,dynamic_lookup -headerpad_max_install_names"
fi

cd libseq
gfortran -fPIC -shared -Wl,${all_load} libmpiseq.a ${libs[@]} -Wl,${noall_load} ${extra[@]} -o libmpiseq.${dlext}
cp libmpiseq.${dlext} ${libdir}

cd ../lib
libs=(-L${libdir} -lmetis -lmpiseq)
gfortran -fPIC -shared -Wl,${all_load} libpord.a ${libs[@]} ${libmkl[@]} -Wl,${noall_load} ${extra[@]} -o libpord.${dlext}
cp libpord.${dlext} ${libdir}

libs+=(-lpord)
gfortran -fPIC -shared -Wl,${all_load} libmumps_common.a ${libs[@]} ${libmkl[@]} -Wl,${noall_load} ${extra[@]} -o libmumps_common.${dlext}
cp libmumps_common.${dlext} ${libdir}

libs+=(-lmumps_common)
for libname in cmumps dmumps smumps zmumps
do
  gfortran -fPIC -shared -Wl,${all_load} lib${libname}.a ${libs[@]} ${libmkl[@]} -Wl,${noall_load} ${extra[@]} -o lib${libname}.${dlext}
done
cp *.${dlext} ${libdir}
cd ..

mkdir -p ${prefix}/include/mumps_seq
cp include/* ${prefix}/include/mumps_seq
cp libseq/*.h ${prefix}/include/mumps_seq
"""

platforms = [
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("x86_64", "macos"),
    Platform("i686", "windows"),
    Platform("x86_64", "windows"),
]
platforms = expand_gfortran_versions(expand_cxxstring_abis(platforms))

# The products that we will ensure are always built
products = [
    LibraryProduct("libsmumps", :libsmumps),
    LibraryProduct("libdmumps", :libdmumps),
    LibraryProduct("libcmumps", :libcmumps),
    LibraryProduct("libzmumps", :libzmumps),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("CompilerSupportLibraries_jll"),
    Dependency("METIS_jll"),
    Dependency("MKL_jll"),
]

# Build the tarballs
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text: my_f(x, y) = 2 * x + 3 * y

ext: jl


text:
function setup(backend::GPUBackend, layer::BernoulliReconLossLayer, inputs::Vector{Blob}, diffs::Vector{Blob})
    data_type = eltype(inputs[1])
    aux_ones = make_blob(backend, data_type, size(inputs[1])...) # for summing
    fill!(aux_ones, 1.0)
    tmp1_blob = make_blob(backend, data_type, size(inputs[1])...)
    tmp2_blob = make_blob(backend, data_type, size(inputs[1])...)
    tmp = @compat(Dict(:aux_ones => aux_ones,
                       :tmp1 => tmp1_blob,
                       :tmp2 => tmp2_blob))
    state = BernoulliReconLossLayerState(layer, zero(data_type), zero(data_type), 0, tmp)
    return state
end

function forward(backend::GPUBackend, state::BernoulliReconLossLayerState, inputs::Vector{Blob})
    data_type = eltype(inputs[1])
    N = get_num(inputs[1])
    ND = length(inputs[1])

    p = inputs[1]
    x = inputs[2]

    aux_ones = state.tmp_blobs[:aux_ones]
    tmp1 = state.tmp_blobs[:tmp1]
    tmp2 = state.tmp_blobs[:tmp2]

    copy!(tmp1, aux_ones)
    CuVec.sub!(backend, tmp1, p)
    CuVec.log!(backend, tmp1)
    copy!(tmp2, aux_ones)
    CuVec.sub!(backend, tmp2, x)
    CuVec.mul!(backend, tmp1, tmp2)

    copy!(tmp2, p)
    CuVec.log!(backend, tmp2)
    CuVec.mul!(backend, tmp2, x)

    CuVec.add!(backend, tmp1, tmp2)
    ll = CuBLAS.dot(backend.cublas_ctx, data_type, ND, tmp1.ptr, 1, aux_ones.ptr, 1)

    state.loss = -ll * state.layer.weight / N

    # accumulate statistics
    state.loss_accum *= state.n_accum
    state.loss_accum += state.loss * N
    state.loss_accum /= state.n_accum + N

    state.n_accum += N
end

function backward(backend::GPUBackend, state::BernoulliReconLossLayerState,
            inputs::Vector{Blob}, diffs::Vector{Blob})
    data_type = eltype(inputs[1])
    N = get_num(inputs[1])
    ND = length(inputs[1])

    p = inputs[1]
    x = inputs[2]

    aux_ones = state.tmp_blobs[:aux_ones]
    tmp1 = state.tmp_blobs[:tmp1]
    tmp2 = state.tmp_blobs[:tmp2]

    c = convert(data_type, state.layer.weight / N)

    if isa(diffs[1], CuTensorBlob)
        copy!(diffs[1], x)
        CuVec.sub!(backend, diffs[1], aux_ones)
        copy!(tmp1, p)
        CuVec.sub!(backend, tmp1, aux_ones)
        CuVec.div!(backend, diffs[1], tmp1)

        copy!(tmp1, x)
        CuVec.div!(backend, tmp1, p)

        CuVec.sub!(backend, diffs[1], tmp1)
        CuVec.mul_scal!(backend, diffs[1], c)
    end

    if isa(diffs[2], CuTensorBlob)
        raise("last blob should be a data blob")
    end
end


text:
using Base.Test
using QuantumOptics

T = [0.:0.1:0.5;]

# Define Spin 1/2 operators
spinbasis = SpinBasis(1//2)
I = dense_identity(spinbasis)
sigmax = full(spin.sigmax(spinbasis))
sigmay = full(spin.sigmay(spinbasis))
sigmaz = full(spin.sigmaz(spinbasis))
sigmap = full(spin.sigmap(spinbasis))
sigmam = full(spin.sigmam(spinbasis))
I_spin = identity(spinbasis)

N = 2
b = tensor([spinbasis for i=1:N]...)

# psi0 = 1./sqrt(2)*(spinup(spinbasis) + spindown(spinbasis))
psi0 = normalize(spinup(spinbasis) + 0.5*spindown(spinbasis))
rho0 = cumulantexpansion.ProductDensityOperator([psi0 for i=1:N]...)

Ω = [1. 2. 3.;
     2. 1. 4.;
     3. 4. 1.]
# Ω = zeros(Float64, N, N)
γ = 1.
δ = 0.2

H = LazyTensor[LazyTensor(b, i, sigmaz, 0.5*δ) for i=1:N]
for i=1:N, j=1:N
    if i==j
        continue
    end
    h = LazyTensor(b, [i, j], [sigmap, sigmam], Ω[i, j])
    push!(H, h)
end
H = LazySum(H...)
J = LazyTensor[LazyTensor(b, i, sigmam, γ) for i=1:N]
# J = LazyTensor[]

tout, rho_t = cumulantexpansion.master(T, rho0, H, J)
tout, rho_t_full = timeevolution.master(T, full(rho0), full(H), map(full, J))

type ProductState
    N::Int
    data::Vector{Float64}
end

ProductState(N::Int) = ProductState(N, zeros(Float64, 3*N))

function ProductState(rho::DenseOperator)
    state = ProductState(N)
    sx, sy, sz = splitstate(state)
    f(ind, op) = real(expect(embed(b, ind, op), rho))
    for k=1:N
        sx[k] = f(k, sigmax)
        sy[k] = f(k, sigmay)
        sz[k] = f(k, sigmaz)
    end
    combinestate(sx, sy, sz, state.data)
    return state
end

splitstate(N::Int, data::Vector{Float64}) = vec(data[0*N+1:1*N]), vec(data[1*N+1:2*N]), vec(data[2*N+1:3*N])
splitstate(state::ProductState) = splitstate(state.N, state.data)
function combinestate(sx::Vector{Float64}, sy::Vector{Float64}, sz::Vector{Float64}, state::Vector{Float64})
    state[0*N+1:1*N] = sx
    state[1*N+1:2*N] = sy
    state[2*N+1:3*N] = sz
    state
end

function f(t, y::Vector{Float64}, dy::Vector{Float64})
    sx, sy, sz = splitstate(N, y)
    dsx, dsy, dsz = splitstate(N, dy)
    @inbounds for k=1:N
        dsx[k] = -δ*sy[k] - 0.5*γ*sx[k]
        dsy[k] = δ*sx[k] - 0.5*γ*sy[k]
        dsz[k] = -γ*(1+sz[k])
        for j=1:N
            if j==k
                continue
            end
            dsx[k] += Ω[k,j]*sy[j]*sz[k]
            dsy[k] += -Ω[k,j]*sx[j]*sz[k]
            dsz[k] += Ω[k,j]*(sx[j]*sy[k] - sy[j]*sx[k])
        end
    end
    combinestate(dsx, dsy, dsz, dy)
end

t_out = Float64[]
state_out = ProductState[]

function fout_(t, y::Vector{Float64})
    push!(t_out, t)
    push!(state_out, ProductState(N, deepcopy(y)))
end

function densityoperator(sx::Real, sy::Real, sz::Real)
    return 0.5*(I + sx*sigmax + sy*sigmay + sz*sigmaz)
end
function densityoperator(state::ProductState)
    sx, sy, sz = splitstate(state)
    rho = densityoperator(sx[1], sy[1], sz[1])
    for i=2:state.N
        rho = tensor(rho, densityoperator(sx[i], sy[i], sz[i]))
    end
    return rho
end

state0 = ProductState(full(rho0))
QuantumOptics.ode_dopri.ode(f, T, state0.data, fout_)

@test tracedistance(densityoperator(state_out[end]), full(rho_t[end])) < 1e-6
# println(tracedistance(densityoperator(state_out[end]), full(rho_t_full[end])))
@test tracedistance(full(rho_t[end]), full(rho_t[1])) > 0.1


text:
for N in [Float64, Rational{Int}, Float32]
    X = Hyperrectangle(zeros(N, 2), ones(N, 2))
    U = underapproximate(X, BoxDirections{N}(2))
    @test U ⊆ X
end


text:
using UnPack, OrderedCollections, LinearAlgebra, JLD2, LabelledArrays, SparseArrays

# Representative agent version, which is essentially Bansal and Yaron (2004)
mutable struct BansalYaron2004{T <: Real}
    p::LArray{T, 1, Array{T, 1}, (:μ_y, :ρ_x, :σ_x, :ρ_σ, :σ_y, :ς, :β, :ψ, :γ)}
    N_approx::LArray{Int64, 1, Array{Int64, 1}, (:q, :ω)}
    S::OrderedDict{Symbol, Int}
    J::OrderedDict{Symbol, Int}
    E::OrderedDict{Symbol, Int}
    SH::OrderedDict{Symbol, Int}
end

# Parameters are based off Schorfheide et al. (2016) and converted from a monthly frequency to quarterly
function BansalYaron2004(; μ_y::T = 0.0016 * 3., ρ_x::T = 0.99^3, σ_x::T = sqrt((0.74 * sqrt(1. - ρ_x^2))^2 * 3.), ρ_σ::T = 0.99^3,
                         σ_y::T = sqrt(0.0021^2 * 3.), ς::T = sqrt(0.0014^2 * 3.),
                         β::T = 0.999^3, ψ::T = 2., γ::T = 9.,
                         N_approx::LArray{Int64, 1, Array{Int64, 1}, (:q, :ω)} = (@LArray [1, 1] (:q, :ω))) where {T <: Real}

    @assert all(N_approx[k] > 0 for k in keys(N_approx)) "N_approx must be at least 1 for all variables."

    # Create indexing dictionaries
    S_init  = [:yy, :x, :σ²_y]                                        # State variables
    J_init  = [:q, :v, :ce, :ω]                                       # Jump variables
    SH_init = [:ε_y, :ε_x, :ε_σ²]                                     # Exogenous shocks
    E_init  = [:value_fnct, :certainty_equiv, :ez_fwd_diff, :cap_ret] # Equations
    for var in [:q, :ω]
        inds = (var == :q) ? (1:N_approx[var]) : (0:(N_approx[var] - 1))
        push!(J_init, [Symbol(:d, var, "$(i)") for i in inds]...)
        push!(J_init, [Symbol(:p, var, "$(i)") for i in 1:N_approx[var]]...)
        push!(E_init, [Symbol(:eq_d, var, "$(i)") for i in inds]...)
        push!(E_init, [Symbol(:eq_p, var, "$(i)") for i in 1:N_approx[var]]...)
    end

    S  = OrderedDict{Symbol, Int}(k => i for (i, k) in enumerate(S_init))
    J  = OrderedDict{Symbol, Int}(k => i for (i, k) in enumerate(J_init))
    E  = OrderedDict{Symbol, Int}(k => i for (i, k) in enumerate(E_init))
    SH = OrderedDict{Symbol, Int}(k => i for (i, k) in enumerate(SH_init))

    para = @LArray [μ_y, ρ_x, σ_x, ρ_σ, σ_y, ς, β, ψ, γ] (:μ_y, :ρ_x, :σ_x, :ρ_σ, :σ_y, :ς, :β, :ψ, :γ)
    return BansalYaron2004{T}(para, N_approx, S, J, E, SH)
end

function bansal_yaron_2004(m::BansalYaron2004{T}; sparse_arrays::Bool = false,
                           sparse_jacobian::Vector{Symbol} = Symbol[]) where {T <: Real}

    # Unpack parameters and indexing dictionaries
    @unpack p, N_approx, S, J, E, SH = m

    @unpack yy, x, σ²_y = S
    @unpack q, v, ce, ω = J
    @unpack ε_y, ε_x, ε_σ² = SH
    @unpack value_fnct, certainty_equiv, ez_fwd_diff, cap_ret = E

    Nz = length(S)
    Ny = length(J)
    Nε = length(SH)

    ## Define nonlinear equations

    # Some helper functions
    m_ξ(z, y) = log(p.β) - (p.ψ - p.γ) * y[ce] - p.γ * p.μ_y
    function m_fwd!(i, Γ₅, Γ₆)
        Γ₅[i, yy] = -p.γ
        Γ₆[i, v]  = (p.ψ  - p.γ)
    end

    function μ(F, z, y) # note that y here refers to jump variables
        F[yy]   = z[x]
        F[x]    = p.ρ_x * z[x]
        F[σ²_y] = (1. - p.ρ_σ) * p.σ_y^2 + p.ρ_σ * z[σ²_y]
    end

    function ξ(F, z, y)

        m_ξv = m_ξ(z, y) # evaluate SDF
        F[value_fnct] = 1. / (1. - p.ψ) * (log(1. - p.β) + y[ω]) - y[v]
        F[certainty_equiv] = 1. / (1. - p.ψ) * (log(1. - p.β) - log(p.β) + log(exp(y[ω]) - 1.)) - y[ce]

        ## Forward-difference equations separately handled b/c recursions
        F[cap_ret]     = y[q] - log(sum([exp(y[J[Symbol("dq$(i)")]]) for i in 1:N_approx[:q]]) +
                                     exp(y[J[Symbol("pq$(N_approx[:q])")]]))
        F[ez_fwd_diff] = y[ω] - log(sum([exp(y[J[Symbol("dω$(i)")]]) for i in 0:(N_approx[:ω] - 1)]) +
                                     exp(y[J[Symbol("pω$(N_approx[:ω])")]]))

        # Set initial boundary conditions
        F[E[:eq_dq1]] = p.μ_y - y[J[:dq1]] + m_ξv
        F[E[:eq_pq1]] = p.μ_y - y[J[:pq1]] + m_ξv
        F[E[:eq_dω0]] = y[J[:dω0]]
        F[E[:eq_pω1]] = p.μ_y - y[J[:pω1]] + m_ξv

        # Recursions for forward-difference equations
        for i in 2:N_approx[:q]
            F[E[Symbol("eq_dq$(i)")]] = p.μ_y - y[J[Symbol("dq$(i)")]] + m_ξv
            F[E[Symbol("eq_pq$(i)")]] = p.μ_y - y[J[Symbol("pq$(i)")]] + m_ξv
        end
        for i in 2:N_approx[:ω]
            F[E[Symbol("eq_dω$(i-1)")]] = p.μ_y - y[J[Symbol("dω$(i-1)")]] + m_ξv
            F[E[Symbol("eq_pω$(i)")]]   = p.μ_y - y[J[Symbol("pω$(i)")]]   + m_ξv
        end
    end

    # The cache is initialized as zeros so we only need to fill non-zero elements
    Λ = sparse_arrays ? spzeros(T, Nz, Ny) : zeros(T, Nz, Ny)

    # The cache is initialized as zeros so we only need to fill non-zero elements
    function Σ(F, z)
        F[yy,   ε_y]  = sqrt(z[σ²_y])
        F[x,    ε_x]  = sqrt(z[σ²_y]) * p.σ_x
        F[σ²_y, ε_σ²] = sqrt(z[σ²_y]) * p.ς
    end

    function ccgf(F, α, z)
        # F .= .5 * RiskAdjustedLinearizations.diag(α * α') # slower but this is the underlying math
        sum!(F, α.^2) # faster implementation
        F .*= .5
    end

    if sparse_arrays
        Γ₅ = spzeros(T, Ny, Nz)
        Γ₆ = spzeros(T, Ny, Ny)
    else
        Γ₅ = zeros(T, Ny, Nz)
        Γ₆ = zeros(T, Ny, Ny)
    end

    # Forward difference equations: boundary conditions
    m_fwd!(E[:eq_dq1], Γ₅, Γ₆)
    Γ₅[E[:eq_dq1], yy] = one(T)

    m_fwd!(E[:eq_pq1], Γ₅, Γ₆)
    Γ₅[E[:eq_pq1], yy] = one(T)
    Γ₆[E[:eq_pq1], q]  = one(T)

    m_fwd!(E[:eq_pω1], Γ₅, Γ₆)
    Γ₅[E[:eq_pω1], yy] = one(T)
    Γ₆[E[:eq_pω1], ω]  = one(T)

    # Forward difference equations: recursions
    for i in 2:N_approx[:q]
        m_fwd!(E[Symbol("eq_dq$(i)")], Γ₅, Γ₆)
        Γ₅[E[Symbol("eq_dq$(i)")], yy] = one(T)
        Γ₆[E[Symbol("eq_dq$(i)")], J[Symbol("dq$(i-1)")]] = one(T)

        m_fwd!(E[Symbol("eq_pq$(i)")], Γ₅, Γ₆)
        Γ₅[E[Symbol("eq_pq$(i)")], yy] = one(T)
        Γ₆[E[Symbol("eq_pq$(i)")], J[Symbol("pq$(i-1)")]] = one(T)
    end
    for i in 2:N_approx[:ω]
        m_fwd!(E[Symbol("eq_dω$(i-1)")], Γ₅, Γ₆)
        Γ₅[E[Symbol("eq_dω$(i-1)")], yy] = one(T)
        Γ₆[E[Symbol("eq_dω$(i-1)")], J[Symbol("dω$(i-2)")]] = one(T)

        m_fwd!(E[Symbol("eq_pω$(i)")], Γ₅, Γ₆)
        Γ₅[E[Symbol("eq_pω$(i)")], yy] = one(T)
        Γ₆[E[Symbol("eq_pω$(i)")], J[Symbol("pω$(i-1)")]] = one(T)
    end

    z, y = create_deterministic_guess(m)
    Ψ = zeros(T, Ny, Nz)

    if sparse_arrays
        return RiskAdjustedLinearization(μ, Λ, Σ, ξ, Γ₅, Γ₆, ccgf, vec(z), vec(y), Ψ, Nε; sparse_jacobian = sparse_jacobian,
                                         Λ_cache_init = dims -> spzeros(dims...), Σ_cache_init = dims -> spzeros(dims...))
    else
        return RiskAdjustedLinearization(μ, Λ, Σ, ξ, Γ₅, Γ₆, ccgf, vec(z), vec(y), Ψ, Nε; sparse_jacobian = sparse_jacobian)
    end
end

function create_deterministic_guess(m::BansalYaron2004{T}) where {T <: Real}

    ## Set up

    # Unpack parameters and indexing dictionaries
    @unpack p, N_approx, S, J, E, SH = m

    @unpack yy, x, σ²_y = S
    @unpack q, v, ce, ω = J
    @unpack ε_y, ε_x, ε_σ² = SH
    @unpack value_fnct, certainty_equiv, ez_fwd_diff, cap_ret = E
    Nz = length(S)
    Ny = length(J)

    # Initialize deterministic steady state guess vectors
    z = Vector{T}(undef, Nz)
    y = Vector{T}(undef, Ny)

    ## Compute guesses

    # Steady state values of state variables known ex-ante
    z[yy]   = 0.
    z[x]    = 0.
    z[σ²_y] = p.σ_y^2

    # Now make guesses for remaining quantities
    Y0  = 1. # long-run value is long-run value of X, which is 1
    Ω0  = 1. / (1. - (p.β * Y0 * exp(p.μ_y)) ^ (1. - p.ψ))
    V0  = ((1. - p.β) * Ω0) ^ (1. / (1. - p.ψ))
    CE0 = ((1. - p.β) / p.β * (Ω0 - 1.)) ^ (1. / (1. - p.ψ))
    M0  = p.β * (p.β * Ω0 / (Ω0 - 1.)) ^ ((p.ψ - p.γ) / (1. - p.ψ)) * (Y0 * exp(p.μ_y)) ^ (-p.γ)
    Q0  = exp(p.μ_y) * M0 * Y0 / (1. - exp(p.μ_y) * M0 * Y0)

    y[q]  = log(Q0)
    y[v]  = log(V0)
    y[ce] = log(CE0)
    y[ω]  = log(Ω0)

    y[J[:dq1]] = convert(T, log(exp(p.μ_y) * M0 * Y0))
    y[J[:pq1]] = convert(T, log(exp(p.μ_y) * M0 * Y0 * Q0))
    y[J[:dω0]] = zero(T)
    y[J[:pω1]] = convert(T, log(exp(p.μ_y) * M0 * Y0 * Ω0))

    for i in 2:N_approx[:q]
        y[J[Symbol("dq$(i)")]] = convert(T, log(M0) + μ_y + log(Y0) + y[J[Symbol("dq$(i-1)")]])
        y[J[Symbol("pq$(i)")]] = convert(T, log(M0) + μ_y + log(Y0) + y[J[Symbol("pq$(i-1)")]])
    end
    for i in 2:N_approx[:ω]
        y[J[Symbol("dω$(i-1)")]] = convert(T, μ_y + log(M0) + log(Y0) + y[J[Symbol("dω$(i-2)")]])
        y[J[Symbol("pω$(i)")]] = convert(T, μ_y + log(M0) + log(Y0) + y[J[Symbol("pω$(i-1)")]])
    end

    return z, y
end


text:
module TestConversions
    using Compat, Compat.Test, DataFrames, Compat.InteractiveUtils
    using DataStructures: OrderedDict, SortedDict
    const ≅ = isequal

    df = DataFrame()
    df[:A] = 1:5
    df[:B] = [:A, :B, :C, :D, :E]
    @test isa(convert(Array, df), Matrix{Any})
    @test isa(convert(Array{Any}, df), Matrix{Any})

    df = DataFrame()
    df[:A] = 1:5
    df[:B] = 1.0:5.0
    @test isa(convert(Array, df), Matrix{Float64})
    @test isa(convert(Array{Any}, df), Matrix{Any})
    @test isa(convert(Array{Float64}, df), Matrix{Float64})

    df = DataFrame()
    df[:A] = Vector{Union{Float64, Missing}}(1.0:5.0)
    df[:B] = Vector{Union{Float64, Missing}}(1.0:5.0)
    a = convert(Array, df)
    aa = convert(Array{Any}, df)
    ai = convert(Array{Int}, df)
    @test isa(a, Matrix{Union{Float64, Missing}})
    @test a == convert(Array, convert(Array{Union{Float64, Missing}}, df))
    @test a == convert(Matrix, df)
    @test isa(aa, Matrix{Any})
    @test aa == convert(Matrix{Any}, df)
    @test isa(ai, Matrix{Int})
    @test ai == convert(Matrix{Int}, df)

    df[1,1] = missing
    @test_throws ErrorException convert(Array{Float64}, df)
    na = convert(Array{Union{Float64, Missing}}, df)
    naa = convert(Array{Union{Any, Missing}}, df)
    nai = convert(Array{Union{Int, Missing}}, df)
    @test isa(na, Matrix{Union{Float64, Missing}})
    @test na ≅ convert(Matrix, df)
    @test isa(naa, Matrix{Union{Any, Missing}})
    @test naa ≅ convert(Matrix{Union{Any, Missing}}, df)
    @test isa(nai, Matrix{Union{Int, Missing}})
    @test nai ≅ convert(Matrix{Union{Int, Missing}}, df)

    a = Union{Float64, Missing}[1.0,2.0]
    b = Union{Float64, Missing}[-0.1,3]
    c = Union{Float64, Missing}[-3.1,7]
    di = Dict("a"=>a, "b"=>b, "c"=>c)

    df = convert(DataFrame, di)
    @test isa(df, DataFrame)
    @test names(df) == Symbol[x for x in sort(collect(keys(di)))]
    @test df[:a] == a
    @test df[:b] == b
    @test df[:c] == c

    od = OrderedDict("c"=>c, "a"=>a, "b"=>b)
    df = convert(DataFrame,od)
    @test isa(df, DataFrame)
    @test names(df) == Symbol[x for x in keys(od)]
    @test df[:a] == a
    @test df[:b] == b
    @test df[:c] == c

    sd = SortedDict("c"=>c, "a"=>a, "b"=>b)
    df = convert(DataFrame,sd)
    @test isa(df, DataFrame)
    @test names(df) == Symbol[x for x in keys(sd)]
    @test df[:a] == a
    @test df[:b] == b
    @test df[:c] == c

    a = [1.0]
    di = Dict("a"=>a, "b"=>b, "c"=>c)
    @test_throws DimensionMismatch convert(DataFrame,di)

end


text:
module PlutoUI

import Base: show, get
import Markdown: htmlesc, withtag

using Reexport

const PKG_ROOT_DIR = normpath(joinpath(@__DIR__, ".."))

include("./Builtins.jl")
include("./Resource.jl")
include("./Terminal.jl")
include("./RangeSlider.jl")
include("./DisplayTricks.jl")

@reexport module MultiCheckBoxNotebook
    include("./MultiCheckBox.jl")
end
@reexport module TableOfContentsNotebook
    include("./TableOfContents.jl")
end
@reexport module ClockNotebook
    include("./Clock.jl")
end
@reexport module ScrubbableNotebook
    include("./Scrubbable.jl")
end

end


text:
# Implemment the mutable arithmetics interface for UnitAffExpr

_MA.mutability(::Type{UnitAffExpr}) = _MA.IsMutable()

#  Minimal support of operate!! to allow @rewrite macro to work on linear constraints

UnitAffOrVar = Union{UnitVariableRef, UnitAffExpr}
AddSub = Union{typeof(_MA.add_mul), typeof(_MA.sub_mul)}
NumQuant = Union{Number, Quantity}

function Base.convert(::Type{UnitAffExpr}, uv::UnitVariableRef)
    return UnitAffExpr(convert(AffExpr, uv.vref), uv.u)
end

function update_expression!(ua::UnitAffExpr, a::NumQuant, x::UnitVariableRef)
    aval = ustrip(uconvert(ua.u, a * Quantity(1, x.u)))
    return UnitAffExpr(add_to_expression!(ua.expr, aval, x.vref), ua.u)
end

function update_expression!(ua::UnitAffExpr, a::Quantity, x::VariableRef)
    aval = ustrip(uconvert(ua.u, a))
    return UnitAffExpr(add_to_expression!(ua.expr, aval, x), ua.u)
end

function update_expression!(ua::UnitAffExpr, a::Quantity)
    aval = ustrip(uconvert(ua.u, a))
    return UnitAffExpr(add_to_expression!(ua.expr, aval), ua.u)
end

function update_expression!(ua::UnitAffExpr, x::UnitAffExpr)
    factor = ustrip(uconvert(ua.u, Quantity(1, x.u)))
    return UnitAffExpr(add_to_expression!(ua.expr, factor, x.expr), ua.u)
end

function update_expression!(ua::UnitAffExpr, a::NumQuant, b::NumQuant)
    aval = ustrip(uconvert(ua.u, a * b))
    return UnitAffExpr(add_to_expression!(ua.expr, aval), ua.u)
end

function create_expression(t::AddSub, z::typeof(_MA.Zero()), a::NumQuant, x::UnitVariableRef)
    val = a * Quantity(1, x.u)
    return UnitAffExpr(_MA.operate!!(t, z, ustrip(val), x.vref), Unitful.unit(val))
end

function create_expression(t::AddSub, z::typeof(_MA.Zero()), a::Quantity, x::VariableRef)
    return UnitAffExpr(_MA.operate!!(t, z, ustrip(a), x), Unitful.unit(a))
end

function create_expression(t::AddSub, z::typeof(_MA.Zero()), a::Quantity)
    return UnitAffExpr(_MA.operate!!(t, z, ustrip(a)), Unitful.unit(a))
end

# Two arguments
_MA.operate!!(::typeof(_MA.add_mul), uav::UnitAffOrVar, x::UnitVariableRef) = update_expression!(convert(UnitAffExpr, uav), 1, x)
_MA.operate!!(::typeof(_MA.sub_mul), uav::UnitAffOrVar, x::UnitVariableRef) = update_expression!(convert(UnitAffExpr, uav), -1, x)

_MA.operate!!(::typeof(_MA.add_mul), uav::UnitAffOrVar, x::UnitAffExpr) = update_expression!(convert(UnitAffExpr, uav), x)
_MA.operate!!(::typeof(_MA.sub_mul), uav::UnitAffOrVar, x::UnitAffExpr) = update_expression!(convert(UnitAffExpr, uav), -1, x)

_MA.operate!!(::typeof(_MA.add_mul), uav::UnitAffOrVar, a::Quantity) = update_expression!(convert(UnitAffExpr, uav), a)
_MA.operate!!(::typeof(_MA.sub_mul), uav::UnitAffOrVar, a::Quantity) = update_expression!(convert(UnitAffExpr, uav), -a)

# Three arguments
_MA.operate!!(::typeof(_MA.add_mul), uav::UnitAffOrVar, a::NumQuant, x::UnitVariableRef) = update_expression!(convert(UnitAffExpr, uav), a, x)
_MA.operate!!(::typeof(_MA.sub_mul), uav::UnitAffOrVar, a::NumQuant, x::UnitVariableRef) = update_expression!(convert(UnitAffExpr, uav), -a, x)

_MA.operate!!(::typeof(_MA.add_mul), uav::UnitAffOrVar, a::Quantity, x::VariableRef) = update_expression!(convert(UnitAffExpr, uav), a, x)
_MA.operate!!(::typeof(_MA.sub_mul), uav::UnitAffOrVar, a::Quantity, x::VariableRef) = update_expression!(convert(UnitAffExpr, uav), -a, x)

_MA.operate!!(::typeof(_MA.add_mul), uav::UnitAffOrVar, a::NumQuant, b::NumQuant) = update_expression!(convert(UnitAffExpr, uav), a, b)
_MA.operate!!(::typeof(_MA.sub_mul), uav::UnitAffOrVar, a::NumQuant, b::NumQuant) = update_expression!(convert(UnitAffExpr, uav), -a, b)

_MA.operate!!(::typeof(_MA.add_mul), uav::UnitAffOrVar, a::Number, b::Unitful.Units) = update_expression!(convert(UnitAffExpr, uav), a * Quantity(1.0, b))
_MA.operate!!(::typeof(_MA.sub_mul), uav::UnitAffOrVar, a::Number, b::Unitful.Units) = update_expression!(convert(UnitAffExpr, uav), - a * Quantity(1.0, b))

# Multiple arguments
function _MA.operate!!(t::AddSub, uav::UnitAffOrVar, x, y, z, other_args...)
    args = [x, y, z, other_args...]
    n = length(args)
    varidx = findall(t -> (typeof(t) <: UnitAffExpr) || (typeof(t) <: UnitVariableRef) || (typeof(t) <: VariableRef), args)
    var = args[varidx[1]]
    val = prod(args[i] for i in setdiff(1:n, varidx))
    return _MA.operate!!(t, uav, val, var)
end

# Zero handling
_MA.operate!!(t::AddSub, z::_MA.Zero, a::Quantity) = create_expression(t, z, a)
_MA.operate!!(t::AddSub, z::_MA.Zero, x::UnitVariableRef) = create_expression(t, z, 1, x)
_MA.operate!!(t::AddSub, z::_MA.Zero, ua::UnitAffExpr) = ua
_MA.operate!!(t::AddSub, z::_MA.Zero, a::NumQuant, x::UnitVariableRef) = create_expression(t, z, a, x)
_MA.operate!!(t::AddSub, z::_MA.Zero, x::UnitVariableRef, a::NumQuant) = create_expression(t, z, a, x)
_MA.operate!!(t::AddSub, z::_MA.Zero, a::Quantity, x::VariableRef) = create_expression(t, z, a, x)
_MA.operate!!(t::AddSub, z::_MA.Zero, x::VariableRef, a::Quantity) = create_expression(t, z, a, x)

function _MA.operate!!(t::AddSub, zero::_MA.Zero, x, y, z, other_args...)
    args = [x, y, z, other_args...]
    n = length(args)
    varidx = findall(t -> (typeof(t) <: UnitAffExpr) || (typeof(t) <: UnitVariableRef) || (typeof(t) <: VariableRef), args)
    var = args[varidx[1]]
    val = prod(args[i] for i in setdiff(1:n, varidx))
    return _MA.operate!!(t, zero, val, var)
end


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
import Base.cat
using PyCall
p = PyCall.pyimport("pyplasm")
using DataStructures

"""
	apply(fun::Function, params::Array)

"""
function apply(fun::Function, params::Array)
	return fun(params)
end
function apply(fun::Function)
	function apply0(params::Array)
		return fun(params)
	end
	return apply0
end

"""
	apply(affineMatrix::Matrix)(larmodel::LAR)::LAR

Apply the `affineMatrix` parameter to the vertices of `larmodel`.

# Example

```
julia> square = LinearAlgebraicRepresentation.cuboid([1,1])
([0.0 0.0 1.0 1.0; 0.0 1.0 0.0 1.0], Array{Int64,1}[[1, 2, 3, 4]])

julia> Plasm.apply(LinearAlgebraicRepresentation.t(1,2))(square)
([1.0 1.0 2.0 2.0; 2.0 3.0 2.0 3.0], Array{Int64,1}[[1, 2, 3, 4]])
```
"""
function apply(affineMatrix)
	function apply0(larmodel)
		return Lar.struct2lar(Lar.Struct([ affineMatrix,larmodel ]))
	end
	return apply0
end

"""
	comp(funs::Array)

Standard mathematical composition.

Pipe execution from right to left on application to actual parameter.
"""
function comp(funs)
    function compose(f,g)
	  return x -> f(g(x))
	end
    id = x->x
    return reduce(compose, funs; init=id)
end

"""
	cons(funs::Array)(x::Any)::Array

*Construction* functional of FL and PLaSM languages.

Provides a *vector* functional that returns the array of
applications of component functions to actual parameter.

# Example

```
julia> Plasm.cons([cos,sin])(0)
2-element Array{Float64,1}:
 1.0
 0.0
```
"""
function cons(funs)
	return x -> [f(x) for f in funs]
end

"""
	k(Any)(x)

*Constant* functional of FL and PLaSM languages.

Gives a constant functional that always returns the actual parameter
when applied to another parameter.

#	Examples

```
julia> Plasm.k(10)(100)
10

julia> Plasm.k(sin)(cos)
sin
```
"""
function k(Any)
	x->Any
end

"""
	aa(fun::Function)(args::Array)::Array

AA applies fun to each element of the args sequence

# Example

```
julia> Plasm.aa(sqrt)([1,4,9,16])
4-element Array{Float64,1}:
 1.0
 2.0
 3.0
 4.0
```
"""
function aa(fun)
	function aa1(args::Array)
		map(fun,args)
	end
	return aa1
end

"""
	id(x::Anytype)

Identity function.  Return the argument.

"""
id = x->x

"""
	distr(args::Union{Tuple,Array})::Array

Distribute right. The parameter `args` must contain a `list` and an element `x`.
Return the `pair` array with the elements of `args` coupled with `x`

# Example

```
julia> Plasm.distr(([1,2,3],10))
3-element Array{Array{Int64,1},1}:
 [1, 10]
 [2, 10]
 [3, 10]
```
"""
function distr(args)
	list,element = args
	return [ [e,element] for e in list ]
end

"""
	distl(args::Union{Tuple,Array})::Array

Distribute right. The parameter `args` must contain an element `x` and a `list`.
Return the `pair` array with `x` coupled with the elements of `args`.

# Example

```
julia> Plasm.distl((10, [1,2,3]))
3-element Array{Array{Int64,1},1}:
 [10, 1]
 [10, 2]
 [10, 3]
```
"""
function distl(args)
	element, list = args
	return [ [element, e] for e in list ]
end


text:
module CommunalHelperDreamFloatySpaceBlock

using ..Ahorn, Maple
using Ahorn.CommunalHelper
using Ahorn.CommunalHelperEntityPresets: CustomDreamBlockData

@mapdefdata Entity "CommunalHelper/DreamFloatySpaceBlock" DreamFloatySpaceBlock CustomDreamBlockData

const placements = Ahorn.PlacementDict(
    "Dream Floaty Space Block (Communal Helper)" => Ahorn.EntityPlacement(
        DreamFloatySpaceBlock,
        "rectangle"
    ),
)

Ahorn.minimumSize(entity::DreamFloatySpaceBlock) = 8, 8
Ahorn.resizable(entity::DreamFloatySpaceBlock) = true, true

Ahorn.selection(entity::DreamFloatySpaceBlock) = Ahorn.getEntityRectangle(entity)

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::DreamFloatySpaceBlock, room::Maple.Room)
    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    renderDreamBlock(ctx, 0, 0, width, height, entity.data)
end

end


text:
## Basic SMC algorithm without rejuvenation

function basic_smc(T::Int, observations::ChoiceMap, n_particles::Int;
                   obs_noise::Float64=5.0)
    # Initialize filter
    state = pf_initialize(bayesian_sir, (0, obs_noise), choicemap(), n_particles)
    # Step through filter
    println("Running basic SMC...")
    @showprogress for t=1:T
        obs = get_selected(observations, Gen.select(:step => t))
        pf_update!(state, (t, obs_noise), (UnknownChange(), NoChange()), obs)
    end
    # Return particles and log ML estimate
    return state
end


text:
using Glutils
using Test

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

max_stars_repo_name: jorge-brito/Glutils

max_issues_repo_name: jorge-brito/Glutils

max_forks_repo_name: jorge-brito/Glutils


text:
@testset "Treebanks" begin

    datadir = joinpath(@__DIR__, "data")

    @testset "Treebanks" begin
        files = [joinpath(datadir, file) for file in readdir(datadir)
                 if endswith(file, ".conllu")]
        @test length(files) == 3

        @test_throws Exception Treebank("not a corpus")

        treebank = Treebank(files)
        @test showstr(treebank) == "Treebank (3 files)"
        trees = collect(treebank)
        @test length(trees) == 4

        np = Treebank(joinpath(datadir, "nonprojective.conll"))
        @test length(collect(np)) == 3
        @test length(filter(is_projective, collect(np))) == 2
    end

    @testset "Oracles & Projectivity" begin
        tb = Treebank(joinpath(datadir, "nonprojective1.conll"))
        np_oracle = Oracle(ListBasedNonProjective(), static_oracle, untyped)
        p_oracle = Oracle(ArcEager(), static_oracle, untyped)
        np_data = collect(Iterators.flatten(np_oracle.(tb)))
        p_data = collect(Iterators.flatten(p_oracle.(tb)))
        @test length(np_data) > 1
        @test length(p_data) == 0
        @test length(p_oracle(first(tb))) == 0

        count = 0
        for (cfg, t) in p_oracle(first(tb))
            count += 1
        end
        @test count == 0
    end
end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

"""Encapsulates the lens assembly and the display for one lenslet."""

struct Display{T} <: Surface{T}
    surface::Rectangle{T}
    pixellattice::LatticeBasis{2,T}
    pixelpitch::Tuple{T,T} #x,y pitch of pixels
    xresolution::T
    yresolution::T
    transform::Geometry.Transform{T}

    function Display(xres::Int64,yres::Int64,xpitch::Unitful.Length,ypitch::Unitful.Length,transform::Geometry.Transform{T}) where{T<:Real}
        # convert xpitch,ypitch to unitless numbers representing μm
        pitchx = ustrip(μm,xpitch)
        pitchy = ustrip(μm,ypitch)
        size = (xres*pitchx,yres*pitchy)
        surface = Rectangle(size[1]/T(2),size[2]/T(2),SVector(T(0),T(0),T(1)),SVector(T(0),T(0),T(0))) #surface normal is +z axis
        pixellattice = Repeat.rectangularlattice(size[2],size[1]) #lattice takes ypitch,xpitch in that order
        return new{T}(surface,pixellattice,(pitchx,pitchy),xres,yres,transform)
    end
end

surfaceintersection(a::Display) = surfaceintersection(a.surface)
normal(a::Display) = normal(a.surface)
interface(a::Display) = interface(a.surface)
makemesh(a::Display) = makemesh(a.surface)
uv(a::Display,p::SVector{3,T}) where{T<:Real} = uv(a.surface,p)
onsurface(a::Display, point::SVector{3,T}) where {T<:Real} = onsurface(a.surface,point)
pixelposition(xindex::Int,yindex::Int) = transform*pixellattice[yindex,xindex]
struct LensletAssembly{T}
    lens::ParaxialLens{T}
    transform::Geometry.Transform{T}
    display::Display{T}
    worldtodisplay::Geometry.Transform{T}

    LensletAssembly(lens::ParaxialLens{T},transform::Geometry.Transform{T},display::Display{T}) where{T} = new{T}(lens,transform,display,display.transform*transform)
end

lens(a::LensletAssembly) = a.lens
display(a::LensletAssembly) = a.display
transform(a::LensletAssembly) = a.transform

worldtolens(a::LensletAssembly, pt::SVector{3})= a.transform*pt
worldtolens(a::LensletAssembly,mat::SMatrix{3}) = a.transform*mat
lenstodisplay(a::LensletAssembly,pt::SVector{3,T}) where{T<:Real} = a.display.transform*pt
worldtodisplay(a::LensletAssembly,pt::SVector{3,T}) where{T<:Real} = a.worldtodisplay*pt

vertices3d(vertices::SMatrix{2,N,T}) where{N,T<:Real} = vcat(vertices,zero(SMatrix{1,N,T}))
export vertices3d

"""Projects vertexpoints in world space onto the lens plane and converts them to two dimensional points represented in the local x,y coordinates of the lens coordinate frame. Used for projecting eye pupil onto lens plane."""
function project(lenslet::LensletAssembly{T},displaypoint::SVector{3,T},vertexpoints::SMatrix{3,N,T}) where{T<:Real,N}
    projectedpoints = MMatrix{2,N,T}(undef)
    locvertices = worldtolens(lenslet,vertexpoints) #transform pupil vertices into local coordinate frame of lens
    virtpoint = point(virtualpoint(lens(lenslet), displaypoint)) #compute virtual point corresponding to physical display point
    for i in 1:N
        ppoint = locvertices[:,i]
        vec = ppoint-virtpoint
        vecdist = ppoint[3]
        virtdist = virtpoint[3]
        scale =  abs(virtdist/(vecdist-virtdist))
        planepoint = scale*vec + virtpoint
        projectedpoints[1:2,i] = SVector{2,T}(planepoint[1],planepoint[2]) #local lens coordinate frame has z axis aligned with the positive normal to the lens plane.
    end

    return SMatrix{2,N,T}(projectedpoints)
end

function convertlazysets(verts::LazySets.VectorIterator)
    M = length(verts)
    T = eltype(eltype(verts))
    temp = MMatrix{3,M,T}(undef)
    i = 1
    for vertex in verts
        temp[1:2,i] = vertex
        temp[3,i] = 0
        i += 1
    end
    return SMatrix{3,M,T}(temp)
 end

"""Returns a number between 0 and 1 representing the ratio of the lens radiance to the pupil radiance. Assume lᵣ is the radiance transmitted through the lens from the display point. Some of this radiance, pᵣ, passes through the pupil. The beam energy is the ratio pᵣ/lᵣ."""
function beamenergy(assy::LensletAssembly{T},displaypoint::AbstractVector{T},pupilpoints::SMatrix{3,N,T}) where{N,T<:Real}
    llens::ParaxialLens{T} = lens(assy)
    virtpoint = point(virtualpoint(llens,displaypoint))
    projectedpoints = project(assy,displaypoint,pupilpoints)
    beampupil = SphericalPolygon(pupilpoints,virtpoint,T(1))

    intsct = LazySets.VPolygon(projectedpoints) ∩ LazySets.VPolygon(vertices(llens)) #this could be slow, especially multithreaded, because it will allocate. Lazysets.vertices returns Vector{SVector}, rather than SMatrix or SVector{SVector}.

    if isempty(intsct)
        return T(0)
    else
        beamintsct = SphericalPolygon(convertlazysets(LazySets.vertices(intsct)),virtpoint,T(1))
        return OpticSim.area(beamintsct)/OpticSim.area(beampupil)
    end
end

"""Compute the bounding box in the display plane of the image of the worldpoly on the display plane. Pixels inside this area, conservatively, need to be turned on because some of their rays may pass through the worldpoly"""
function activebox(lenslet::LensletAssembly{T},worldpoly::SVector{N,SVector{3,T}}) where{T<:Real,N}
    maxx = typemin(T)
    maxy = typemin(T)
    minx = typemax(T)
    miny = typemax(T)
    lens = lens(lenslet)
    display = display(lenslet)

    for point in worldpoly
        locpoint = transform(lenslet)*point

        for lenspoint in vertices(lens)
            ray = Ray(locpoint,lenspoint-locpoint)
            refractedray,_,_ = processintersection(interface(lens),lenspoint,ray, T(OpticSim.GlassCat.TEMP_REF), T(OpticSim.GlassCat.PRESSURE_REF))
            displaypoint = surfaceintersection(display,refractedray)
            maxx = max(maxx,displaypoint[1])
            minx = min(minx,displaypoint[1])
            maxy = max(maxy,displaypoint[2])
            miny = min(miny,displaypoint[2])
        end
    end
    return minx,miny,maxx,maxy
end


text:
using ConfParser, Compat
using Compat.Test

conf = ConfParse("confs/config.simple")
parse_conf!(conf)

@test retrieve(conf, "protocol") == "kreberos"
@test retrieve(conf, "port")     == "6643"
@test erase!(conf, "protocol")   == true
save!(conf, "confs/out.conf")

conf = ConfParse("confs/out.conf")
parse_conf!(conf)
@test retrieve(conf, "port")     == "6643"


text:
"""Perform the dataset `xs` separation by thresholding a dataset `bins`ed histogram based on the algorith `T`"""
function separation(::Type{T}, xs::Vector{S};
                    bins = 20, tol::Real = 1.0e-5) where {T <: Thresholding, S <: Real}
    minX, maxX = extrema(xs)

    # get histogram
    r = range(minX, stop=maxX, length=bins+1)
    H = fit(Histogram, xs, r, closed=:left)
    N = length(H.weights)
    @logmsg DEBUG_SEPARATION "Separation histogram" H

    return try
        thresh = fit(T, H; tol=tol)

        minIdx = minimum(thresh)
        s = stats(thresh)
        threshold = minX + (minIdx)*(maxX - minX)/N
        # println(s[minIdx,:])
        discriminability = (abs(s[minIdx,3]-s[minIdx,4]))/(sqrt(s[minIdx,5]+s[minIdx,6]))

        Separation(convert(S, depth(thresh)), convert(S, discriminability), threshold, minIdx, minX, maxX, bins)
    catch ex
        # rethrow(ex)
        @logmsg TRACE "Separation error" ex
        Separation(zero(S), zero(S), maxX, 0, minX, maxX, bins) # threshold is set to maximal element
    end
end

"""Kittler algorithm thresholding type"""
struct Kittler <: Thresholding
    statistics::Matrix{Float64}
    depth::Float64
    minindex::Int
end
depth(t::Kittler) = t.depth
stats(t::Kittler) = t.statistics
Base.minimum(t::Kittler) = t.minindex

"""Performs Kittler's minimal thresholding algorithm

J. Kittler & J. Illingworth: "Minimum Error Thresholding", Pattern Recognition, Vol 19, nr 1. 1986, pp. 41-47.
"""
function fit(::Type{Kittler}, H::Histogram; tol::Real=1.0e-5)
    N = length(H.weights)

    # calculate stats
    S = histstats(H.weights./sum(H.weights))

    # Compute criterion function
    J = fill(2*log(eps()), N-1)
    for t=1:(N-1)
        @logmsg TRACE "Statistics" P₁=S[t,1] P₂=S[t,2] σ₁=S[t,5] σ₂=S[t,6]
        if S[t,1]!=0 && S[t,2]!=0 && S[t,5]>0 && S[t,6]>0
            σ₁ = sqrt(max(0, S[t,5]))
            logσ₁ = log(σ₁ == 0 ? eps() : σ₁)
            σ₂ = sqrt(max(0, S[t,6]))
            logσ₂ = log(σ₂ == 0 ? eps() : σ₂)
            ses = S[t,1]*logσ₁ + S[t,2]*logσ₂
            se = S[t,1]*log(S[t,1]) + S[t,2]*log(S[t,2])
            J[t] = 1 + 2*ses - 2*se
        end
    end
    @logmsg DEBUG_SEPARATION "Kittler criterion function" J=J' extrema=extrema(J)

    # Global minimum parameters
    depth, global_min = find_global_min(J, tol)

    return Kittler(S, depth, round(Int, global_min))
end

"""Otsu algorith thresholding type"""
struct Otsu <: Thresholding
    statistics::Matrix{Float64}
    depth::Float64
    minindex::Int
end
depth(t::Otsu) = t.depth
stats(t::Otsu) = t.statistics
Base.minimum(t::Otsu) = t.minindex

"""Performs Otsu thresholding algorithm

N. Otsu: "A threshold selection method from gray-level histograms", Automatica, 1975, 11, 285-296.
"""
function fit(::Type{Otsu}, H::Histogram; tol=1.0e-5)
    N = length(H.weights)

    # calculate stats
    S = histstats(H.weights./sum(H.weights))

    # Compute criterion function
    J = zeros(N-1)
    for t in 1:N-1
        σb = S[t,1]*S[t,2]*(S[t,4]-S[t,3])^2
        σw = S[t,1]*S[t,5]^2 - S[t,2]*S[t,6]^2
        J[t] = σb/σw
    end
    @logmsg DEBUG_SEPARATION "Otsu criterion function" J=J'

    # Global minimum parameters
    depth, global_min = find_global_min(J, tol)

    return Otsu(S, depth, round(Int, global_min))
end


text:
using Roots

export waveNumber, waveNumberHunt

# Dispersion relation for water waves (propagating mode)
dispersion_rel(k,m0,h)=m0-k*tanh(k*h);

# First derivative of dispersion relation for water waves (propagating mode)
dispersio_rel_1st_der(k,m0,h)=-(m0+h*(k^2-m0^2))/k;

# Dispersion relation for water waves (evanescent modes)
dispersion_rel_evan(k,m0,h) = m0 + k*sin(k*h)/cos(k*h);

# First derivative of dispersion relation for water waves (evanescent modes)
dispersion_rel_evan_1st_der(k,m0,h) = -( m0 - h*(k*k + m0*m0) )/k;

function waveNumber(imode,m0::T,h::T) where T
    if imode==0
       f(x) = dispersion_rel(x,m0,h)
       df(x) = dispersio_rel_1st_der(x,m0,h)
       ki=Roots.newton(f, df, m0);
       #ki = Roots.find_zero(f,m0,Order0())
    else
       fa(x) = dispersion_rel_evan(x,m0,h)
       dfa(x) = dispersion_rel_evan_1st_der(x,m0,h)
       ki = Roots.newton(fa, dfa,imode*pi/h);
    end
    return ki;
end

function waveNumberHunt(m0::T,h::T) where T
   d=[0.6666666666,0.3555555555,0.1608465608,0.0632098765,0.0217540484,0.0065407983]
   y=m0*h
   temp=0.
   for i in range(1,stop=6)
      temp+=d[i]*y^i
   end
   return sqrt(y^2+y/(1+temp))/h
end


text:
# ---
# title: 1105. Filling Bookcase Shelves
# id: problem1105
# author: Indigo
# date: 2022-03-03
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/filling-bookcase-shelves/description/>
# hidden: true
# ---
#
# We have a sequence of `books`: the `i`-th book has thickness `books[i][0]` and
# height `books[i][1]`.
#
# We want to place these books **in order**  onto bookcase shelves that have
# total width `shelf_width`.
#
# We choose some of the books to place on this shelf (such that the sum of their
# thickness is `<= shelf_width`), then build another level of shelf of the
# bookcase so that the total height of the bookcase has increased by the maximum
# height of the books we just put down.  We repeat this process until there are
# no more books to place.
#
# Note again that at each step of the above process, _the order of the books we
# place is the same order as the given sequence of books_.  For example, if we
# have an ordered list of 5 books, we might place the first and second book onto
# the first shelf, the third book on the second shelf, and the fourth and fifth
# book on the last shelf.
#
# Return the minimum possible height that the total bookshelf can be after
# placing shelves in this manner.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2019/06/24/shelves.png)
#
#
#
#     Input: books = [[1,1],[2,3],[2,3],[1,1],[1,1],[1,1],[1,2]], shelf_width = 4
#     Output: 6
#     Explanation:
#     The sum of the heights of the 3 shelves are 1 + 3 + 2 = 6.
#     Notice that book number 2 does not have to be on the first shelf.
#
#
#
#
# **Constraints:**
#
#   * `1 <= books.length <= 1000`
#   * `1 <= books[i][0] <= shelf_width <= 1000`
#   * `1 <= books[i][1] <= 1000`
#
#
## @lc code=start
using LeetCode

function min_height_shelves(books::Vector{Vector{Int}}, shelf_width::Int)
    n = length(books)
    dp = OffsetArray(fill(typemax(Int) >> 1, n + 1), -1)
    dp[0] = 0
    for i in 1:n
        idx, width, height = i - 1, books[i][1], books[i][2]
        while width <= shelf_width
            dp[i] = min(dp[i], dp[idx] + height)
            idx == 0 && break
            height = max(height, books[idx][2])
            width += books[idx][1]
            idx -= 1
        end
    end
    dp[end]
end
## @lc code=end


text:
using Primes;
using DelimitedFiles;

function rand1024()
  x::BigInt=BigInt(1)
  for i=1:16
    x=x<<64+rand(UInt64)
  end
  return x
end

function randRSAPrime()
  x::BigInt=0
  while true
    x = rand1024()
    if !isprime(x)
      continue
    end
    if isprime(div(x,2))
      break
    end
  end
  return x
end

writedlm("divisor.txt", [randRSAPrime()*randRSAPrime()])


text:
using ReinforcementLearningEnvironmentAtari, Test, ReinforcementLearningBase

test_envinterface(AtariEnv("pong"))

env = AtariEnv("pong", frame_skip = 10000)
@test ReinforcementLearningEnvironmentAtari.interact!(env, 1).isdone == true
env = AtariEnv("pong", frame_skip = 1)
@test ReinforcementLearningEnvironmentAtari.interact!(env, 1).isdone == false


text:
using Plots

theme(:mydark)   # or any other dark theme

f(x) = sin(10-x) * cos((10-x)/10) * exp(-(10-x)/10)
f1(x) = 0.5*sin(2x) + 0.3*sin(1.4x+5) - 0.27
f2(x) = 1.27*sinc(1.3x-3.5) - 0.42

let p = plot(f, 0, 10; size=(1300,500), color=RGBA(1,0,0,0.5), legend=false)
  for a ∈ 1.0:-0.001:0.7
    plot!(p, x -> a*f(x), 0, 10; color=RGBA(1,0,0,0.03), legend=false)
  end
  plot!(p, f2, 0, 10; color=RGBA(0,1,0,0.25), legend=false)
  for a ∈ 1.0:-0.001:0.7
    plot!(p, x -> a*f2(x), 0, 10; color=RGBA(0,1,0,0.01), legend=false)
  end
  plot!(p, f1, 0, 10; color=RGBA(0,0,1,0.5), legend=false)
  for a ∈ 1.0:-0.001:0.7
    plot!(p, x -> a*f1(x), 0, 10; color=RGBA(0,0,1,0.03), legend=false)
  end
  scatter!(p, [2.7], [0.5]; markersize=10, markercolor=RGB(0.220,0.596,0.149), markeralpha=0.5, markerstrokecolor=:black, legend=false)
  scatter!(p, [2.8], [0.4]; markersize=10, markercolor=RGB(0.584,0.345,0.698), markeralpha=0.5, markerstrokecolor=:black, legend=false)
  scatter!(p, [2.6], [0.4]; markersize=10, markercolor=RGB(0.796,0.235,0.200), markeralpha=0.5, markerstrokecolor=:black, legend=false)
  p
end

max_stars_repo_name: org-arl/unet-blog

max_issues_repo_name: org-arl/unet-blog

max_forks_repo_name: org-arl/unet-blog


text:
relaxed_tols = (default_tol_relax = 1000,)
insts = OrderedDict()
insts["minimal"] = [
    ((2, 4, true, true),),
    ((2, 4, true, false),),
    ((2, 2, false, true),),
    ((2, 2, false, false),),
    ]
insts["fast"] = [
    ((3, 6, true, true),),
    ((3, 6, true, false),),
    ((5, 5, false, true),),
    ((5, 5, false, false),),
    ((10, 20, true, true),),
    ((10, 20, true, false),),
    ((15, 15, false, true),),
    ((15, 15, false, false),),
    ((30, 30, false, false),),
    ((30, 30, false, true),),
    ]
insts["various"] = [
    ((6, 6, true, true), nothing, relaxed_tols),
    ((6, 6, true, false), nothing, relaxed_tols),
    ((6, 6, false, true),),
    ((6, 6, false, false),),
    ((12, 12, true, true), nothing, relaxed_tols),
    ((12, 12, false, true),),
    ((12, 12, false, false),),
    ((24, 24, true, true), nothing, relaxed_tols),
    ((24, 24, false, true), nothing, relaxed_tols),
    ((24, 24, false, false), nothing, relaxed_tols),
    ]
return (LyapunovStabilityJuMP, insts)


text:
"""
$(TYPEDEF)

The solution to a `DataDrivenProblem` derived via a certain algorithm.
The solution is represented via an `AbstractBasis`, which makes it callable.

# Fields
$(FIELDS)

# Note

The L₂ norm error, AIC and coefficient of determinantion get only computed, if eval_expression is set to true or
if the solution can be interpreted as a linear regression result.
"""
struct DataDrivenSolution{L, A, O} <: AbstractDataDrivenSolution
    "The basis representation of the solution"
    basis::AbstractBasis
    "Parameters of the solution"
    parameters::AbstractVecOrMat
    "Returncode"
    retcode::Symbol
    "Algorithm"
    alg::A
    "Original output of the solution algorithm"
    out::O
    "Problem"
    prob::AbstractDataDrivenProblem
    "L₂ norm error"
    l2_error::AbstractVector
    "AIC"
    aic::AbstractVector
    "Coefficient of determinantion"
    rsquared::AbstractVector

    function DataDrivenSolution(linearity::Bool,b::AbstractBasis, p::AbstractVector, retcode::Symbol, alg::A, out::O, prob::AbstractDataDrivenProblem; kwargs...) where {A,O}
        return new{linearity, A,O}(
            b, p, retcode, alg, out, prob
        )
    end

    function DataDrivenSolution(linearity::Bool,b::AbstractBasis, p::AbstractVector, retcode::Symbol, alg::A, out::O, prob::AbstractDataDrivenProblem, l2::AbstractVector; kwargs...) where {A,O}
        return new{linearity, A,O}(
            b, p, retcode, alg, out, prob, l2
        )
    end

    function DataDrivenSolution(linearity::Bool,b::AbstractBasis, p::AbstractVector, retcode::Symbol, alg, out, prob::AbstractDataDrivenProblem, l2::AbstractVector, aic::AbstractVector; kwargs...)
        return new{linearity, typeof(alg), typeof(out)}(
            b, p, retcode, alg, out, prob, l2, aic
        )
    end

    function DataDrivenSolution(b::AbstractBasis, p::AbstractVector, retcode::Symbol, alg, out, prob::AbstractDataDrivenProblem, l2::AbstractVector, aic::AbstractVector, rsquared::AbstractVector; kwargs...)
        return new{true, typeof(alg), typeof(out)}(
            b, p, retcode, alg, out, prob, l2, aic, rsquared
        )
    end

end

function DataDrivenSolution(b::AbstractBasis, p::AbstractVector, retcode::Symbol, alg::A, out::O, prob::AbstractDataDrivenProblem, linearity::Bool;
    eval_expression = false, kwargs...) where {A,O}

    if !eval_expression
        # Compute the errors

        x, _, t, u = get_oop_args(prob)
        y = get_target(prob)

        if is_implicit(b)
            e = b([x; y], p, t, u)
        else
            e = get_target(prob) - b(x, p, t, u)
        end

        l2 = sum(abs2, e, dims = 2)[:,1]
        aic = 2*(-size(e, 2) .* log.(l2 / size(e, 2)) .+ length(p))

        if linearity
            rsquared = 1 .- mean(e, dims = 2)[:,1] ./ var(get_target(prob), dims = 2)[:,1]
            #return l2, aic, rsquared
            return DataDrivenSolution(
                b, p, retcode, alg, out, prob, l2, aic, rsquared
            )
        end

        return DataDrivenSolution(
            linearity, b, p, retcode, alg, out, prob, l2, aic
        )
    end

    return DataDrivenSolution(
        linearity, b, p, retcode, alg, out, prob
    )
end

## Make it callable
(r::DataDrivenSolution)(args...) = r.basis(args...)

function Base.show(io::IO, ::DataDrivenSolution{linearity}) where linearity
    if linearity
        print(io, "Linear Solution")
    else
        print(io, "Nonlinear Solution")
    end
end

function Base.print(io::IO, r::DataDrivenSolution{linearity, a, o}) where {linearity, a, o}
    show(io, r)
    print(io, " with $(length(r.basis)) equations and $(length(r.parameters)) parameters.\n")
    print(io, "Returncode: $(r.retcode)\n")
    isdefined(r, :l2_error) && print(io, "L₂ Norm error : $(r.l2_error)\n")
    isdefined(r, :aic) && print(io,"AIC : $(r.aic)\n")
    isdefined(r, :rsquared) && print(io, "R² : $(r.rsquared)\n")
    return
end

function Base.print(io::IO, r::DataDrivenSolution, fullview::DataType)

    fullview != Val{true} && return print(io, r)

    print(io, r)

    if length(r.parameters) > 0
        x = parameter_map(r)
        println(io, "Parameters:")
        for v in x
            println(io, "   $(v[1]) : $(v[2])")
        end
    end

    return
end

##
"""
$(SIGNATURES)

Returns the `Basis` of the result.
"""
result(r::DataDrivenSolution) = r.basis

"""
$(SIGNATURES)

Returns the AIC of the result.
"""
aic(r::DataDrivenSolution) = begin
    isdefined(r, :aic) && return r.aic
    return NaN
end

"""
$(SIGNATURES)

Returns the L₂ norm error of the result.
"""
error(r::DataDrivenSolution) = begin
    isdefined(r, :l2_error) && return r.l2_error
    return NaN
end

"""
$(SIGNATURES)

Returns the coefficient of determinantion of the result, if the result has been
derived via a linear regression, e.g. sparse regression or koopman.
"""
determination(r::DataDrivenSolution{l, o, a}) where {l,o, a} = begin
    if l
        return r.rsquared
    else
        return NaN
    end
end

"""
$(SIGNATURES)

Returns the original `DataDrivenProblem`.
"""
get_problem(r::DataDrivenSolution) = getfield(r, :prob)

"""
$(SIGNATURES)

Returns the estimated parameters in form of an `Vector`.
"""
ModelingToolkit.parameters(r::DataDrivenSolution) = r.parameters

"""
$(SIGNATURES)

Generate a mapping of the parameter values and symbolic representation useable
to `solve` and `ODESystem`.
"""
function parameter_map(r::DataDrivenSolution)
    return [
        ps_ => p_ for (ps_, p_) in zip(parameters(r.basis), r.parameters)
    ]
end

"""
$(SIGNATURES)

Returns the algorithm used to derive the solution.
"""
algorithm(r::DataDrivenSolution) = r.alg

"""
$(SIGNATURES)

Returns the original output of the algorithm.
"""
output(r::DataDrivenSolution) = r.out

"""
$(SIGNATURES)

Returns all applicable metrics of the solution.
"""
function metrics(r::DataDrivenSolution)
    fnames_ = (:l2_error, :aic, :rsquared)
    names_ = (:L₂, :AIC, :R²)
    m = Dict()
    for i in 1:length(fnames_)
        if isdefined(r, fnames_[i])
            push!(m, names_[i] => getfield(r, fnames_[i]))
        end
    end
    m
end
## Helper for the solution

# Check linearity

function assert_linearity(eqs::AbstractVector{Equation}, x::AbstractVector)
    return assert_linearity(map(x->Num(x.rhs), eqs), x)
end

# Returns true iff x is not in the arguments of the jacobian of eqs
function assert_linearity(eqs::AbstractVector{Num}, x::AbstractVector)
    j = Symbolics.jacobian(eqs, x)
    # Check if any of the variables is in the jacobian
    v = unique(reduce(vcat, map(get_variables, j)))
    for xi in x, vi in v
        isequal(xi, vi) && return false
    end
    return true
end

function construct_basis(X, b, implicits = Num[]; dt = one(eltype(X)), lhs::Symbol = :continuous, is_implicit = false, eval_expression = false)

    # Create additional variables
    sp = sum(.! iszero.(X))
    sps = norm.(eachcol(X), 0)
    inds = sps .> zero(eltype(X))
    pl = length(parameters(b))

    p = [Symbolics.variable(:p, i) for i in (pl+1):(pl+sp)]
    p = collect(p)
    ps = zeros(eltype(X), sp)

    eqs = zeros(Num, sum(inds))
    eqs_ = [e.rhs for e in equations(b)]
    cnt = 1
    for j in 1:size(X, 2)
        if sps[j] == zero(eltype(X))
            continue
        end
        for i in 1:size(X, 1)
            if iszero(X[i,j])
                continue
            end
            ps[cnt] = X[i,j]
            eqs[j] += p[cnt]*eqs_[i]
            cnt += 1
        end
    end

    # Build the lhs
    xs = states(b)
    if isempty(implicits) || !is_implicit
        if length(eqs) == length(states(b))
            if lhs == :continuous
                d = Differential(get_iv(b))
                eqs = [d(xs[i]) ~ eq for (i,eq) in enumerate(eqs)]
            elseif lhs == :discrete
                d = Difference(get_iv(b), dt = dt)
                eqs = [d(xs[i]) ~ eq for (i,eq) in enumerate(eqs)]
            end
        end
    else
        if !isempty(implicits)
            eqs = eqs .~ 0
            if assert_linearity(eqs, implicits)
                try
                    # Try to solve the eq for the implicits
                    eqs = ModelingToolkit.solve_for(eqs, implicits)
                    eqs = implicits .~ eqs
                    implicits = []
                catch
                    @warn "Failed to solve recovered equations for implicit variables. Returning implicit equations."
                end
            end
        end
    end

    Basis(
        eqs, xs,
        parameters = [parameters(b); p], iv = get_iv(b),
        controls = controls(b), observed = observed(b),
        implicits = implicits,
        name = gensym(:Basis),
        eval_expression = eval_expression
    ), ps
end

function _round!(x::AbstractArray{T, N}, digits::Int) where {T, N}
    for i in eachindex(x)
        x[i] = round(x[i], digits = digits)
    end
    return x
end

function assert_lhs(prob)
    dt = mean(diff(prob.t))
    lhs = :direct
    if isa(prob, AbstracContProb)
        lhs = :continuous
    elseif isa(prob, AbstractDiscreteProb)
        lhs = :discrete
    else
        lhs = :direct
    end
    return lhs, dt
end

function DataDrivenSolution(prob, s, b::B, opt::O, implicits = Num[];
    eval_expression = false, digits::Int = 10, by = :min, kwargs...) where {B <: AbstractBasis, O <: AbstractOptimizer}
    Ξ, _... = select_by(by, s)
    # Build a basis and returns a solution
    if all(iszero.(Ξ))
        @warn "Sparse regression failed! All coefficients are zero."
        return DataDrivenSolution(
        b, [], :failed, opt, Ξ, prob)
    end

    # Assert continuity
    lhs, dt = assert_lhs(prob)

    sol , ps = construct_basis(round.(Ξ, digits = digits), b, implicits,
        lhs = lhs, dt = dt,
        is_implicit = isa(opt, AbstractSubspaceOptimizer) ,eval_expression = eval_expression
        )

    ps = isempty(parameters(b)) ? ps : vcat(prob.p, ps)

    return DataDrivenSolution(
        sol, ps, :solved, opt, s, prob, true, eval_expression = eval_expression
    )
end

function DataDrivenSolution(prob, k, b::BS, alg::KA;
    digits::Int = 10, by = :min, eval_expression = false, operator_only = false, kwargs...) where {BS <: AbstractBasis, KA <: AbstractKoopmanAlgorithm}
    k_, _... = select_by(by, k)
    @unpack inds = k
    K, B, C, P, Q = k_
    operator_only && return (K = K, B = B, C = C, P = P, Q = Q)
    # Build parameterized equations, inds indicate the location of basis elements containing an input
    Ξ = zeros(eltype(C), size(C,2), length(b))

    Ξ[:, inds] .= real.(Matrix(K))
    if !isempty(B)
        Ξ[:, .! inds] .= B
    end

    # Assert continuity
    lhs, dt = assert_lhs(prob)

    bs, ps = construct_basis(round.(C*Ξ, digits = digits)', b,
        lhs = lhs, dt = dt,
        eval_expression = eval_expression)

    res_ = Koopman(equations(bs), states(bs),
        parameters = parameters(bs),
        controls = controls(bs), iv = get_iv(bs),
        K = K, C = C, Q = Q, P = P, lift = get_f(b),
        is_discrete = is_discrete(prob),
        eval_expression = eval_expression)

    ps = isempty(parameters(b)) ? ps : vcat(prob.p, ps)

    return DataDrivenSolution(
        res_, ps, :solved, alg, k, prob, true, eval_expression = eval_expression
    )
end


text:
# This code changes behavior in LinearAlgebra.Diagonal
# It can't go to FillArrays

function sum(f::Function, D::LinearAlgebra.Diagonal{T}) where T
    m = size(D, 1)
    return (m * (m - 1)) * f(zero(T)) + sum(f, D.diag)
end
sum(D::LinearAlgebra.Diagonal) = sum(identity, D)

function prod(f::Function, D::LinearAlgebra.Diagonal{T}) where T
    m = size(D, 1)
    return f(zero(T))^(m * (m - 1)) * prod(f, D.diag)^m
end
prod(D::LinearAlgebra.Diagonal) = prod(identity, D)

function minimum(D::Diagonal{T}) where T <: Number
    mindiag = minimum(D.diag)
    size(D, 1) > 1 && return (min(zero(T), mindiag))
    return mindiag
end

function maximum(D::Diagonal{T}) where T <: Number
    maxdiag = maximum(D.diag)
    size(D, 1) > 1 && return (max(zero(T), maxdiag))
    return maxdiag
end

any(f::Function, x::Diagonal{T}) where T <: Number = size(x, 1) == 1 ? f(x[1]) :
    f(zero(T)) || any(f, x.diag)

all(f::Function, x::Diagonal{T}) where T <: Number = size(x, 1) == 1 ? f(x[1]) :
    f(zero(T)) && all(f, x.diag)

function kron(A::AbstractMatrix{T}, B::Diagonal{S}) where {T<:Number, S<:Number}
    @assert ! has_offset_axes(A)
    (mA, nA) = size(A); (mB, nB) = size(B)
    R = zeros(promote_op(*, T, S), mA * mB, nA * nB)
    m = 1
    for j = 1:nA
        for l = 1:mB
            Bll = B[l,l]
            for k = 1:mA
                R[m] = A[k,j] * Bll
                m += nB
            end
            m += 1
        end
        m -= nB
    end
    return R
end

# Broadcast.broadcasted(f::typeof(*), x::Number, D::Diagonal{<:Number}) = Diagonal(Broadcast.broadcast(f, x, D.diag))
# Broadcast.broadcasted(f::typeof(*), D::Diagonal{<:Number}, x::Number) = Broadcast.broadcasted(f, x, D)

# for ft in (round, floor, abs, sqrt, sin, sind, cbrt, tan)
#     @eval Broadcast.broadcasted(f::typeof($ft), D::Diagonal{<:Number}) = Diagonal(broadcast(f, D.diag))
# end

# function Broadcast.broadcasted(f::Function, D::Diagonal{T}) where T <: Number
#     fz = f(zero(T))
#     return iszero(fz) ? broadcasted_diag_zero(f, D) : broadcasted_diag_full(f, fz, D)
# end

# @inline function broadcasted_diag_zero(f, D)
#     return Diagonal(broadcast(f, D.diag))
# end

# @inline function broadcasted_diag_full(f, fz, D::Diagonal{T}) where T
#     R = fill(fz, size(D))
#     for i in 1:size(D, 1)
#         R[i,i] = f(D.diag[i])
#     end
#     return R
# end


text:
function free(cur_player::AbstractPlayer)
  cur_player.action = ""

  cur_player.clock.hold =
    Nullable{Base.Random.UUID}()

  return
end


text:
using ArgParse
import Reactive, GeometryTypes, Chipmunk, LinearAlgebra
const re = Reactive
const gt = GeometryTypes
const Cp = Chipmunk
const la = LinearAlgebra

"""
    Config

Saves all data, that is need through the whole programm. Will be initialized at
startup.
"""
mutable struct Config{WS<:WorldState, T<:Number}
    ws::WS
    """
    Field:

    Rectangular, (0,0) is lower left corner
    (9,6) is the right upper corner
    all in meters
    """
    field_size::Vector{T}
    field_corners::Vector{<:gt.Point2{T}}
    active_strategy::String
    strategies_table::Dict{String, Module}
    maximum_player_per_side::Int
    space::Cp.Space
    # set_progress_bar::Function
"""
    printables

For values that shall be visible in the gui. The name will be a first label and the value will be second label, as a signal. If it changes, the value in the second label will also be updated.
"""
    printables::Dict{String, re.Signal}

"""
    drawables

A dictionary for drawing something on tha canvas in the gui (where the field is displayed).
The name uses the folling conventions for the drawing:
- if the name beginns with '_' it will be drawn befor anything on the field is drawn and it will be assumed, that it is an image (matrix), the position is the top left corner and the botel right corner
- in all other cases it must be a element from Graphics2D, as function, that recives the first argument (the position). In the gui the posiiton will be transformed from world coordinates to screen coordinates
"""
    drawables::Dict{String, <:Tuple{<:AbstractVector{<:AbstractVector{<:Real}}, Function}} # the tuples saves the coordinate in the first place and a drawable function in the second

"""
    manipulators

Saves Gtk gui-elements. The name will be printed in a label above the elemente.
Name conventions:
- if the name starts with a D: instead of a widget as the element a tuple with a widget and a function is assumed. The function will be added as "button-pressed" event of the widget and after the execution of the function, the canvas will be redrawn.
- if the name starts with DN: the same as D, but the widget will not be shown in the widget area
- if the widget is a Gtk.GtkRadioButtonGroupLeaf and the are less than two RadioButtons in it, it will not be shown
- all others are names+widgets and the widget will be shown
"""
    manipulators::Dict{String, <:Tuple}

    global_communication::Dict{String, Any}
end
Config(ws::WorldState, field_size::Vector{<:AbstractFloat},
    field_corners::Vector{<:gt.Point2{<:Number}},
    active_strategy::String, space::Cp.Space) = Config(ws, field_size, field_corners,
    active_strategy, Dict{String, Module}(), 6, space,
    Dict{String, re.Signal}(),
    Dict{String, Tuple{Vector{Vector{Float64}}, Function}}(),
    Dict{String, Tuple}(),
    Dict{String, Any}())
function Config(config::Config)
    ws = WorldState(config.ws)
    space = createspace(ws)

    config = Config(ws, config.field_size, config.field_corners, config.active_strategy, space)

    setstrategies!(config)

    config
end

"""
    ArgParse.parse_item(::Type{Robot{Float64}}, x::AbstractString)

A Robot has its own Arguments. A string `x` containing the arguments are parsed to a `Robot` object.
"""
function ArgParse.parse_item(::Type{Robot{Float64}}, x::AbstractString)
    vec = [substring for substring in split(x[2:end-1], [' ']) if substring != ""]
    default = Nao([0.0, 0.0], 0.0)
    s = ArgParseSettings(commands_are_required = false)
    @add_arg_table s begin
        "--position"
            required = true
            help = "Position in world coordinate system (metric)."
            # range_tester = x -> x > 0
            nargs = 2
            arg_type = Float64
        "--body_orientation"
            help = "Robot orientation relative to the x-axis in rad."
            range_tester = x -> -π ≤ x ≤ π
            arg_type = Float64
            default = default.body_orientation
        "--head_orientation"
            help = "Head orientation relative to the `body_orientation` in rad."
            range_tester = x -> deg2rad(-119.5) ≤ x ≤ deg2rad(119.5)
            arg_type = Float64
            default = default.head_orientation
        "--velocity"
            help = "The current velocity of the robot in m/s."
            # range_tester = x -> -π ≤ x ≤ π
            nargs = 2
            arg_type = Float64
            default = get_velocity(default)
        "--angularvelocity"
            help = "The current angular velocity of the robot rad/s."
            # range_tester = x -> x ≥ 0
            arg_type = Float64
            default = get_angularvelocity(default)
        "--size"
            help = "Robot size in m."
            range_tester = x -> x > 0
            arg_type = Float64
            default = get_size(default)
        "--mass"
            help = "Robot mass in kg."
            range_tester = x -> x ≥ 0
            arg_type = Float64
            default = get_mass(default)
        "--hfov"
            help = "Horizontal field of view in rad."
            range_tester = x -> x ≤ 2π
            arg_type = Float64
            default = default.hfov
        "--max_velocity"
            help = "The maximum velocity of the robot in m/s."
            range_tester = x -> x ≥ 0
            arg_type = Float64
            default = default.max_velocity
        "--max_angularvelocity"
            help = "The maximum angular velocity of the robot in rad/s."
            range_tester = x -> x ≥ 0
            arg_type = Float64
            default = default.max_angularvelocity
        # "--shape"
        #     help = "Shape of the robot."
        #     arg_type = Symbol
        #     default = default.shape
    end
    res = parse_args(vec, s, as_symbols=true)
    val = Robot(;res...)
    return val
end

"""
    ArgParse.parse_item(::Type{Obstacle{Float64}}, x::AbstractString)

A Robot has its own Arguments. A string `x` containing the arguments are parsed to a `Robot` object.
"""
function ArgParse.parse_item(::Type{Obstacle{Float64}}, x::AbstractString)
    vec = [substring for substring in split(x[2:end-1], [' ']) if substring != ""]
    default = Ball([0.0, 0.0])
    s = ArgParseSettings(commands_are_required = false)
    @add_arg_table s begin
        "--position"
            required = true
            help = "Position in world coordinate system (metric)."
            # range_tester = x -> x > 0
            nargs = 2
            arg_type = Float64
        "--velocity"
            help = "The current velocity of the robot in m/s."
            # range_tester = x -> -π ≤ x ≤ π
            nargs = 2
            arg_type = Float64
            default = get_velocity(default)
        "--size"
            help = "Obstacle size in m."
            range_tester = x -> x > 0
            arg_type = Float64
            default = get_size(default)
        "--mass"
            help = "Obstacle mass in kg."
            range_tester = x -> x ≥ 0
            arg_type = Float64
            default = get_mass(default)
        "--max_velocity"
            help = "The maximum velocity of the robot in m/s."
            range_tester = x -> x ≥ 0
            arg_type = Float64
            default = default.max_velocity
        "--friction"
            help = "Friction depending on the ground."
            range_tester = x -> x ≥ 0
            arg_type = Float64
            default = get_friction(default)
        "--fixed"
            help = "If the obstacle can't be moved."
            nargs = 0
            action = :store_true
        # "--shape"
        #     help = "Shape of the robot. Allowed is circle, poly, box and segment.0"
        #     arg_type = Symbol
        #     default = default.shape
        "--color"
            help = "Color of the object."
            range_tester = x -> 0 ≤ x ≤ 1
            nargs = 4
            arg_type = Float64
            default = default.color
    end
    res = parse_args(vec, s, as_symbols=true)
    val = Obstacle(;res...)
    return val
end

"""
    parse_commandline(args::AbstractVector{<:AbstractString})

Parses a string to a list with all command line arguments. If the command line
argument --file is given, the file at this location is parsed with this
function. All other command line arguments in the first place are rejected.
"""
function parse_commandline(args::AbstractVector{<:AbstractString})
    s = ArgParseSettings(commands_are_required = false)

    @add_arg_table s begin
        "--field_size"
            help = "Defines the width and height of the field in meters."
            nargs = 2
            arg_type = Float64
            default = [9.0, 6.0]
        "--add_obstacle"
            help = "Add as many obstacles as you want."
            action = :append_arg
            arg_type = Obstacle{Float64}
            default = Obstacle[]
        "--add_myPlayer", "-p"
            help = "Add as many myPlayers as you want."
            action = :append_arg
            arg_type = Robot{Float64}
            default = Robot[]
        "--add_opponent", "-o"
            help = "Add as many opponents as you want."
            action = :append_arg
            arg_type = Robot{Float64}
            default = Robot[]
        "--gravity"
            help = "Gravity for your location. Default should match for the most areas on earth."
            # range_tester = x -> x ≥ 0
            arg_type = Float64
            default = 9.81
        "--strategy", "-s"
            help = "The strategy, which shall be used."
            arg_type = String
            default = "ShowShadow"
        "--file", "-f"
            help = "Path to a file with all parameters."
            arg_type = String
    end
    arguments = parse_args(args, s)
    if haskey(arguments, "file") && arguments["file"] != nothing
        arguments = parse_args(collect(eachline(open(arguments["file"]))), s)
    end
    return arguments
end

function createspace(ws::WorldState)
    space = Cp.Space()

    for object in vcat(ws.myPlayers, ws.opponents, ws.obstacles)
        Cp.set_userdata(object.body, pointer_from_objref(object))
        Cp.add_body(space, object.body)
        Cp.add_shape(space, object.shape)
    end

    # Cp.set_damping(space, 0.1)

    for player in ws.myPlayers
        Cp.set_velocity_update_func(player.body, (b, g, d, t) -> nothing)
    end
    for obs in ws.obstacles

        # μr = get_friction(obs)
        # g = arguments["gravity"]

        function ball_linear_drag(body_ptr::Ptr{Nothing}, gravity::Cp.Vect, damping::Cdouble, dt::Cdouble) # TODO: needs too long
        # // Skip kinematic bodies.
        # if(cpBodyGetType(body) == CP_BODY_TYPE_KINEMATIC) return;
            # deceleration = 0.09
            body = Cp.Body(body_ptr) # TODO: tooks way too much memory
            # ud = unsafe_pointer_to_objref(Cp.get_userdata(body))
            # μr = get_friction(ud)
            # g = 9.81
            current_vel_vect = Cp.get_velocity(body)
            # current_vel = [current_vel_vect.x, current_vel_vect.y]
            cur_vel_len = sqrt(current_vel_vect.x^2 + current_vel_vect.y^2)
            mul(x::T, y::T) where T<:Number = x*y
            vel_len = cur_vel_len - mul(9.81, dt)
            new_vel_x = vel_len > 0 ? current_vel_vect.x/cur_vel_len*vel_len : 0.0
            new_vel_y = vel_len > 0 ? current_vel_vect.y/cur_vel_len*vel_len : 0.0
            Cp.set_velocity(body, Cp.Vect(new_vel_x, new_vel_y))
            return nothing
        end

        Cp.set_velocity_update_func(obs.body, ball_linear_drag)

    end

    return space
end

function setstrategies!(config::Config)
    config.strategies_table["CompareShadow"] = CompareShadow
    config.strategies_table["ShowShadow"] = ShowShadow
    config.strategies_table["TestStrategy"] = TestStrategy
    config.strategies_table["GlobalComShadow"] = GlobalComShadow
    config.strategies_table["GrowingRegions"] = GrowingRegions
    # config.strategies_table["BaMichelle"] = BaMichelle
    # config.strategies_table["BaJakob"] = BaJakob
    config.strategies_table["StartPositions01"] = StartPositions01
    config.strategies_table["PushballNN"] = PushballNN
end

"""
    get_config(args::AbstractVector{<:AbstractString})

Wraps the parse_commandline function and creates a `config` from the command
line list.
"""
function get_config(args::AbstractVector{<:AbstractString})
    arguments = parse_commandline(args) #String["--file", config_file_path]

    player = vcat(arguments["add_myPlayer"], arguments["add_opponent"])
    for i in 1:length(player), j in i+1:length(player)
        if distance(player[i], player[j]) < get_size(player[i]) + get_size(player[j]) # HACK: asumes circles
            a, b = player[i], player[j]
            error("Player on position $a and $b are too close.")
        end
    end

    field_size = arguments["field_size"]
    field_corners = [gt.Point2(-field_size[1]/2, -field_size[2]/2), gt.Point2(field_size[1]/2, -field_size[2]/2),
        gt.Point2(field_size[1]/2, field_size[2]/2), gt.Point2(-field_size[1]/2, field_size[2]/2)]
    ws = WorldState(arguments["add_myPlayer"], arguments["add_opponent"],
        arguments["add_obstacle"])
    active_strategy = arguments["strategy"]

    space = createspace(ws)

    config = Config(ws, field_size, field_corners, active_strategy, space)

    setstrategies!(config)

    config
end

function add_object!(config::Config, object::Object)
    Cp.set_userdata(object.body, pointer_from_objref(object))
    Cp.add_body(config.space, object.body)
    Cp.add_shape(config.space, object.shape)
end
function remove_object!(config::Config, object::Object)
    Cp.remove_body(config.space, object.body)
    Cp.remove_shape(config.space, object.shape)
end

"""
    get_strategy(config::Config)

Wrapper for easier access of the current strategy run function in the config.
"""
function get_strategy(config::Config)
    config.strategies_table[config.active_strategy]
end

function free!(config::Config)
    for obj in vcat(config.ws.myPlayers, config.ws.opponents, config.ws.obstacles)
        Cp.free(obj.shape)
        Cp.free(obj.body)
    end
    Cp.free(config.space)
    return
end


text:
#=
    Exploring Corner transfer matrices ...
    p 16
=#

struct fonedirCTMIterable{T,
        TA <: AbstractTensor{T,4},
        TC <: AbstractTensor{T,2},
        TT <: AbstractTensor{T,3}} <: AbstractCTMIterable
    A::TA
    χ::Int
    Csinit::Union{NTuple{2,TC},Nothing}
    Tsinit::Union{NTuple{3,TT},Nothing}
end

function fonedirctmiterable(A::DTensor{T,4}, χ::Int, Csinit = nothing, Tsinit = nothing) where T
    fonedirCTMIterable{T,DTensor{T,4},DTensor{T,2},DTensor{T,3}}(A, χ, Csinit, Tsinit)
end

function fonedirctmiterable(A::DASTensor{T,4,SYM,CHS,SS,CH}, χ::Int,
        Csinit = nothing, Tsinit = nothing) where {T,N,SYM,CHS,SS,CH}
    TC = DASTensor{T,2,SYM,CHS,SS,CH}
    TT = DASTensor{T,3,SYM,CHS,SS,CH}
    TA = DASTensor{T,4,SYM,CHS,SS,CH}
    fonedirCTMIterable{T,TA,TC,TT}(A, χ, Csinit, Tsinit)
end

struct fonedirCTMState{S,
        TA <: AbstractTensor,
        TC <: AbstractTensor,
        TT <: AbstractTensor} <: AbstractCTMState
    Cs::NTuple{2,TC}
    Ts::NTuple{3,TT}
    oldsvdvals::Vector{S}
    diffs::Vector{S}
    n_it::Ref{Int}
end

function iterate(iter::fonedirCTMIterable{S,TA,TC,TT}) where {S,TA,TC,TT}
    @unpack A, χ, Tsinit, Csinit = iter
    Cs = Csinit == nothing ? ntuple(i -> initializeC(A,χ), 2) : deepcopy(Csinit)
    Ts = Tsinit == nothing ? ntuple(i -> initializeT(A,χ), 3) : deepcopy(Tsinit)

    l = ifelse(TA <: DASTensor, 2χ, χ)
    oldsvdvals = zeros(S,l)
    state = fonedirCTMState{S,TA,TC,TT}(Cs, Ts, oldsvdvals, [], Ref(0))
    return state, state
end

function iterate(iter::fonedirCTMIterable, state::fonedirCTMState{S,TA,TC,TT}) where {S,TA,TC,TT}
    @unpack A, χ = iter
    @unpack Cs, Ts, oldsvdvals, diffs, n_it = state
    C1, C2 = Cs
    T1, T2, T4 = Ts
    #=
        [C1] - [T1] - [C2]
          |      |      |
        [T4] - [ A] - [T2]
          |      |      |
        [C1]*- [T1]*- [C2]*
    =#

    #xmove
    fonedirxmove!(Cs, Ts, A, χ)

    #ymove
    fonedirymove!(Cs, Ts, A, χ)

    @tensor C[1,2] := C1[1,-1] * C2[-1,2]
    vals  = diag(tensorsvd(C)[2])
    normalize!(vals,1)

    #compare
    push!(diffs, sum(abs, oldsvdvals - vals))
    oldsvdvals[:] = vals
    n_it[] += 1

    return state, state
end

function fonedirxmove!((C1,C2), (T1,T2,T4), A, χ)
    #=
        [C1] - [T1] - [C2]
          |      |      |
        [T4] - [ A] - [T2]
          |      |      |
        [C1]*- [T1]*- [C2]*
    =#
    @tensor begin
        C1p[1,2,3]     := C1[1,-1]   * T1[-1,2,3]
        T4p[1,2,3,4,5] := T4[1,-1,5] * A[2,3,4,-1]
        C2p[1,2,3]     := T1[1,2,-1]  * C2[-1,3]
        T2p[1,2,3,4,5] := A[4,-1,2,3] * T2[1,-1,5]
    end

    #renormalize
    @tensor begin
        C1pC2p[1,2,3,4] := C1p[1,2,-1] * C2p[-1,3,4]
        C1pC2pC2pC1p[1,2,3,4] := C1pC2p[1,2,-1,-2] * C1pC2p'[3,4,-1,-2]
        C2pC1pC1pC2p[1,2,3,4] := C1pC2p[-1,-2,1,2] * C1pC2p'[-1,-2,3,4]
    end

    U, = tensoreig(C1pC2pC2pC1p, ((1,2),(3,4)), truncfun = svdtrunc_maxχ(χ))
    V, = tensoreig(C1pC2pC2pC1p, ((1,2),(3,4)), truncfun = svdtrunc_maxχ(χ))'

    @tensor begin
        C1[1,2]   = U'[-1,-2,1] * C1p[-1,-2,2]
        T4[1,2,3] = U'[-1,-2,1] * T4p[-1,-2,2,-4,-3] * U[-3,-4,3]
        C2[1,2]   = C2p[1,-1,-2] * V[-2,-1,2]
        T2[1,2,3] = V'[-1,-2,1] * T2p[-1,-2,2,-4,-3] * V[-3,-4,3]
    end

    mval = maximum(diag(tensorsvd(C1)[2]))
    apply!(C1, x -> x .= x ./ mval)
    mval = maximum(diag(tensorsvd(C2)[2]))
    apply!(C2, x -> x .= x ./ mval)
    apply!(T4, x -> x .= x ./ sqrt(norm(T4)))
    apply!(T2, x -> x .= x ./ sqrt(norm(T2)))
end

function fonedirymove!((C1,C2), (T1,T2,T4), A, χ)
    #=
        [C1] - [T1] - [C2]
          |      |      |
        [T4] - [ A] - [T2]
          |      |      |
        [C1]*- [T1]*- [C2]*
    =#
    @tensor begin
        C1p[1,2,3]     := T4[1,2,-1]  * C1[-1,3]
        T1p[1,2,3,4,5] := A[3,4,-1,2] * T1[1,-1,5]
        C2p[1,2,3]     := C2[1,-1]    * T2[-1,2,3]
    end

    #renormalize
    @tensor begin
        C2pC2pC1pC1p[3,4,1,2] := C1p[1,2,-1] * C1p'[-2,-3,-1] * C2p'[-2,-3,-4] * C2p[3,4,-1]
    end

    P, = tensoreig(C1pC1pC2pC2p, ((1,2),(3,4)), truncfun = svdtrunc_maxχ(χ))'
    Pi = pinv(P)

    @tensor begin
        C1[1,2]   = C1p[1,-1,-2] * P[-2,-1,2]
        T1[1,2,3] = Pi[-1,-2,1] * T1p[-1,-2,2,-4,-3] * P[-3,-4,3]
        C2[1,2]   = C2p[-1,-2,2] * Pi[-1,-2,1]
        #FUUUUU
    end

    mval = maximum(diag(tensorsvd(C1)[2]))
    apply!(C1, x -> x .= x ./ mval)
    mval = maximum(diag(tensorsvd(C2)[2]))
    apply!(C2, x -> x .= x ./ mval)
    apply!(T1, x -> x .= x ./ sqrt(norm(T1)))
end


text:
using Test, SeisModels

@testset "Conversion" begin
    @testset "To LinearLayeredModel" begin
        # PREMPolyModel
        let m = LinearLayeredModel(PREM)
            @test m isa LinearLayeredModel
            # Default spacing
            @test maximum(diff(m.r)) == 20
            @test minimum(diff(m.r)) == 0
            # Specified spacing
            @test maximum(diff(LinearLayeredModel(IASP91, 30).r)) == 30
            @test vp(m, 1000) ≈ vp(PREM, 1000) rtol=0.0001
            @test vs(m, 3600) ≈ vs(PREM, 3600) rtol=0.0001
        end
        # SteppedLayeredModel
        let m = MOON_WEBER_2011, m′ = LinearLayeredModel(m)
            @test m′ isa LinearLayeredModel
            @test vp(m, 100) ≈ vp(m′, 100) rtol=0.0001
        end
    end

    @testset "To PREMPolyModel" begin
        # SteppedLayeredModel
        @test_throws ArgumentError PREMPolyModel(MOON_WEBER_2011, -1)
        @test reffrequency(PREMPolyModel(MOON_WEBER_2011, fref=10.0)) == 10.0
        let m = MOON_WEBER_2011, m′ = PREMPolyModel(MOON_WEBER_2011)
            for property in (:vp, :vs, :density)
                @test all(r -> evaluate(m, property, r) ≈ evaluate(m′, property, r),
                    range(0, stop=surface_radius(m), length=10_000))
            end
            @test mass(m) ≈ mass(m′)
            @test pressure(m, 0) ≈ pressure(m′, 0)
        end
        # LinearLayeredModel
        @test_throws ArgumentError PREMPolyModel(AK135, 0)
        @test_throws ArgumentError PREMPolyModel(AK135, -1)
        @test reffrequency(PREMPolyModel(AK135, fref=0.1)) == 0.1
        let m = AK135, m′ = PREMPolyModel(AK135)
            for property in (:vp, :vs, :density)
                @test all(r -> evaluate(m, property, r) ≈ evaluate(m′, property, r),
                    range(0, stop=surface_radius(m), length=10_000))
            end
            @test mass(m) ≈ mass(m′)
            @test pressure(m, 0) ≈ pressure(m′, 0)
        end
    end

    @testset "To SteppedLayeredModel" begin
        @test_throws ArgumentError SteppedLayeredModel(IASP91, 0)
        for m in (PREM, AK135)
            m′ = SteppedLayeredModel(m)
            @test m′ isa SteppedLayeredModel
            radii = range(0, stop=surface_radius(m), length=10_000)
            for property in (:vp, :vs, :density)
                @test all(r -> ≈(evaluate.((m, m′), property, r)..., rtol=0.01),
                    radii)
            end
            @test mass(m) ≈ mass(m′) rtol=0.001
            @test pressure(m, 0) ≈ pressure(m′, 0) rtol=0.001
        end
        @test maximum(diff(SteppedLayeredModel(PREM, 20).r)) == 20
    end

    @testset "Roundtrip" begin
        let rtol = 0.01
            for T in (PREMPolyModel, LinearLayeredModel, SteppedLayeredModel)
                for m in (PREM, AK135, MOON_WEBER_2011)
                    m isa T && continue # Skip identity conversions
                    m′ = T(m)
                    m″ = typeof(m)(m′)
                    models = (m, m″)
                    radii = range(0, stop=surface_radius(m), length=10_000)
                    for property in (:vp, :vs, :density)
                        @test all(r -> ≈(evaluate.(models, property, r)..., rtol=rtol),
                            radii)
                    end
                    @test ≈(mass.(models)..., rtol=rtol)
                    @test ≈(gravity.(models, surface_radius(m))..., rtol=rtol)
                    @test ≈(pressure.(models, 0)..., rtol=rtol)
                end
            end
        end
    end
end


text:
# # Подготовка данных для тестирования алгоритма выбора начального приближения

# Здесь мы готовим данные, для тестирования алгоритмов.
# Подразумевается, что скрипт будет исполнятся очень редко. Все свои результаты
# он запишет в соответствующие файлы в текущей папке.

using NonLinearReactionAdvectionDiffusionWithFrontData;
using DelimitedFiles;
using Plots; pyplot();

# Выберем параметры таким образом, чтобы решение точно было хорошим.
# Здесь мы генерируем экспериментальные данные, поэтому решим прямую задачу
# всего раз для каждого случая.

x_tp    = 0.05;         # Стартовое местоположение фронта
T_end   = 0.58         # Регулируем конечное местоположение фронта
Nx      = 15000;         # Число интервалов по ``X``
Mt      = 30000;        # Число интервалов по ``T``
α       = 0.005;        # Параметр регуляризации в нахождении начального приближения

# Генерируем данные для набора малых параметров

plts = Any[];           # Вектор для хранения графиков Plots

for (eps, x_tp) in [(0.03, 0.04), (0.01, 0.03), (0.001,0.02)]
    @info "Генерируем данные для ε = $(eps), x_tp = $(x_tp)"
    ε       = eps;          # Крутизна фронта
#+
    a, b, t₀, T, N, M, ε, Xₙ, Tₘ, qₙ, ulₘ, urₘ, u₀ = dparams(x_tp = x_tp,
                                                            ε = ε,
                                                            Nx = Nx,
                                                            Mt = Mt,
                                                            T_end = T_end,
                                                            qfunc = x -> sin(3pi*x)
                                                            );
    u, XX, TP = solve(u₀, Xₙ, N, Tₘ, M, ε, ulₘ, urₘ, qₙ);
    ϕl, ϕr, ϕ, f1_data, f2_data = generate_obs_data(u, Xₙ, N, Tₘ, M, qₙ, ulₘ, urₘ);
#+
    # Убеждаемся в корректности прямого решения на глаз и пишем данные в txt
    directP = draft(u, Xₙ, N, Tₘ, M, title = "Эскиз прямого решения, \\varepsilon = $(ε)")
    savefig("direct_eps_$(ε)_sin.png");
    writedlm("tm_f1_m_eps_$(ε)_sin.txt",[Tₘ f1_data])
    writedlm("xn_qn_eps_$(ε).txt",[Xₙ qₙ])
#+
    # Находим начальные приближения
    #  q_guess = initial_guess(f1_data, Xₙ, N, Tₘ, M, ulₘ, urₘ, α);
    #  p = plot(Xₙ, qₙ, label="Истинное",  title = "Начальное приближение, \\varepsilon = $(ε)");
    #  p = plot!(Xₙ, q_guess, label="Найденное")
    #  savefig("guess_eps_$(ε)_sin.png");
    #  push!(plts, p);
#+
    # Два гаусиана разной высоты, которые частично перекрываются
    gauss_init(x) = 2 * exp(-((x-5.0/13.0)^2)/(1.5/13.0^2)) +
        1.2 * exp((-(x-8.0/13.0)^2)/(3.0/13.0^2));
    qₙ = [gauss_init(x) for x in Xₙ]
    plot(Xₙ, qₙ)
#+
    u, XX, TP = solve(u₀, Xₙ, N, Tₘ, M, ε, ulₘ, urₘ, qₙ);
    ϕl, ϕr, ϕ, f1_data, f2_data = generate_obs_data(u, Xₙ, N, Tₘ, M, qₙ, ulₘ, urₘ);
#+
    # Убеждаемся в корректности прямого решения на глаз и пишем данные в text
    directP = draft(u, Xₙ, N, Tₘ, M, title = "Эскиз прямого решения, \\varepsilon = $(ε)")
    savefig("direct_eps_$(ε)_gauss.png");
    writedlm("tm_f1_m_eps_$(ε)_gauss.txt",[Tₘ f1_data])
    writedlm("xn_qn_eps_$(ε)_gauss.txt",[Xₙ qₙ])
#+
    # Находим начальные приближения
    #  q_guess = initial_guess(f1_data, Xₙ, N, Tₘ, M, ulₘ, urₘ, α);
    #  p = plot(Xₙ, qₙ, label="Истинное",  title = "Начальное приближение, \\varepsilon = $(ε)");
    #  p = plot!(Xₙ, q_guess, label="Найденное")
    #  savefig("guess_eps_$(ε)_gauss.png");
    #  push!(plts, p);

end

# График
plts[1]

# График
plts[2]

# График
plts[3]

# График
plts[4]

# График
plts[5]

# График
plts[6]


text:
function extract_NN(FILE_PATH, OUTPUT_PATH, type)
    @info "Opening file"
    if type == "NDE"
        train_files, train_parameters, losses, loss_scalings, uw_NN, vw_NN, wT_NN, optimizer = jldopen(FILE_PATH, "r") do file
            train_files = file["training_info/train_files"]
            N_stages = length(keys(file["training_data/neural_network/uw"]))
            @show N_stages
            N_data = length(keys(file["training_data/neural_network/uw/$N_stages"]))

            if "parameters" in keys(file["training_info"])
                train_parameters = file["training_info/parameters"]
            else
                train_parameters = nothing
            end

            total_losses = Array{Float32}(undef, N_data)

            if "loss_scalings" in keys(file["training_info"])
                loss_scalings = file["training_info/loss_scalings"]
                profile_losses = similar(total_losses)
                gradient_losses = similar(total_losses)
                u_losses = similar(total_losses)
                v_losses = similar(total_losses)
                T_losses = similar(total_losses)
                ∂u∂z_losses = similar(total_losses)
                ∂v∂z_losses = similar(total_losses)
                ∂T∂z_losses = similar(total_losses)

                @info "Loading Losses"

                for i in 1:length(total_losses)
                    total_losses[i] = file["training_data/loss/total/$N_stages/$i"]
                    profile_losses[i] = file["training_data/loss/profile/$N_stages/$i"]
                    gradient_losses[i] = file["training_data/loss/gradient/$N_stages/$i"]
                    u_losses[i] = file["training_data/loss/u/$N_stages/$i"]
                    v_losses[i] = file["training_data/loss/v/$N_stages/$i"]
                    T_losses[i] = file["training_data/loss/T/$N_stages/$i"]
                    ∂u∂z_losses[i] = file["training_data/loss/∂u∂z/$N_stages/$i"]
                    ∂v∂z_losses[i] = file["training_data/loss/∂v∂z/$N_stages/$i"]
                    ∂T∂z_losses[i] = file["training_data/loss/∂T∂z/$N_stages/$i"]
                end

            else
                loss_scalings = nothing
                @info "Loading Total Loss"
                for i in 1:length(total_losses)
                    total_losses[i] = file["training_data/loss/$N_stages/$i"]
                    profile_losses = nothing
                    gradient_losses = nothing
                    u_losses = nothing
                    v_losses = nothing
                    T_losses = nothing
                    ∂u∂z_losses = nothing
                    ∂v∂z_losses = nothing
                    ∂T∂z_losses = nothing
                end
            end

            losses = (
                total = total_losses,
              profile = profile_losses,
             gradient = gradient_losses,
                    u = u_losses,
                    v = v_losses,
                    T = T_losses,
                 ∂u∂z = ∂u∂z_losses,
                 ∂v∂z = ∂v∂z_losses,
                 ∂T∂z = ∂T∂z_losses
             )

            @info "Loading NN"
            NN_index = argmin(total_losses)
            uw_NN = file["training_data/neural_network/uw/$(N_stages)/$NN_index"]
            vw_NN = file["training_data/neural_network/vw/$(N_stages)/$NN_index"]
            wT_NN = file["training_data/neural_network/wT/$(N_stages)/$NN_index"]

            if "optimizer" in keys(file["training_data"])
                @info "Loading Optimizer State"
                η = file["training_data/optimizer/η/$(N_stages)/$NN_index"]
                β = file["training_data/optimizer/β/$(N_stages)/$NN_index"]
                state = file["training_data/optimizer/state/$(N_stages)/$NN_index"]
                optimizer = (; η, β, state)
            else
                optimizer = nothing
            end
            return train_files, train_parameters, losses, loss_scalings, uw_NN, vw_NN, wT_NN, optimizer
        end
    else
        train_files, train_parameters, losses, NN = jldopen(FILE_PATH, "r") do file
            train_files = file["training_info/train_files"]
            train_parameters = file["training_info/parameters"]
            N_data = length(keys(file["training_data/loss"]))
            losses = zeros(N_data)

            @info "Loading Loss"
            for i in 1:length(losses)
                losses[i] = file["training_data/loss/$i"]
            end

            @info "Loading NN"
            NN_index = argmin(losses)
            NN = file["training_data/neural_network/$NN_index"]
            return train_files, train_parameters, losses, NN
        end
    end

    @info "Writing file"
    if type == "NDE"
        jldopen(OUTPUT_PATH, "w") do file
            @info "Writing Training Info"
            file["training_info/train_files"] = train_files
            file["training_info/parameters"] = train_parameters
            file["training_info/loss_scalings"] = loss_scalings

            file["losses/total"] = losses.total
            file["losses/profile"] = losses.profile
            file["losses/gradient"] = losses.gradient
            file["losses/u"] = losses.u
            file["losses/v"] = losses.v
            file["losses/T"] = losses.T
            file["losses/∂u∂z"] = losses.∂u∂z
            file["losses/∂v∂z"] = losses.∂v∂z
            file["losses/∂T∂z"] = losses.∂T∂z

            file["neural_network/uw"] = uw_NN
            file["neural_network/vw"] = vw_NN
            file["neural_network/wT"] = wT_NN

            file["optimizer/η"] = optimizer.η
            file["optimizer/β"] = optimizer.β
            file["optimizer/state"] = optimizer.state
        end
    else
        jldopen(OUTPUT_PATH, "w") do file
            @info "Writing Training Info"
            file["training_info/train_files"] = train_files
            file["training_info/parameters"] = train_parameters

            @info "Writing Loss"
            file["losses"] = losses

            @info "Writing NN"
            file["neural_network"] = NN
        end
    end

    @info "End"
end

function extract_parameters_modified_pacanowski_philander_optimisation(FILE_PATH, OUTPUT_PATH)
    @info "Opening file"
    train_files, train_parameters, losses, loss_scalings, parameters = jldopen(FILE_PATH, "r") do file
        train_files = file["training_info/train_files"]

        train_parameters = Dict()
        train_parameters["train_epochs"] = file["training_info/train_epochs"]
        train_parameters["train_tranges"] = file["training_info/train_tranges"]
        train_parameters["parameters"] = file["training_info/parameters"]

        if haskey(file, "training_info/loss_scalings")
            loss_scalings = file["training_info/loss_scalings"]
        end

        N_data = length(keys(file["training_data/loss/total"]))

        total_losses = Array{Float32}(undef, N_data)

        u_losses = similar(total_losses)
        v_losses = similar(total_losses)
        T_losses = similar(total_losses)

        ∂u∂z_losses = similar(total_losses)
        ∂v∂z_losses = similar(total_losses)
        ∂T∂z_losses = similar(total_losses)

        @info "Loading Loss"
        for i in 1:length(total_losses)
            total_losses[i] = file["training_data/loss/total/$i"]

            u_losses[i] = file["training_data/loss/u/$i"]
            v_losses[i] = file["training_data/loss/v/$i"]
            T_losses[i] = file["training_data/loss/T/$i"]

            ∂u∂z_losses[i] = file["training_data/loss/∂u∂z/$i"]
            ∂v∂z_losses[i] = file["training_data/loss/∂v∂z/$i"]
            ∂T∂z_losses[i] = file["training_data/loss/∂T∂z/$i"]
        end

        profile_losses = u_losses .+ v_losses .+ T_losses
        gradient_losses = ∂u∂z_losses .+ ∂v∂z_losses .+ ∂T∂z_losses

        losses = (
            total = total_losses,
            profile = profile_losses,
            gradient = gradient_losses,
            u = u_losses,
            v = v_losses,
            T = T_losses,
            ∂u∂z = ∂u∂z_losses,
            ∂v∂z = ∂v∂z_losses,
            ∂T∂z = ∂T∂z_losses,
        )

        @info "Loading Best Parameters"
        NN_index = argmin(total_losses)
        parameters = file["training_data/parameters/$NN_index"]
        return train_files, train_parameters, losses, loss_scalings, parameters
    end

    @info "Writing file"
    jldopen(OUTPUT_PATH, "w") do file
        @info "Writing Training Info"
        file["training_info/train_files"] = train_files
        file["training_info/parameters"] = train_parameters
        file["training_info/loss_scalings"] = loss_scalings

        @info "Writing Loss"

        file["losses/total"] = losses.total
        file["losses/profile"] = losses.profile
        file["losses/gradient"] = losses.gradient

        file["losses/u"] = losses.u
        file["losses/v"] = losses.v
        file["losses/T"] = losses.T

        file["losses/∂u∂z"] = losses.∂u∂z
        file["losses/∂v∂z"] = losses.∂v∂z
        file["losses/∂T∂z"] = losses.∂T∂z

        @info "Writing Best Parameters"
        file["parameters"] = parameters
    end

    @info "End"
end