text:
## 3D Panels
#==========================================================================================#

abstract type AbstractPanel3D <: AbstractPanel end

"""
    Panel3D(p1, p2, p3, p4)

A composite type consisting of 4 Cartesian coordinates `p1, p2, p3, p4` representing corners of a panel in 3 dimensions. The following commutative diagram depicts the order:

```
z → y
↓
x
        p1 —→— p4
        |       |
        ↓       ↓
        |       |
        p2 —→— p3
```
"""
struct Panel3D{T <: Real} <: AbstractPanel3D
    p1 :: SVector{3,T}
    p2 :: SVector{3,T}
    p3 :: SVector{3,T}
    p4 :: SVector{3,T}
end

Panel3D(p1 :: FieldVector{3,T}, p2 :: FieldVector{3,T}, p3 :: FieldVector{3,T}, p4 :: FieldVector{3,T}) where T <: Real = Panel3D{T}(p1, p2, p3, p4)

average_chord(panel :: Panel3D) = (p2(panel) - p1(panel) + p3(panel) - p4(panel)) / 2
average_width(panel :: Panel3D) = (p4(panel) - p1(panel) + p3(panel) - p2(panel)) / 2

"""
    panel_area(panel :: Panel3D)

Compute the area of a Panel3D.
"""
panel_area(panel :: Panel3D) = (norm ∘ cross)(average_chord(panel), average_width(panel))

wetted_area(panels :: Matrix{<: Panel3D}) = sum(panel -> panel_area(panel), panels)

"""
    panel_coords(panel :: Panel3D)

Compute the coordinates of a Panel3D.
"""
panel_coords(panel :: Panel3D) = structtolist(panel)

"""
    transform(panel :: Panel3D, rotation, translation)

Perform an affine transformation on the coordinates of a Panel3D given a rotation matrix and translation vector.
"""
function transform(panel :: Panel3D, rotation, translation)
    p1, p2, p3, p4 = (Translation(translation) ∘ LinearMap(rotation)).(panel_coords(panel))
    Panel3D(p1, p2, p3, p4)
end

"""
    midpoint(panel :: Panel3D)

Compute the midpoint of a Panel3D.
"""
midpoint(panel :: Panel3D) = (p1(panel) + p2(panel) + p3(panel) + p4(panel)) / 4

"""
    panel_normal(panel :: Panel3D)

Compute the normal vector of a Panel3D.
"""
panel_normal(panel :: Panel3D) = let p31 = p3(panel) - p1(panel), p42 = p4(panel) - p2(panel); p31 × p42 end

"""
    transform_normal(panel :: Panel3D, h_l, g_l)

Compute the normal ``n_l``, the normal ``n_0`` of a `Panel3D` perturbed by the control gain ``\\delta_l`` about the hinge axis ``h_l``.
"""
transform_normal(panel :: Panel3D, h_l, g_l) = g_l * cross(h_l, panel_normal(panel))

max_forks_repo_name: GodotMisogi/AeroMDAO


text:
using Plots
using LinearAlgebra
using DelimitedFiles

M_spotify_alpha = readdlm("./outputs/optimization/spotify_estimator_alpha.txt")

features_spotify = readdlm("./features/txt/features.txt")

println(size(M_spotify_alpha))

println(size(features_spotify[1, :]))

x_axis = 1
y_axis = 18

# Procesamiento del Alpha

a = 1.5141161212919718e-6
b = 1.7208871806760954e-7

f = features_spotify[y_axis, :]

x = range(-1, stop = 1, length = 100)
y = range(-1, stop = 1, length = 100)

z = ((a^2) * (x.^2)+(b^2) * (y.^2)).^(1/2)

plotd = plot(x, z)

ylims!((-10, 100))
scatter!(f, [i for i=1:100])

savefig(plotd, string("./outputs/cone_plot_alpha_", x_axis, "_", y_axis, ".png"))

ext: jl

max_stars_repo_name: RodrigoNazar/SuccessPredictor

max_issues_repo_name: RodrigoNazar/SuccessPredictor

max_forks_repo_name: RodrigoNazar/SuccessPredictor


text:
using TestFunctions
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

using NBInclude
@nbinclude("../notebooks/Gallery.ipynb")


text:
using AMQPClient
const VIRTUALHOST = "/"
const HOST = "127.0.0.1"

function send()
    # 1. Create a connection to the localhost or 127.0.0.1 of virtualhost '/'
    connection(; virtualhost=VIRTUALHOST, host=HOST) do conn
        # 2. Create a channel to send messages
        channel(conn, AMQPClient.UNUSED_CHANNEL, true) do chan
            # 3. Declare exchange
            exchange = "direct_logs"
            exchange_declare(chan, exchange, EXCHANGE_TYPE_DIRECT)
            # 4. Get severity and message
            if length(Base.ARGS) >= 3
                severity = Base.ARGS[1]
                received = join(Base.ARGS[2:end], ' ')
            else
                severity = "info"
                received = "Hello World"
            end
            data = convert(Vector{UInt8}, codeunits(received))
            msg = Message(data, content_type="text/plain", delivery_mode=PERSISTENT)
            # 5. Publish the message
            basic_publish(chan, msg; exchange=exchange, routing_key=severity)
            println("Message sent: $received, Severity: $severity")
        end
    end
end

send()

max_forks_repo_name: tomashq/rabbitmq-tutorials


text:
"""
    BurrExpert(k, c, λ)

PDF:

```math
f(x; k, c, \\lambda) = \\frac{kc}{\\lambda} \\left( \\frac{x}{\\lambda} \\right)^{c-1} \\left( 1+ \\left( \\frac{x}{\\lambda} \\right)^{c} \\right)^{-k-1},
\\quad x \\geq 0
```

See also: [Burr Distribution](https://www.mathworks.com/help/stats/burr-type-xii-distribution.html) (Mathworks, implemented in this package),
    [Burr Distribution](https://en.wikipedia.org/wiki/Burr_distribution) (Wikipedia, with λ = 1)

"""
struct BurrExpert{T<:Real} <: NonZIContinuousExpert
    k::T
    c::T
    λ::T
    BurrExpert{T}(k::T, c::T, λ::T) where {T<:Real} = new{T}(k, c, λ)
end

function BurrExpert(k::T, c::T, λ::T; check_args=true) where {T <: Real}
    check_args && @check_args(BurrExpert, k >= zero(k) && c >= zero(c) && λ > zero(λ))
    return BurrExpert{T}(k, c, λ)
end

## Outer constructors
BurrExpert(k::Real, c::Real, λ::Real) = BurrExpert(promote(k, c, λ)...)
BurrExpert(k::Integer, c::Integer, λ::Integer) = BurrExpert(float(k), float(c), float(λ))
BurrExpert() = BurrExpert(1.0, 1.0, 1.0)

## Conversion
function convert(::Type{BurrExpert{T}}, k::S, c::S, λ::S) where {T <: Real, S <: Real}
    BurrExpert(T(k), T(c), T(λ))
end
function convert(::Type{BurrExpert{T}}, d::BurrExpert{S}) where {T <: Real, S <: Real}
    BurrExpert(T(d.k), T(d.c), T(d.λ), check_args=false)
end
copy(d::BurrExpert) = BurrExpert(d.k, d.c, d.λ, check_args=false)

## Loglikelihood of Expert
logpdf(d::BurrExpert, x...) = Distributions.logpdf.(LRMoE.Burr(d.k, d.c, d.λ), x...)
pdf(d::BurrExpert, x...) = Distributions.pdf.(LRMoE.Burr(d.k, d.c, d.λ), x...)
logcdf(d::BurrExpert, x...) = Distributions.logcdf.(LRMoE.Burr(d.k, d.c, d.λ), x...)
cdf(d::BurrExpert, x...) = Distributions.cdf.(LRMoE.Burr(d.k, d.c, d.λ), x...)

## expert_ll, etc
expert_ll_exact(d::BurrExpert, x::Real) = LRMoE.logpdf(d, x)
function expert_ll(d::BurrExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    expert_ll = (yl == yu) ? logpdf.(d, yl) : logcdf.(d, yu) + log1mexp.(logcdf.(d, yl) - logcdf.(d, yu))
    expert_ll = (tu == 0.) ? -Inf : expert_ll
    return expert_ll
end
function expert_tn(d::BurrExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    expert_tn = (tl == tu) ? logpdf.(d, tl) : logcdf.(d, tu) + log1mexp.(logcdf.(d, tl) - logcdf.(d, tu))
    expert_tn = (tu == 0.) ? -Inf : expert_tn
    return expert_tn
end
function expert_tn_bar(d::BurrExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    expert_tn_bar = (tl == tu) ? 0.0 : log1mexp.(logcdf.(d, tu) + log1mexp.(logcdf.(d, tl) - logcdf.(d, tu)))
    return expert_tn_bar
end

exposurize_expert(d::BurrExpert; exposure = 1) = d

## Parameters
params(d::BurrExpert) = (d.k, d.c, d.λ)
function params_init(y, d::BurrExpert)
    pos_idx = (y .> 0.0)

    function init_obj(logparams, y)
        n = length(y)
        c_tmp = exp(logparams[1])
        λ_tmp = exp(logparams[2])

        k_tmp = n / (sum(log1p.((y ./ λ_tmp).^c_tmp)))
        return -1 * ( n*log(c_tmp*k_tmp) - n*(c_tmp-1)*log(λ_tmp) + c_tmp*sum(log.(y)) - (k_tmp+1)*sum(log1p.((y ./ λ_tmp).^c_tmp))
                      + log(c_tmp) - c_tmp/100 + log(λ_tmp) - λ_tmp/100)  # to avoid spurious models
    end

    logparams_init = Optim.minimizer( Optim.optimize(x -> init_obj(x, y[pos_idx]),  [0.0, 0.0] ))

    c_init = exp(logparams_init[1])
    λ_init = exp(logparams_init[2])
    k_init = length(y[pos_idx]) / (sum(log1p.((y[pos_idx] ./ λ_init).^c_init)))

    try
        return BurrExpert(k_init, c_init, λ_init)
    catch;
        BurrExpert()
    end
end

## KS stats for parameter initialization
function ks_distance(y, d::BurrExpert)
    p_zero = sum(y .== 0.0) / sum(y .>= 0.0)
    return max(abs(p_zero-0.0), (1-0.0)*HypothesisTests.ksstats(y[y .> 0.0], LRMoE.Burr(d.k, d.c, d.λ))[2])
end

## Simululation
sim_expert(d::BurrExpert) = Distributions.rand(LRMoE.Burr(d.k, d.c, d.λ), 1)[1]

## penalty
penalty_init(d::BurrExpert) = [2.0 10.0 2.0 10.0 2.0 10.0]
no_penalty_init(d::BurrExpert) = [1.0 Inf 1.0 Inf 1.0 Inf]
penalize(d::BurrExpert, p) = (p[1]-1)*log(d.k) - d.k/p[2] + (p[3]-1)*log(d.c) - d.c/p[4] + (p[5]-1)*log(d.λ) - d.λ/p[6]

## statistics
mean(d::BurrExpert) = mean(LRMoE.Burr(d.k, d.c, d.λ))
var(d::BurrExpert) = var(LRMoE.Burr(d.k, d.c, d.λ))
quantile(d::BurrExpert, p) = quantile(LRMoE.Burr(d.k, d.c, d.λ), p)
function lev(d::BurrExpert, u)
    uu = 1/(1+(u/d.λ)^d.c)
    return d.λ*gamma(float(1+1/d.c))*gamma(float(d.k-1/d.c)) / gamma(float(d.k)) * beta_inc(1+1/d.c, d.k-1/d.c, 1-uu)[1] + u*uu^d.k
end
excess(d::BurrExpert, u) = mean(d) - lev(d, u)

## Misc functions for E-Step
function _int_logy_func(d::BurrExpert, x)
    return (iszero(x) || isinf(x)) ? 0.0 : log(x) * pdf.(d, x)
end

function _int_obs_logY_raw(d::BurrExpert, yl, yu)
    if yl == yu
        return log(yl) * pdf.(d, yl)
    else
        return quadgk.(x -> _int_logy_func(d, x), yl, yu, rtol = 1e-8)[1]
    end
end

function _int_lat_logY_raw(d::BurrExpert, tl, tu)
    return (tl==0 ? 0.0 : quadgk.(x -> _int_logy_func(d, x), 0.0, tl, rtol = 1e-8)[1]) + (isinf(tu) ? 0.0 : quadgk.(x -> _int_logy_func(d, x), tu, Inf, rtol = 1e-8)[1])
end

function _int_lpow_func(c_new, λ_new, d::BurrExpert, x)
    return (iszero(x) || isinf(x)) ? 0.0 : log1p((x/λ_new)^c_new) * pdf.(d, x)
end

function _int_obs_lpow_raw(c_new, λ_new, d::BurrExpert, yl, yu)
    if yl == yu
        return  log1p((yl/λ_new)^c_new) * pdf.(d, yl) # log(1+(yl/λ_new)^c_new) * pdf.(d, yl)
    else
        return quadgk.(x -> _int_lpow_func(c_new, λ_new, d, x), yl, yu, rtol = 1e-8)[1]
    end
end

function _int_lat_lpow_raw(c_new, λ_new, d::BurrExpert, tl, tu)
    return (tl==0 ? 0.0 : quadgk.(x -> _int_lpow_func(c_new, λ_new, d, x), 0.0, tl, rtol = 1e-8)[1]) + (isinf(tu) ? 0.0 : quadgk.(x -> _int_lpow_func(c_new, λ_new, d, x), tu, Inf, rtol = 1e-8)[1])
end

function _burr_lpow_to_k(sum_term_zkz, sum_term_zkzlpow;
    penalty = true, pen_pararms_jk = [1.0 Inf 1.0 Inf 1.0 Inf])
    if penalty
        return (sum_term_zkz + (pen_pararms_jk[1]-1))/(sum_term_zkzlpow + 1/pen_pararms_jk[2])
    else
        return (sum_term_zkz / sum_term_zkzlpow)
    end
end

function _burr_optim_params(lognew,
                        d_old,
                        tl, yl, yu, tu,
                        expert_ll_pos, expert_tn_bar_pos,
                        z_e_obs, z_e_lat, k_e,
                        sum_term_zkzlogy;
                        penalty = true, pen_pararms_jk = [1.0 Inf 1.0 Inf])
    # Optimization in log scale for unconstrained computation
    c_tmp = exp(lognew[1])
    λ_tmp = exp(lognew[2])

    yl_yu_unique = unique_bounds(yl, yu)
    int_obs_lpow_tmp = _int_obs_lpow_raw.(c_tmp, λ_tmp, d_old, yl_yu_unique[:,1], yl_yu_unique[:,2])
    lpow_e_obs = exp.(-expert_ll_pos) .* int_obs_lpow_tmp[match_unique_bounds(hcat(vec(yl), vec(yu)), yl_yu_unique)]
    nan2num(lpow_e_obs, 0.0) # get rid of NaN

    tl_tu_unique = unique_bounds(tl, tu)
    int_lat_lpow_tmp = _int_lat_lpow_raw.(c_tmp, λ_tmp, d_old, tl_tu_unique[:,1], tl_tu_unique[:,2])
    lpow_e_lat = exp.(-expert_tn_bar_pos) .* int_lat_lpow_tmp[match_unique_bounds(hcat(vec(tl), vec(tu)), tl_tu_unique)]
    nan2num(lpow_e_lat, 0.0) # get rid of NaN

    term_zkz = z_e_obs .+ (z_e_lat .* k_e)
    term_zkzlpow = (z_e_obs .* lpow_e_obs) .+ (z_e_lat .* k_e .* lpow_e_lat)

    sum_term_zkz = sum(term_zkz)[1]
    sum_term_zkzlpow = sum(term_zkzlpow)[1]

    k_tmp = _burr_lpow_to_k(sum_term_zkz, sum_term_zkzlpow, penalty = penalty, pen_pararms_jk = pen_pararms_jk)

    obj = sum_term_zkz*(log(k_tmp)+log(c_tmp)-c_tmp*log(λ_tmp)) + sum_term_zkzlogy*c_tmp - (k_tmp+1)*sum_term_zkzlpow
    p = penalty ? (pen_pararms_jk[1]-1)*log(k_tmp) - k_tmp/pen_pararms_jk[2] + (pen_pararms_jk[3]-1)*log(c_tmp) - c_tmp/pen_pararms_jk[4] + (pen_pararms_jk[5]-1)*log(λ_tmp) - λ_tmp/pen_pararms_jk[6] : 0.0
    return (obj + p) * (-1.0)
end

## EM: M-Step
function EM_M_expert(d::BurrExpert,
                     tl, yl, yu, tu,
                     exposure,
                     z_e_obs, z_e_lat, k_e;
                     penalty = true, pen_pararms_jk = [1.0 Inf 1.0 Inf 1.0 Inf])

    expert_ll_pos = expert_ll.(d, tl, yl, yu, tu)
    expert_tn_bar_pos = expert_tn_bar.(d, tl, yl, yu, tu)

    # Further E-Step
    yl_yu_unique = unique_bounds(yl, yu)
    int_obs_logY_tmp = _int_obs_logY_raw.(d, yl_yu_unique[:,1], yl_yu_unique[:,2])
    logY_e_obs = exp.(-expert_ll_pos) .* int_obs_logY_tmp[match_unique_bounds(hcat(vec(yl), vec(yu)), yl_yu_unique)]
    nan2num(logY_e_obs, 0.0) # get rid of NaN

    tl_tu_unique = unique_bounds(tl, tu)
    int_lat_logY_tmp = _int_lat_logY_raw.(d, tl_tu_unique[:,1], tl_tu_unique[:,2])
    logY_e_lat = exp.(-expert_tn_bar_pos) .* int_lat_logY_tmp[match_unique_bounds(hcat(vec(tl), vec(tu)), tl_tu_unique)]
    nan2num(logY_e_lat, 0.0) # get rid of NaN

    # Update parameters
    pos_idx = (yu .!= 0.0)
    term_zkz = z_e_obs[pos_idx] .+ (z_e_lat[pos_idx] .* k_e[pos_idx])
    term_zkz_logY = (z_e_obs[pos_idx] .* logY_e_obs[pos_idx]) .+ (z_e_lat[pos_idx] .* k_e[pos_idx] .* logY_e_lat[pos_idx])

    sum_term_zkzlogy = sum(term_zkz_logY)[1]

    logparams_new = Optim.minimizer( Optim.optimize(x -> _burr_optim_params(x, d,
                                                tl[pos_idx], yl[pos_idx], yu[pos_idx], tu[pos_idx],
                                                expert_ll_pos[pos_idx], expert_tn_bar_pos[pos_idx],
                                                z_e_obs[pos_idx], z_e_lat[pos_idx], k_e[pos_idx],
                                                sum_term_zkzlogy,
                                                penalty = penalty, pen_pararms_jk = pen_pararms_jk),
                                                [log(d.c), log(d.λ)] ))

    c_new = exp(logparams_new[1])
    λ_new = exp(logparams_new[2])

    # Find k_new
    # yl_yu_unique = unique_bounds(yl, yu)
    int_obs_lpow_tmp = _int_obs_lpow_raw.(c_new, λ_new, d, yl_yu_unique[:,1], yl_yu_unique[:,2])
    lpow_e_obs = exp.(-expert_ll_pos) .* int_obs_lpow_tmp[match_unique_bounds(hcat(vec(yl), vec(yu)), yl_yu_unique)]
    nan2num(lpow_e_obs, 0.0) # get rid of NaN

    # tl_tu_unique = unique_bounds(tl, tu)
    int_lat_lpow_tmp = _int_lat_lpow_raw.(c_new, λ_new, d, tl_tu_unique[:,1], tl_tu_unique[:,2])
    lpow_e_lat = exp.(-expert_tn_bar_pos) .* int_lat_lpow_tmp[match_unique_bounds(hcat(vec(tl), vec(tu)), tl_tu_unique)]
    nan2num(lpow_e_lat, 0.0) # get rid of NaN

    term_zkz = z_e_obs[pos_idx] .+ (z_e_lat[pos_idx] .* k_e[pos_idx])
    term_zkzlpow = (z_e_obs[pos_idx] .* lpow_e_obs[pos_idx]) .+ (z_e_lat[pos_idx] .* k_e[pos_idx] .* lpow_e_lat[pos_idx])

    sum_term_zkz = sum(term_zkz)[1]
    sum_term_zkzlpow = sum(term_zkzlpow)[1]

    k_new = _burr_lpow_to_k(sum_term_zkz, sum_term_zkzlpow, penalty = penalty, pen_pararms_jk = pen_pararms_jk)

    return BurrExpert(k_new, c_new, λ_new)
end

## EM: M-Step, exact observations
function _burr_lpow_to_k_exact(sum_term_zkz, sum_term_zkzlpow;
                                penalty = true, pen_pararms_jk = [1.0 Inf 1.0 Inf 1.0 Inf])
    if penalty
        return (sum_term_zkz + (pen_pararms_jk[1]-1))/(sum_term_zkzlpow + 1/pen_pararms_jk[2])
    else
        return (sum_term_zkz / sum_term_zkzlpow)
    end
end

function _burr_optim_params_exact(lognew,
                        ye,
                        z_e_obs,
                        sum_term_zkzlogy;
                        penalty = true, pen_pararms_jk = [1.0 Inf 1.0 Inf 1.0 Inf])
    # Optimization in log scale for unconstrained computation
    c_tmp = exp(lognew[1])
    λ_tmp = exp(lognew[2])

    lpow_e_obs = log1p.((ye./λ_tmp).^c_tmp)
    # nan2num(lpow_e_obs, 0.0) # get rid of NaN

    term_zkz = z_e_obs
    term_zkzlpow = z_e_obs .* lpow_e_obs

    sum_term_zkz = sum(term_zkz)[1]
    sum_term_zkzlpow = sum(term_zkzlpow)[1]

    k_tmp = _burr_lpow_to_k_exact(sum_term_zkz, sum_term_zkzlpow, penalty = penalty, pen_pararms_jk = pen_pararms_jk)

    obj = sum_term_zkz*(log(k_tmp)+log(c_tmp)-c_tmp*log(λ_tmp)) + sum_term_zkzlogy*c_tmp - (k_tmp+1)*sum_term_zkzlpow
    p = penalty ? (pen_pararms_jk[1]-1)*log(k_tmp) - k_tmp/pen_pararms_jk[2] + (pen_pararms_jk[3]-1)*log(c_tmp) - c_tmp/pen_pararms_jk[4] + (pen_pararms_jk[5]-1)*log(λ_tmp) - λ_tmp/pen_pararms_jk[6] : 0.0
    return (obj + p) * (-1.0)
end

function EM_M_expert_exact(d::BurrExpert,
                            ye, exposure,
                            z_e_obs;
                            penalty = true, pen_pararms_jk = [Inf 1.0 Inf])

    # Further E-Step
    logY_e_obs = log.(ye)

    # Update parameters
    pos_idx = (ye .!= 0.0)
    term_zkz = z_e_obs[pos_idx]
    term_zkz_logY = z_e_obs[pos_idx] .* logY_e_obs[pos_idx]

    sum_term_zkzlogy = sum(term_zkz_logY)[1]

    logparams_new = Optim.minimizer( Optim.optimize(x -> _burr_optim_params_exact(x,
                                                ye[pos_idx],
                                                z_e_obs[pos_idx],
                                                sum_term_zkzlogy,
                                                penalty = penalty, pen_pararms_jk = pen_pararms_jk),
                                                [log(d.c), log(d.λ)] ))

    c_new = exp(logparams_new[1])
    λ_new = exp(logparams_new[2])

    # Find k_new
    lpow_e_obs = log1p.((ye ./ λ_new).^c_new)
    nan2num(lpow_e_obs, 0.0)

    term_zkz = z_e_obs[pos_idx]
    term_zkzlpow = z_e_obs[pos_idx] .* lpow_e_obs[pos_idx]

    sum_term_zkz = sum(term_zkz)[1]
    sum_term_zkzlpow = sum(term_zkzlpow)[1]

    k_new = _burr_lpow_to_k_exact(sum_term_zkz, sum_term_zkzlpow, penalty = penalty, pen_pararms_jk = pen_pararms_jk)

    return BurrExpert(k_new, c_new, λ_new)
end

max_stars_repo_name: sparktseung/JuliaLRMoE

max_issues_repo_name: sparktseung/JuliaLRMoE

max_forks_repo_name: sparktseung/JuliaLRMoE


text:
using Hyper
using Test

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


text:
using TriplePendulums
using Test

@testset "TriplePendulums.jl" begin
    @test 1+1 == 2
end


text:
function contour_estimate_eig(A::AbstractMatrix, contour::Contour, B=I;
                samples::Integer=min(100, size(A,1)), ϵ=1e-12, debug=false,
                mixed_prec=false, factorizer=lu, left_divider=ldiv!)

    N, m₀ = size(A,1), samples
    Ctype = if mixed_prec ComplexF32 else ComplexF64 end
    # X = if typeof(A) <: AbstractSparseMatrix
    #     sprandn(Ctype, N, m₀, min(1.0, m₀^2/N))
    # else
    #     randn(Ctype, N, m₀)
    # end
    ### Using sparse X does not work, possibly due to not being i.i.d. mean zero RV
    X = randn(Ctype, N, m₀)

    temp = zeros(Ctype, N, m₀)
    P = zeros(ComplexF64, N, m₀)
    ZmA = similar(A, Ctype)
    nodes = size(contour.nodes, 1)
    est = 0.0

    for i=1:nodes
        ZmA .= Ctype.(B*contour.nodes[i] - A)
		F = factorizer(ZmA)
		est += estimate_trace!(P, temp, X, F, B, left_divider=left_divider) * contour.weights[i]
        if debug print(".") end
    end
    if debug println() end

    return real(est)
end

# function contour_estimate_eig(A::AbstractMatrix, contour::Contour, B=I;
#                 samples::Integer=min(100, size(A,1)), ϵ=1e-12, debug=false,
#                 mixed_prec=false, factorizer=lu, left_divider=ldiv!)
#
#
#     N, m₀ = size(A,1), samples
#     Ctype = if mixed_prec ComplexF32 else ComplexF64 end
#     # X = if typeof(A) <: AbstractSparseMatrix
#     #     sprandn(Ctype, N, m₀, min(1.0, m₀^2/N))
#     # else
#     #     randn(Ctype, N, m₀)
#     # end
#     ### Using sparse X does not work, possibly due to not being i.i.d. mean zero RV
#     X = randn(Ctype, N, m₀)
#
#     temp = zeros(Ctype, N, m₀)
#     P = zeros(ComplexF64, m₀, m₀)
#     ZmA = similar(A, Ctype)
#     nodes = size(contour.nodes, 1)
#     est = 0.0
#
#     for i=1:nodes
#         ZmA .= Ctype.(B*contour.nodes[i] - A)
#         linsolve!(temp, ZmA, X, factorizer, left_divider)
#         mul!(P, X', temp)
#         est += tr(P)*contour.weights[i]/samples
#         if debug print(".") end
#     end
#     if debug println() end
#
#     return real(est)
# end

estimate_trace(A, N) = estimate_trace(convert(typeof(A), sparse(I, size(A)...)), A, N)

estimate_trace(A::AbstractMatrix, B::AbstractMatrix, N) = begin
	@assert size(A,1) == size(A,2)
	@assert size(B,1) == size(B,2)
	M = size(A,1)
	T = promote_type(eltype(A), eltype(B))
	y = zeros(T,M,N)
	x, b = copy(y), copy(y)
	res = estimate_trace!(b, x, y, lu(A), B)
	return res
end

function estimate_trace!(b::AbstractMatrix, x::AbstractMatrix, y::AbstractMatrix, Af::F, B::AbstractMatrix; left_divider::F1=ldiv!, randfill::Bool=true) where {F, F1}
	# Estimate trace of inv(A)*B where Af is a factorization of type F
	@assert size(B,1) == size(B,2)
	M = size(B,2)
	@assert size(y,1) == M
	@assert size(x) == size(y) == size(b)
	N = size(y,2)
	T = eltype(y)
	# Fill sample vectors with iid random numbers
	if randfill @inbounds for i = 1:length(y) y[i] = randn(T) end end
	# Compute b = B*y
	mul!(b, B, y)
	# Compute x = inv(A)*B*y
	left_divider(x, Af, b)
	res = zero(T)
	# Sum over sample vectors
	@inbounds for j = 1:N
		# Compute the dot-product
		@simd for i = 1:M
			res += adjoint(y[i,j]) * x[i,j]
		end
	end
	return res / N
end


text:
export solve_modelARCDiag_HO

function solve_modelARCDiag_HO(nlp_stop, PData :: PDataFact, α:: T; ho_correction :: Symbol = :Shamanskii, λfact = 1.0) where T
    # Solve the ARC subproblem once diagonalized into Δ
    # printstyled("On est dans solve_modelARCDiag ⇊ \n", color = :red)
    nlp_at_x = nlp_stop.current_state
    # Setup the problem
    # M = ones(PData.Δ)
    # M = fill(1, size(PData.Δ))
    M = ones(T, size(PData.Δ))
    # printstyled("On a M ↢ \n", color = :red)

    # ϵ = 1.0e-10 * (1.0 + PData.λ)
    ϵ1 = sqrt(eps(T)) / 100.0
    # printstyled("On a ϵ1 ↢  \n", color = :red)
    ϵ = T.(ϵ1 * (1.0 + PData.λ))
    # printstyled("On a ϵ ↢  \n", color = :red)
    λ = max(ϵ,PData.λ+ϵ)
    # printstyled("On a λ ↢  \n", color = :red)

    d̃ = -(PData.Δ+λ*M) .\ PData.g̃
    # printstyled("On a d̃ ↢  \n", color = :red)
    normd̃ = sqrt(d̃⋅d̃)
    # printstyled("On a normd̃ ↢  \n", color = :red)
    seuil = normd̃/α
    # printstyled("On a seuil ↢  \n", color = :red)

    # Solve the subproblem (Δ + λ I) d̃ = -g̃ such that λ = ||d̃||/α
    d̃,λ = solve_diag(λ,PData.Δ,PData.g̃,seuil,α,ϵ)
    # printstyled("On a d̃, λ ↢  \n", color = :red)

    PData.λ = λ
    # Transform back d̃ into d

    d = AInv(PData,d̃)

    if PData.λ <= λfact
        # println("ici ! 🍆")
        # @show λfact
        # @show ho_correction
        dHO = eval(ho_correction)(nlp_stop, PData, d, PData.g̃)
        # @show dHO
        # @show (norm(dHO) < 2.0 .* α)
        # @show ((-(nlp_at_x.gx + 0.5 .* nlp_at_x.Hx * dHO)⋅dHO) > 0.0)
        # @show nlp_at_x.gx
        # @show nlp_at_x.Hx
        # @show -(nlp_at_x.gx + 0.5 .* nlp_at_x.Hx * dHO)
        # @show (-(nlp_at_x.gx + 0.5 .* nlp_at_x.Hx * dHO)⋅dHO)

        if (norm(dHO) < 2.0 .* α) && ((-(nlp_at_x.gx + 0.5 * nlp_at_x.Hx * dHO)⋅dHO) > 0.0)
            # printstyled("on prend dHO 🐣\n", color = :green)
            return dHO, dHO, PData.λ
        else
            # printstyled("on prend d 🐲 \n", color = :green)
            return d, dHO, PData.λ
        end
    end

    return d, NaN * rand(length(d)), λ
end


text:
import DelayEmbeddings:
    Dataset

import ..TransferOperators: get_binvisits
import CausalityToolsBase

import StaticArrays:
    SVector,
    MVector

include("average_invariant_measure.jl")

export AverageRectangularInvariantMeasure, averagerectangularinvariantmeasure

"""
    average_measure(induced_measures::Vector{InducedRectangularInvariantMeasure})

Compute the average measure from a set of precomputed induced measures.
Assumes that all induced measures were computed at the same final resolution
`ϵF`.
"""
function average_measure(induced_measures::Vector{InducedRectangularInvariantMeasure{T}}) where T
    n_measures = length(induced_measures)
    n_states = length(induced_measures[1].measure_induced.dist)
    measures = zeros(Float64, n_measures, n_states)

    for i in 1:n_measures
        measures[i, :] = induced_measures[i].measure_induced.dist
    end

    avg_measure = zeros(Float64, n_states)
    for i in 1:n_states
        avg_measure[i] = (1/n_measures) * sum(measures[:, i])
    end

    return InvariantDistribution(avg_measure, findall(avg_measure .> 0))
end

export average_measure

"""
    struct AverageRectangularInvariantMeasure{T} <: AbstractRectangularInvariantMeasure where {T}
        points::AbstractArray{T, 2}
        ϵF::Union{Int, Float64, Vector{Int}, Vector{Float64}}
        ϵs::Union{Int, Float64, Vector{Int}, Vector{Float64}}
        induced_measures::Vector{InducedRectangularInvariantMeasure}
        averaged_measure::RectangularInvariantMeasure
    end

An average rectangular invariant measure induced at resolution
`ϵF` by the invariant measure computed at multiple resolutions `ϵs`.

## Fields

- **`points`**: The points for which to estimate the invariant measure.

- **`ϵs`**: The binning schemes. Must be an iterable of `ϵ`s of the following
    types:
        1. `ϵ::Int` divides each axis into `ϵ` intervals of the same size.
        2. `ϵ::Float` divides each axis into intervals of size `ϵ`.
        3. `ϵ::Vector{Int}` divides the i-th axis into `ϵᵢ` intervals of the same size.
        4. `ϵ::Vector{Float64}` divides the i-th axis into intervals of size `ϵᵢ`.

- **`induced_measures`**: The measures induced by each of the resolutions in `ϵ`.

- **`measure_ϵF`**: The invariant measure computed at resolution `ϵF`.

- **`averaged_measure_ϵF`**: The average rectangular invariant measure at resolution `ϵF`
    induced by the measures computed at resolutions `ϵs[1], ϵs[2], …, ϵs[end]`.
"""
struct AverageRectangularInvariantMeasure{T} <: AbstractRectangularInvariantMeasure where {T}
    points::AbstractArray{T, 2}
    ϵF::RectangularBinning
    ϵs::AbstractArray{RectangularBinning, 1}
    induced_measures::Vector{InducedRectangularInvariantMeasure{T}}
    measure_ϵF::RectangularInvariantMeasure
    measure_ϵF_averaged::RectangularInvariantMeasure
end

"""
    AverageRectangularInvariantMeasure(points, ϵF, ϵⱼ)

Compute the average invariant measure the boxes of a hyperrectangular subdivision
of the hypercube spanned by `points` at resolution
`ϵF` induced by the invariant measure computed for the same points at multiple
resolutions provided by `ϵs`.

This constructor returns an `AverageRectangularInvariantMeasure` instance.

## Arguments

- **`points`**: The points for which to estimate the invariant measure.

- **`ϵF`**: The final resolution at which to induce the measure. Must be a valid
    `RectangularBinning` instances. See docs for `RectangularBinning` for details.

- **`ϵs`**: The multiple resolutions from which the measure at resolution `ϵF`
    is induced. Must be an iterable of valid `RectangularBinning` instances.
    See docs for `RectangularBinning` for details.
"""
function averagerectangularinvariantmeasure(points::AbstractArray{T, 2},
        ϵF::RectangularBinning,
        ϵs::AbstractArray{RectangularBinning, 1};
        estimator::Symbol = :TransferOperatorEstimatorRectangularBinning,
        kwargs...) where {T}

    ptsvec = [points[:, i] for i = 1:size(points, 2)]
    edgelengths_ϵF = CausalityToolsBase.get_edgelengths(ptsvec, ϵF)
    edgelengths_ϵs = [CausalityToolsBase.get_edgelengths(ptsvec, ϵi) for (k, ϵi) in enumerate(ϵs)]

    induced_measures = [inducedrectangularinvariantmeasure(points, edgelengths_ϵF, edgelengths_ϵs) for ϵᵢ in ϵs]
    avg_measure = average_measure(induced_measures)

    # Check that measure sums to 1! With too few input points and too fine partitions,
    # this might not be the case
    if sum(avg_measure.dist[isfinite.(avg_measure.dist)]) < 0.99 # some tolerance allowed
        @warn "The number of input points is too low for these resolutions. The average measure will not sum to 1."
    end

    measure_ϵF::RectangularInvariantMeasure = rectangularinvariantmeasure(points, ϵF)
    avg_measure_ϵF::RectangularInvariantMeasure = RectangularInvariantMeasure(
        measure_ϵF.points,
        measure_ϵF.binning_scheme,
        measure_ϵF.axisminima,
        measure_ϵF.edgelengths,
        hcat(measure_ϵF.encoded_points...,),
        measure_ϵF.visited_bins_coordinates,
        measure_ϵF.binvisits,
        measure_ϵF.transfermatrix,
        avg_measure
    )

    AverageRectangularInvariantMeasure(
        points,
        ϵF,#get_edgelengths(points, ϵF.ϵ),#
        ϵs,#[get_edgelengths(points, ϵs[i]) for i = 1:length(ϵs)],
        induced_measures,
        measure_ϵF,
        avg_measure_ϵF
    )
end

averagerectangularinvariantmeasure(points::AbstractArray{T, 2}, ϵF, ϵs) where T =
    averagerectangularinvariantmeasure(
        points,
        RectangularBinning(ϵF),
        [RectangularBinning(ϵi) for (i, ϵi) in enumerate(ϵs)]
        )

function averagerectangularinvariantmeasure(points::Vector{Vector{T}}, ϵF, ϵs) where {T}
    averagerectangularinvariantmeasure(hcat(points...,), ϵF, ϵs)
end

function averagerectangularinvariantmeasure(points::Vector{SVector{T}}, ϵF, ϵs) where {T}
    averagerectangularinvariantmeasure(Array(hcat(points...,)), ϵF, ϵs)
end

function averagerectangularinvariantmeasure(points::Vector{MVector{T}}, ϵF, ϵs) where {T}
    averagerectangularinvariantmeasure(Array(hcat(points...,)), ϵF, ϵs)
end

function averagerectangularinvariantmeasure(points::Dataset, ϵF, ϵs) where {T}
    averagerectangularinvariantmeasure(transpose(Matrix(points)), ϵF, ϵs)
end

function summarise(invm::AverageRectangularInvariantMeasure)
    ϵF = invm.ϵF
    ϵs = invm.ϵs
    ϵF_str = " partition resulting from binning scheme  ϵF = $ϵF"
    ϵs_str = " partitions resulting from binning schemes ϵs = $ϵs"
    D = size(invm.points, 1)
    npts = size(invm.points, 2)

    measure_type = typeof(invm)

    pts_str = "$npts $D-dimensional points\n"
    return join([measure_type, " from ", pts_str, " AVERAGED at", ϵF_str,
                "\n INDUCED by", ϵs_str])

end

Base.show(io::IO, invm::AverageRectangularInvariantMeasure) = println(io, summarise(invm))


text:
toSMatrix() = SMatrix{0,0,Float64}()
toSMatrix(s) = toSMatrix(tuple(s))
toSMatrix(ss::NTuple{M,NTuple{N,Number}}) where {N,M} = toSMatrix(SVector{N}.(ss))
toSMatrix(ss::NTuple{M,SVector{N}}) where {N,M} = hcat(ss...)

toSMatrix(::Type{T}, ss) where {T<:Number} = _toSMatrix(T, toSMatrix(ss))
_toSMatrix(::Type{T}, s::SMatrix{N,M}) where {N,M,T} = convert(SMatrix{N,M,T}, s)
# Dynamic dispatch
toSMatrix(s::AbstractMatrix) = SMatrix{size(s,1), size(s,2)}(s)
toSMatrix(s::AbstractVector) = toSMatrix(Tuple(s))

toSVector(::Tuple{}) = SVector{0,Float64}()
toSVector(v::SVector) = v
toSVector(v::NTuple{N,Any}) where {N} = SVector(v)
toSVector(x::Number) = SVector{1}(x)
toSVector(::Type{T}, v) where {T} = T.(toSVector(v))
toSVector(::Type{T}, ::Tuple{}) where {T} = SVector{0,T}()
# Dynamic dispatch
toSVector(v::AbstractVector) = SVector(Tuple(v))

unitvector(::Type{SVector{L,T}}, i) where {L,T} = SVector{L,T}(unitvector(NTuple{L,T}, i))
unitvector(::Type{CartesianIndex{L}}, i) where {L} = CartesianIndex(unitvector(NTuple{L,Int}, i))
unitvector(::Type{NTuple{L,T}}, i) where {L,T} =ntuple(j -> j == i ? one(T) : zero(T), Val(L))

ensuretuple(s::Tuple) = s
ensuretuple(s) = (s,)

indstopair(s::Tuple) = Pair(last(s), first(s))

filltuple(x, L) = ntuple(_ -> x, L)
filltuple(x, ::NTuple{N,Any}) where {N} = ntuple(_ -> x, Val(N))

# # toSVector can deal with the L=0 edge case, unlike SVector
# frontsvec(x, ::Val{L}) where {L} = toSVector(ntuple(i->x[i], Val(L)))
# tailtuple(x::SVector{N}, ::Val{L}) where {N,L} = ntuple(i->x[N+i-L], Val(L))

@inline tuplejoin() = ()
@inline tuplejoin(x) = x
@inline tuplejoin(x, y) = (x..., y...)
@inline tuplejoin(x, y, z...) = (x..., tuplejoin(y, z...)...)

tuplesplice(s::NTuple{N,T}, ind, el) where {N,T} = ntuple(i -> i === ind ? T(el) : s[i], Val(N))

shiftleft(s::NTuple{N,Any}, n = 1) where {N} = ntuple(i -> s[mod1(i+n, N)], Val(N))
shiftright(s, n = 1) = shiftleft(s, -n)

tupleproduct(p1, p2) = tupleproduct(ensuretuple(p1), ensuretuple(p2))
tupleproduct(p1::NTuple{M,Any}, p2::NTuple{N,Any}) where {M,N} =
    ntuple(i -> (p1[1+fld(i-1, N)], p2[1+mod(i-1, N)]), Val(M * N))

tupleswapfront(tup::NTuple{L}, (i, j)) where {L} =
    i < j ? swap(swap(tup, i => 1), j => 2) : swap(swap(tup, j => 2), i => 1)

swap(tup::NTuple{L}, (i, i´)) where {L} =
    ntuple(l -> tup[ifelse(l == i´, i, ifelse(l == i, i´, l))], Val(L))

tuplepairs(::Val{V}) where {V} = tuplepairs((), ntuple(identity, Val(V)))
tuplepairs(c::Tuple, ::Tuple{}) = c

function tuplepairs(c::Tuple, r::NTuple{V}) where {V}
    t = Base.tail(r)
    c´ = (c..., tuple.(first(r), t)...)
    return tuplepairs(c´, t)
end

# Base.tail(t) .- first(t) but avoiding rounding errors in difference
tuple_minus_first(t::Tuple{T,Vararg{T,D}}) where {D,T} =
    ntuple(i -> ifelse(t[i+1] ≈ t[1], zero(T), t[i+1] - t[1]), Val(D))

firsttail(t::Tuple) = first(t), Base.tail(t)

frontSVector(s::SVector) = SVector(Base.front(Tuple(s)))

mergetuples(ts...) = keys(merge(tonamedtuple.(ts)...))

tonamedtuple(ts::Val{T}) where {T} = NamedTuple{T}(filltuple(0,T))

function deletemultiple_nocheck(dn::SVector{N}, axes::NTuple{M,Int}) where {N,M}
    ind = first(axes)
    dn´ = deleteat(dn, ind)
    taxes = Base.tail(axes)
    axes´ = taxes .- (taxes .> ind)
    return deletemultiple_nocheck(dn´, axes´)
end
deletemultiple_nocheck(dn::SVector, axes::Tuple{}) = dn

_rdr(r1, r2) = (0.5 * (r1 + r2), r2 - r1)

# zerotuple(::Type{T}, ::Val{L}) where {T,L} = ntuple(_ -> zero(T), Val(L))

function padright!(v::Vector, x, n::Integer)
    n0 = length(v)
    resize!(v, max(n, n0))
    for i in (n0 + 1):n
        v[i] = x
    end
    return v
end

padright(sv::StaticVector{E,T}, x::T, ::Val{E}) where {E,T} = sv
padright(sv::StaticVector{E1,T1}, x::T2, ::Val{E2}) where {E1,T1,E2,T2} =
    (T = promote_type(T1,T2); SVector{E2,T}(padright(Tuple(sv), x, Val(E2))))
padright(sv::StaticVector{E,T}, ::Val{E2}) where {E,T,E2} = padright(sv, zero(T), Val(E2))
padright(sv::StaticVector{E,T}, ::Val{E}) where {E,T} = sv
padright(t::NTuple{N´,Any}, x, ::Val{N}) where {N´,N} = ntuple(i -> i > N´ ? x : t[i], Val(N))
padright(t::NTuple{N´,Any}, ::Val{N}) where {N´,N} = ntuple(i -> i > N´ ? 0 : t[i], Val(N))

padright(v, ::Type{<:Number}) = first(v)
padright(v, ::Type{S}) where {T,S<:SVector{<:Any,T}} = padright(v, zero(T), S)
padright(v, x::T, ::Type{S}) where {E,T,S<:SVector{E,T}} =
    SVector{E,T}(ntuple(i -> i > length(v) ? x : convert(T, v[i]), Val(E)))

# Pad element type to a "larger" type
padtotype(s::SMatrix{E,L}, ::Type{S}) where {E,L,E2,L2,T,S<:SMatrix{E2,L2,T}} =
    S(SMatrix{E2,E}(I) * s * SMatrix{L,L2}(I))
padtotype(s::StaticVector, ::Type{S}) where {N,T,S<:SVector{N,T}} =
    padright(T.(s), Val(N))
padtotype(s::SVector{1}, ::Type{T}) where {T<:Number} = T(first(s))
padtotype(x::Number, ::Type{S}) where {E,L,S<:SMatrix{E,L}} =
    S(x * (SMatrix{E,1}(I) * SMatrix{1,L}(I)))
padtotype(s::Number, ::Type{S}) where {N,T,S<:SVector{N,T}} =
    padright(SA[T(s)], Val(N))
padtotype(x::Number, ::Type{T}) where {T<:Number} = T(x)
padtotype(u::UniformScaling, ::Type{T}) where {T<:Number} = T(u.λ)
padtotype(u::UniformScaling, ::Type{S}) where {S<:SMatrix} = S(u)
padtotype(a::AbstractArray, t::Type{<:SVector}) = padright(a, t)

function padtotype(a::AbstractMatrix, ::Type{S}) where {N,M,T,S<:SMatrix{N,M,T}}
    t = ntuple(Val(N*M)) do i
        n, m = mod1(i, N), fld1(i, N)
        @inbounds n > size(a, 1) || m > size(a, 2) ? zero(T) : T(a[n,m])
    end
    return S(t)
end

## Work around BUG: -SVector{0,Int}() isa SVector{0,Union{}}
negative(s::SVector{L,<:Number}) where {L} = -s
negative(s::SVector{0,<:Number}) = s

display_as_tuple(v, prefix = "") = isempty(v) ? "()" :
    string("(", prefix, join(v, string(", ", prefix)), ")")

displayvectors(mat::SMatrix{E,L,<:AbstractFloat}; kw...) where {E,L} =
    ntuple(l -> round.(Tuple(mat[:,l]); kw...), Val(L))
displayvectors(mat::SMatrix{E,L,<:Integer}; kw...) where {E,L} =
    ntuple(l -> Tuple(mat[:,l]), Val(L))

chop(x::T) where {T} = ifelse(abs2(x) < eps(real(T)), zero(T), x)

# pseudoinverse of supercell s times an integer n, so that it is an integer matrix (for accuracy)
pinvmultiple(s::SMatrix{L,0}) where {L} = (SMatrix{0,0,Int}(), 0)
function pinvmultiple(s::SMatrix{L,L´}) where {L,L´}
    L < L´ && throw(DimensionMismatch("Supercell dimensions $(L´) cannot exceed lattice dimensions $L"))
    qrfact = qr(s)
    # Cannot check det(s) ≈ 0 because s can be non-square
    det(qrfact.R) ≈ 0 && throw(ErrorException("Supercell appears to be singular"))
    pinverse = inv(qrfact.R) * qrfact.Q'
    n = round.(Int, det(s's))
    npinverse = round.(Int, n * pinverse)
    return npinverse, n
end

pinverse(::SMatrix{E,0,T}) where {E,T} = SMatrix{0,E,T}() # BUG: workaround StaticArrays bug SMatrix{E,0,T}()'

function pinverse(m::SMatrix)
    @static if VERSION >= v"1.7.0-beta2"
        qrm = qr(m, NoPivot())
    else
        qrm = qr(m)
    end
    return inv(qrm.R) * qrm.Q'
end

issquare(a::AbstractMatrix) = size(a, 1) == size(a, 2)

# normalize_axis_directions(q::SMatrix{M,N}) where {M,N} = hcat(ntuple(i->q[:,i]*sign(q[i,i]), Val(N))...)
padprojector(::Type{S}, ::NTuple{N,Any}) where {S,N} = padprojector(S, N)
padprojector(::Type{S}, n::Integer) where {M,T,S<:SMatrix{M,M,T}} =
    S(Diagonal(SVector(ntuple(i -> ifelse(i > n, zero(T), one(T)), Val(M)))))

_blockdiag(s1::SMatrix{E1,L1,T1}, s2::SMatrix{E2,L2,T2}) where {E1,L1,T1,E2,L2,T2} = hcat(
    ntuple(j->vcat(s1[:,j], zero(SVector{E2,T2})), Val(L1))...,
    ntuple(j->vcat(zero(SVector{E1,T1}), s2[:,j]), Val(L2))...)

function isgrowing(vs::AbstractVector, i0 = 1)
    i0 > length(vs) && return true
    vprev = vs[i0]
    for i in i0 + 1:length(vs)
        v = vs[i]
        v <= vprev && return false
        vprev = v
    end
    return true
end

function ispositive(ndist)
    result = false
    for i in ndist
        i == 0 || (result = i > 0; break)
    end
    return result
end

chop(x::T, x0 = one(T)) where {T<:Real} = ifelse(abs2(x) < eps(T(x0)), zero(T), x)
chop(x::C, x0 = one(R)) where {R<:Real,C<:Complex{R}} = chop(real(x), x0) + im*chop(imag(x), x0)

# function chop!(A::AbstractArray{T}, atol = default_tol(T)) where {T}
#     for (i, a) in enumerate(A)
#         # if abs(a) < atol
#         #     A[i] = zero(T)
#         if !iszero(a) #&& (abs(real(a)) < atol || abs(imag(a)) < atol)
#             A[i] = chop(a)
#         elseif abs(a) > 1/atol || isnan(a)
#             A[i] = T(Inf)
#         end
#     end
#     return A
# end

default_tol(::Type{T}) where {T} = sqrt(eps(real(float(T))))
default_tol(::T) where {T<:Number} = default_tol(T)

function unique_sorted_approx!(v::AbstractVector{T}) where {T}
    i = 1
    xprev = first(v)
    for j in 2:length(v)
        if v[j] ≈ xprev
            xprev = v[j]
        else
            i += 1
            xprev = v[i] = v[j]
        end
    end
    resize!(v, i)
    return v
end

function normalize_columns!(kmat::AbstractMatrix, norm = 1)
    for col in axes(kmat, 2)
        v = view(kmat, :, col)
        iszero(v) || normalize!(v)
        norm == 1 || (v *= norm)
    end
    return kmat
end

normalize_columns(s::SMatrix{L,0}) where {L} = s
normalize_columns(s::SMatrix{0,L}) where {L} = s
normalize_columns(s) = mapslices(v -> v/norm(v), s, dims = 1)

eltypevec(::AbstractMatrix{T}) where {T<:Number} = T
eltypevec(::AbstractMatrix{S}) where {N,T<:Number,S<:SMatrix{N,N,T}} = SVector{N,T}

tuplesort((a,b)::Tuple{<:Number,<:Number}) = a > b ? (b, a) : (a, b)
tuplesort(t::Tuple) = t
tuplesort(::Missing) = missing

# Gram-Schmidt but with column normalization only when norm^2 >= threshold (otherwise zero!)
function orthonormalize!(m::AbstractMatrix, threshold = 0)
    @inbounds for j in axes(m, 2)
        col = view(m, :, j)
        for j´ in 1:j-1
            col´ = view(m, :, j´)
            norm2´ = dot(col´, col´)
            iszero(norm2´) && continue
            r = dot(col´, col)/norm2´
            col .-= r .* col´
        end
        norm2 = real(dot(col, col))
        factor = ifelse(norm2 < threshold, zero(norm2), 1/sqrt(norm2))
        col .*= factor
    end
    return m
end

# Like copyto! but with potentially different tensor orders (adapted from Base.copyto!, see #33588)
function copyslice!(dest::AbstractArray{T1,N1}, Rdest::CartesianIndices{N1},
                    src::AbstractArray{T2,N2}, Rsrc::CartesianIndices{N2}, by = identity) where {T1,T2,N1,N2}
    isempty(Rdest) && return dest
    if length(Rdest) != length(Rsrc)
        throw(ArgumentError("source and destination must have same length (got $(length(Rsrc)) and $(length(Rdest)))"))
    end
    checkbounds(dest, first(Rdest))
    checkbounds(dest, last(Rdest))
    checkbounds(src, first(Rsrc))
    checkbounds(src, last(Rsrc))
    src′ = Base.unalias(dest, src)
    @inbounds for (Is, Id) in zip(Rsrc, Rdest)
        @inbounds dest[Id] = by(src′[Is])
    end
    return dest
end

function append_slice!(dest::AbstractArray, src::AbstractArray{T,N}, Rsrc::CartesianIndices{N}) where {T,N}
    checkbounds(src, first(Rsrc))
    checkbounds(src, last(Rsrc))
    Rdest = (length(dest) + 1):(length(dest) + length(Rsrc))
    resize!(dest, last(Rdest))
    src′ = Base.unalias(dest, src)
    for (Is, Id) in zip(Rsrc, Rdest)
        @inbounds dest[Id] = src′[Is]
    end
    return dest
end

permutations(ss::NTuple) = permutations!(typeof(ss)[], ss, ())

function permutations!(p, s1, s2)
    for (i, s) in enumerate(s1)
        permutations!(p, delete(s1, i), (s2..., s))
    end
    return p
end

permutations!(p, ::Tuple{}, s2) = push!(p, s2)

delete(t::NTuple{N,Any}, i) where {N} = ntuple(j -> j < i ? t[j] : t[j+1], Val(N-1))

function one!(m::StridedMatrix{T}) where {T}
    @inbounds for c in CartesianIndices(m)
        m[c] = ifelse(c[1] == c[2], one(T), zero(T))
    end
    return m
end

######################################################################
# convert a matrix/number block to a matrix/inlinematrix string
######################################################################

_isreal(x) = all(o -> imag(o) ≈ 0, x)
_isimag(x) = all(o -> real(o) ≈ 0, x)

matrixstring(row, x) = string("Onsite[$row] : ", _matrixstring(x))
matrixstring(row, col, x) = string("Hopping[$row, $col] : ", _matrixstring(x))

matrixstring_inline(row, x, digits) = string("Onsite[$row] : ", _matrixstring_inline(round.(x, sigdigits = digits)))
matrixstring_inline(row, col, x, digits) = string("Hopping[$row, $col] : ", _matrixstring_inline(round.(x, sigdigits = digits)))

_matrixstring(x::Number) = numberstring(x)
_matrixstring_inline(x::Number) = numberstring(x)
function _matrixstring(s::SMatrix)
    ss = repr("text/plain", s)
    pos = findfirst(isequal('\n'), ss)
    return pos === nothing ? ss : ss[pos:end]
end

function _matrixstring_inline(s::SMatrix{N}) where {N}
    stxt = numberstring.(transpose(s))
    stxt´ = vcat(stxt, SMatrix{1,N}(ntuple(_->";", Val(N))))
    return string("[", stxt´[1:end-1]..., "]")
end

numberstring(x) = _isreal(x) ? string(" ", real(x)) : _isimag(x) ? string(" ", imag(x), "im") : string(" ", x)

############################################################################################
######## fast sparse copy #  Revise after #33589 is merged #################################
############################################################################################

# Using broadcast .+= instead allocates unnecesarily
function _plain_muladd!(dst, src, α)
    @boundscheck checkbounds(dst, axes(src)...)
    for i in eachindex(src)
        @inbounds dst[i] += α * src[i]
    end
    return dst
end

# Only needed for dense <- sparse (#33589), copy!(sparse, sparse) is fine in v1.5+
fast_sparse_copy!(dst::AbstractMatrix, src::AbstractSparseMatrix, (rows, cols) = axes(src)) =
    _fast_sparse_copy!(dst, src, rows, cols)

_fast_sparse_copy!(dst, src, rows, cols) =
    copy!(dst, view(src, rows, cols))
_fast_sparse_copy!(dst, src, ::Colon, cols) =
    _fast_sparse_copy!(dst, src, axes(src, 1), cols)
_fast_sparse_copy!(dst, src, rows::AbstractRange, ::Colon) =
    _fast_sparse_copy!(dst, src, rows, axes(src, 2))

function _fast_sparse_copy!(dst, src, rows::AbstractRange, cols)
    @boundscheck checkbounds(src, rows, cols)
    @boundscheck checkbounds(dst, eachindex(rows), eachindex(cols))
    rowoffset = first(rows) - 1
    fill!(dst, zero(eltype(src)))
    for (coldst, col) in enumerate(cols)
        for p in nzrange(src, col)
            row = rowvals(src)[p]
            row in rows || continue
            rowdst = row - rowoffset
            @inbounds dst[rowdst, coldst] = nonzeros(src)[p]
        end
    end
    return dst
end

function fast_sparse_muladd!(dst::AbstractMatrix, src::SparseMatrixCSC, α = I, cols = axes(src, 2))
    @boundscheck checkbounds(src, axes(dst, 1), cols)
    @boundscheck checkbounds(dst, axes(src, 1), eachindex(cols))
    @inbounds for (coldst, col) in enumerate(cols)
        for p in nzrange(src, col)
            @inbounds dst[rowvals(src)[p], coldst] += α * nonzeros(src)[p]
        end
    end
    return dst
end

rclamp(r1::UnitRange, r2::UnitRange) = isempty(r1) ? r1 : clamp(minimum(r1), extrema(r2)...):clamp(maximum(r1), extrema(r2)...)

iclamp(minmax, r::Missing) = minmax
iclamp((x1, x2), (xmin, xmax)) = (max(x1, xmin), min(x2, xmax))

############################################################################################
# QR/SVD utils
############################################################################################
# Sparse QR from SparseSuite is also pivoted
pqr(a::SparseMatrixCSC) = qr(a)
pqr(a::Adjoint{T,<:SparseMatrixCSC}) where {T} = qr(copy(a))
pqr(a::Transpose{T,<:SparseMatrixCSC}) where {T} = qr(copy(a))
pqr(a) = qr!(copy(a), Val(true))

getQ(qr::Factorization, cols = :) = qr.Q * Idense(size(qr, 1), cols)
getQ(qr::SuiteSparse.SPQR.QRSparse, cols = :) =  Isparse(size(qr, 1), :, qr.prow) * sparse(qr.Q * Idense(size(qr, 1), cols))
getQ_dense(qr::SuiteSparse.SPQR.QRSparse, cols = :) =  Isparse(size(qr, 1), :, qr.prow) * (qr.Q * Idense(size(qr, 1), cols))

getQ´(qr::Factorization, cols = :) = qr.Q' * Idense(size(qr, 1), cols)
getQ´(qr::SuiteSparse.SPQR.QRSparse, cols = :) = sparse((qr.Q * Idense(size(qr, 1), cols))') * Isparse(size(qr,1), qr.prow, :)

getRP´(qr::Factorization) = qr.R * qr.P'
getRP´(qr::Factorization, rows) = qr.R[rows, :] * qr.P'
getRP´(qr::SuiteSparse.SPQR.QRSparse) = qr.R * Isparse(size(qr, 2), qr.pcol, :)
getRP´(qr::SuiteSparse.SPQR.QRSparse, rows) = Isparse(size(qr, 1), rows, :) * qr.R * Isparse(size(qr, 2), qr.pcol, :)
getRP´thin(qr::Factorization) = getRP´(qr, 1:nonzero_rows(qr.R))

getPR´(qr::Factorization) = qr.P * qr.R'
getPR´(qr::SuiteSparse.SPQR.QRSparse) = Isparse(size(qr, 2), :, qr.pcol) * qr.R'

Idense(n, ::Colon) = Matrix(I, n, n)

function Idense(n, cols)
    m = zeros(Bool, n, length(cols))
    for (j, col) in enumerate(cols)
        m[col, j] = true
    end
    return m
end

# Equivalent to I(n)[rows, cols], but faster
Isparse(n, rows, cols) = Isparse(inds(rows, n), inds(cols, n))

function Isparse(rows, cols)
    rowval = Int[]
    nzval = Bool[]
    colptr = Vector{Int}(undef, length(cols) + 1)
    colptr[1] = 1
    for (j, col) in enumerate(cols)
        push_rows!(rowval, nzval, rows, col)
        colptr[j+1] = length(rowval) + 1
    end
    return SparseMatrixCSC(length(rows), length(cols), colptr, rowval, nzval)
end

inds(::Colon, n) = 1:n
inds(is, n) = is

function push_rows!(rowval, nzval, rows::AbstractUnitRange, col)
    if col in rows
        push!(rowval, 1 + rows[1 + col - first(rows)] - first(rows))
        push!(nzval, true)
    end
    return nothing
end

function push_rows!(rowval, nzval, rows, col)
    for (j, row) in enumerate(rows)
        if row == col
            push!(rowval, j)
            push!(nzval, true)
        end
    end
    return nothing
end

function nonzero_rows(m::AbstractMatrix{T}, atol = default_tol(T)) where {T}
    nrows = size(m, 1)
    zerorows = 0
    for j in 0:nrows-1
        maximum(abs, view(m, nrows - j, :)) > atol && break
        zerorows += 1
    end
    return nrows - zerorows
end

nonzero_rows(m::AbstractSparseMatrix{T}, atol = default_tol(T)) where {T} =
    isempty(rowvals(m)) ? 0 : maximum(rowvals(m))

function svd_sparse(s::AbstractSparseMatrix{T}, atol = default_tol(T)) where {T}
    iszero(s) && return spzeros(T, size(s, 1), 0), spzeros(T, size(s, 1), 0)
    q = qr(s)
    RP´ = getRP´thin(q)
    r = size(RP´, 1)
    U = getQ(q, 1:r)
    q´ = qr(copy(RP´'))
    RP´´ = getRP´(q´)
    r = size(RP´´, 1)
    V = getQ(q´, 1:r)
    U´, S´, V´ = svd(Matrix(RP´´'))
    US = sparse(U´ .* sqrt.(S´)')
    VS = sparse(V´ .* sqrt.(S´)')
    return U * US, V * VS
end

max_stars_repo_name: pablosanjose/Quantica

max_issues_repo_name: pablosanjose/Quantica

max_forks_repo_name: pablosanjose/Quantica


text:
using Pkg
Pkg.activate(".")
using Revise
using Soss
using NamedTupleTools
# using BenchmarkTools
# using Traceur
using TransformVariables
using Plots
using Lazy
using PositiveFactorizations
using LinearAlgebra
using PDMats

# register_primitive(logpdf)

# Kish's effective sample size,
# $n _ { \mathrm { eff } } = \frac { \left( \sum _ { i = 1 } ^ { n } w _ { i } \right) ^ { 2 } } { \sum _ { i = 1 } ^ { n } w _ { i } ^ { 2 } }$

function n_eff(ℓ)
    logw = ℓ.particles
    exp(2 * logsumexp(logw) - logsumexp(2 .* logw))
end

asmatrix(ps) = Matrix([ps...])'

using Distributions: MvNormalStats
function Base.:+(ss1::MvNormalStats,ss2::MvNormalStats)
    tw = ss1.tw + ss2.tw
    s = ss1.s .+ ss2.s
    m = s .* inv(tw)
    s2 = cholesky(Positive, ss1.s2 + ss2.s2) |> Matrix
    MvNormalStats(s, m, s2, tw)
end

using Distributions: MvNormalStats
function Base.:*(k::Real,ss::MvNormalStats)
    tw = k * ss.tw
    s = k * ss.s
    m = s .* inv(tw)
    s2 = cholesky(Positive, k * ss.s2) |> Matrix
    MvNormalStats(s, m, s2, tw)
end

Base.:*(ss::MvNormalStats, k::Real) = k * ss

function logwavg(a,wa,b,wb)
    w1 = logistic(wb-wa)
    w2 = 1-w1
    w1 > 0.9 && return a
    w2 > 0.9 && return b
    w1 * a + w2 * b
end

function makeplot(ℓ,θ,dat)
    @unpack α,β,σ = θ
    @unpack x,y=dat
    ℓmax = maximum(ℓ.particles)
    xs = range(extrema(x)...,length=100)
    plt = scatter(x,y, legend=false)
    # plot!(plt,xs, α .+ β .* xs, legend=false)
    @inbounds for j in eachindex(ℓ.particles)
        alph = exp(ℓ[j]-ℓmax)
        if alph > 0.5
            plot!(plt,xs, α[j] .+ β[j] .* xs, alpha=exp(ℓ[j]-ℓmax), legend=false)
        end
    end
    plt
end

# Expected value of x log-weighted by ℓ
function expect(x,ℓ)
    w = exp(ℓ - maximum(ℓ))
    mean.(x*w) / mean(w)
end

function self_outer(v)
    reshape(v, :, 1) * reshape(v, 1, :)
end

@inline function visStep(N,logp,q)
    x = Particles(N,q)
    ℓ = logp(x) - logpdf(q,x)
    μ = expect(x,ℓ)
    Σ = cholesky(Positive, expect(0.5 * self_outer(x-μ),ℓ)) |> PDMat
    return (x,ℓ,μ,Σ)
end

function runInference(m; kwargs...)
    ℓp_pre = sourceLogdensity(m) |> eval
    ℓp(θ) = Base.invokelatest(ℓp_pre, θ)
    t = xform(m; kwargs...)
    d = t.dimension
    tnames = propertynames(t.transformations)

    # make kwargs a NamedTuple
    kwargs = (;zip(keys(kwargs), values(kwargs))...)

    function logp(x::Vector{T} where T)
        prep(θ) = merge(kwargs, θ)
        f = transform_logdensity(t, ℓp ∘ prep, x)
    end

    logp(x) = logp([x])

    θ = particles(prior(m))

    x = @as xx θ begin
        mapslices(namedtuple(tnames),asmatrix(xx),dims=1)
        inverse(t).(xx)
        hcat(xx...)
        mapslices(Particles,xx,dims=2)
        vec(xx)
    end
    N = 1000
    q = fit_mle(MvNormal, asmatrix(x))
    x = Particles(N,q)
    ℓ = logp(x) - logpdf(q,x)

    plts = []
    neff = [n_eff(ℓ)]
    numiters = 100
    elbo = [max(mean(ℓ),expect(ℓ,ℓ))]
    for j in 1:numiters
        (x,ℓ,μ,Σ) = visStep(N,logp,q)
        @inbounds @show j, elbo[end], neff[end]
        push!(elbo, max(mean(ℓ),expect(ℓ,ℓ)))
        push!(neff, n_eff(ℓ))
        @inbounds neff[end] > 950 && break

        # η = (neff[end] )/(neff[end] +neff[end-1])
        η = 0.8 # learning rate
        μ = η * μ + (1-η) * q.μ
        Σ = η * Σ + (1-η) * q.Σ
        q = MvNormal(μ,1.5 * Σ)
    end

    x = Particles(N,q)
    θ = t(x)
    ℓ = logp(x) - logpdf(q,x)
    (θ,q,ℓ,elbo,neff)
end

m = @model x begin
        μ ~ Cauchy()
        σ ~ HalfCauchy()
        N = length(x)
        x ~ Normal(μ, σ) |> iid(N)
    end

thedata = (x=randn(1000),)

(θ,q,ℓ,elbo,neff) = runInference(m; thedata...)
plot(elbo, label="ELBO")
plot(neff, label="Effective Sample Size")

θ
# (α = thedata.α, β = thedata.β)#, σ = thedata.σ)

@unpack μ,σ = θ

scatter(μ.particles,σ.particles, alpha=exp(ℓ - maximum(ℓ)).particles)


text:
module TestSettable
using Test
using Setfield
import ConstructionBase
# If no constructor is defined explicitly, don't generate any
# inner-consturctor; let Julia generate the default constructor; i.e.,
# @settable should be a no-op.
@settable struct NoConstructor{A,B}
    b::B
end
ConstructionBase.constructorof(::Type{T}) where {T <: NoConstructor} = T

@testset "NoConstructor" begin
    s1 = NoConstructor{:a,Int}(1)
    s2 = @set s1.b = 2
    @test s2 === NoConstructor{:a,Int}(2)
end

# If the position-only constructor is defined explicitly, don't
# generate any inner-consturctor; i.e., @settable should be a no-op.
@settable struct ExplicitConstructor{A,B}
    b::B
    ExplicitConstructor{A,B}(b::B) where {A,B} = new{A,B}(b)
end
ConstructionBase.constructorof(::Type{T}) where {T <: ExplicitConstructor} = T

@testset "ExplicitConstructor" begin
    s1 = ExplicitConstructor{:a,Int}(1)
    s2 = @set s1.b = 2
    @test s2 === ExplicitConstructor{:a,Int}(2)
end

# If non- position-only constructor is given, generate the
# position-only constructor from it.
@settable struct TypedConstructor{A,B}
    a::A
    b::B
    function TypedConstructor{A,B}(a::A; b::B=0) where {A,B}
        # This assertion has to be executed in the generated inner
        # constructor as well:
        @assert a > b
        return new{A,B}(a, b)
    end
end
ConstructionBase.constructorof(::Type{T}) where {T <: TypedConstructor} = T

@testset "TypedConstructor" begin
    s1 = TypedConstructor{Int,Int}(1)
    s2 = @set s1.b = -2
    @test s2 === TypedConstructor{Int,Int}(1,-2)
    @test_throws AssertionError @set s1.b = 2
end

@settable struct UntypedConstructor{A,B}
    a::A
    b::B
    function UntypedConstructor(a::A; b::B=0) where {A,B}
        # This assertion has to be executed in the generated inner
        # constructor as well:
        @assert a > b
        return new{A,B}(a, b)
    end
end

@testset "UntypedConstructor" begin
    s1 = UntypedConstructor(1, 0)
    s2 = @set s1.b = -2
    @test s2 === UntypedConstructor(1, b=-2)
    @test_throws AssertionError @set s1.b = 2
end

# @settable must (1) choose a correct constructor for generating
# positional-only constructor and (2) do so without overriding
# existing constructor.
# https://github.com/jw3126/Setfield.jl/pull/18#discussion_r172649307
@settable struct ManyConstructors
    a
    b
    ManyConstructors(a) = new(a,1)
    ManyConstructors(;a=1,b=2) = new(a,b)
end

@testset "ManyConstructors" begin
    s1 = ManyConstructors(a=0)
    s2 = @set s1.b = -2
    @test s2 === ManyConstructors(a=0, b=-2)
end

# Mimic @add_kwonly from Reconstructables.jl and DiffEqBase.jl.
macro kwonly()
    esc(quote
        WithKwOnly(a; b=2) = new(a, b)
        WithKwOnly(; a=error("`a` is mandatory"), b=2) = new(a, b)
        end)
end

@settable struct WithKwOnly
    a
    b
    @kwonly
end

@testset "WithKwOnly" begin
    # Keyword-only constructor works:
    @test WithKwOnly(0) === WithKwOnly(a=0)
    @test_throws ErrorException WithKwOnly()

    s1 = WithKwOnly(0)
    s2 = @set s1.b = -2
    @test s2 === WithKwOnly(0, b=-2)
end
end


text:
# From https://github.com/fonsp/Pluto.jl/issues/1237

export latexify_md, wrap_tex

latexify_md(args...; kwargs...) = Markdown.LaTeX(repr(MIME"text/latex"(), latexify(args...; kwargs...)))

wrap_tex(x) = Markdown.LaTeX(repr(MIME"text/latex"(), x))

#=
Issue also includes a suggestion to use
begin
	using Latexify
	using LaTeXStrings

	Markdown.html(io::IO, ls::LaTeXString) =
		Markdown.html(io, Markdown.LaTeX(
			repr(MIME"text/latex"(), ls)
		))

	Markdown.htmlinline(io::IO, ls::LaTeXString) =
		Markdown.htmlinline(io, Markdown.LaTeX(
			repr(MIME"text/latex"(), ls)
		))
end
=#


text:
function next(w, η)
	return (w + η*sin(2π*w) + 1) % 1
end
function dnext(w, η)
	return 1.0 + 2π*η*cos(2π*w)
end
function d2next(w, η)
	return -2π*2π*η*sin(2π*w)
end

max_stars_repo_name: nishaChandramoorthy/neuralOMET

max_issues_repo_name: nishaChandramoorthy/neuralOMET

max_forks_repo_name: nishaChandramoorthy/neuralOMET


text:
import DataFrames

using Dance.Router

function get_df(headers::Dict{String, String}) :: DataFrames.DataFrame
    return DataFrames.DataFrame(A = 1:4, B = ["M", "F", "F", "M"])
end

function get_dict(headers::Dict{String, String}) :: Dict{Symbol, Int64}
    return Dict(:a => 123)
end

function post_df(df::DataFrames.DataFrame, headers::Dict{String, String}) :: Tuple{DataFrames.DataFrame, Dict{String, String}}
    return df, Dict("foo" => "bar")
end

function post_dict(params_dict::Dict{Symbol, Union{Float64, Int64, String}}, dict::Dict, headers::Dict{String, String}) :: Dict{String, Int64}
    for key in keys(dict)
        dict[key] = params_dict[:value]
    end
    return dict
end

route("/dict", get_dict; method=GET)
route("/dataframe", get_df; method=GET)
route(r"/post/dict/(?<value>\d{2})", post_dict)
route("/post/dataframe", post_df)


text:
using IntervalArithmetic, IntervalRootFinding, StaticArrays, ForwardDiff
using Test

function all_unique(rts)
    all(root_status.(rts) .== :unique)
end

function roots_dist(rt1::Root{T}, rt2::Root{T}) where {T <: Interval}
    d = dist(interval(rt1), interval(rt2))
    return sum(d)
end

function roots_dist(rt1::Root{T}, rt2::Root{T}) where {T <: IntervalBox}
    return sum(dist.(interval(rt1), interval(rt2)))
end

function roots_dist(rt1::Root{Complex{T}}, rt2::Root{Complex{T}}) where {T <: Interval}
    dreal = dist(real(interval(rt1)), real(interval(rt2)))
    dimag = dist(imag(interval(rt1)), imag(interval(rt2)))

    return sum(dreal) + sum(dimag)
end

function test_newtonlike(f, deriv, X, method, nsol, tol=1e-10)
    rts = roots(f, X, method)
    @test length(rts) == nsol
    @test all_unique(rts)
    @test sum(roots_dist.(rts, roots(f, deriv, X, method))) < tol
end

newtonlike_methods = [Newton, Krawczyk]

@testset "1D roots" begin
    # Default
    rts = roots(sin, -5..5)
    @test length(rts) == 3
    @test all_unique(rts)

    # Bisection
    rts = roots(sin, -5..6, Bisection, 1e-3)
    @test length(rts) == 3

    # Refinement
    rts = roots(sin, rts, Newton)
    @test all_unique(rts)

    for method in newtonlike_methods
        test_newtonlike(sin, cos, -5..5, method, 3)
    end

    # Infinite interval
    rts = roots(x -> x^2 - 2, -∞..∞)
    @test length(rts) == 2
end

in_solution_set(point, solution_intervals) =
    any(interval -> point in interval, solution_intervals)

@testset "2D roots" begin
    f(x, y) = SVector(x^2 + y^2 - 1, y - 2x)
    f(X) = f(X...)
    X = (-6..6) × (-6..6)

    # Bisection
    rts = roots(f, X, Bisection, 1e-3)
    exact_sol = [sqrt(1/5), 2sqrt(1/5)]
    @test in_solution_set(exact_sol, interval.(rts))
    @test in_solution_set(-exact_sol, interval.(rts))

    for method in newtonlike_methods
        deriv = xx -> ForwardDiff.jacobian(f, xx)
        test_newtonlike(f, deriv, X, method, 2)
    end

    # Infinite interval
    X = IntervalBox(-∞..∞, 2)
    rts = roots(f, X, Newton)
    @test length(rts) == 2
end

# From R docs:
# https://www.rdocumentation.org/packages/pracma/versions/1.9.9/topics/broyden

@testset "3D roots" begin
    function g(x)
        (x1, x2, x3) = x
        SVector(    x1^2 + x2^2 + x3^2 - 1,
                    x1^2 + x3^2 - 0.25,
                    x1^2 + x2^2 - 4x3
                )
    end

    X = (-5..5)
    XX = IntervalBox(X, 3)

    for method in newtonlike_methods
        rts = roots(g, XX, method)
        @test length(rts) == 4
        @test all_unique(rts)
        deriv = xx -> ForwardDiff.jacobian(g, xx)
        @test rts == roots(g, deriv, XX, method)
    end
end

@testset "Out of domain" begin
    for method in newtonlike_methods
        @test length(roots(log, -100..2, method)) == 1
        @test length(roots(log, -100..(-1), method)) == 0
    end
end

@testset "Infinite domain" begin
    for method in newtonlike_methods
        rts = roots(x -> x^2 - 2, -∞..∞, method)
        @test length(filter(isunique, rts)) == 2
    end
end

@testset "NaN return value" begin
    f(xx) = ( (x, y) = xx; SVector(log(y/x) + 3x, y - 2x) )
    X = IntervalBox(-100..100, 2)
    for method in newtonlike_methods
        rts = roots(f, X, method)
        @test length(filter(isunique, rts)) == 1
        @test length(filter(x -> contains_zero(x.interval), rts)) == 1
    end
end

@testset "Stationary points" begin
    f(xx) = ( (x, y) = xx; sin(x) * sin(y) )
    gradf = ∇(f)
    XX = IntervalBox(-5..6, 2)
    tol = 1e-5

    for method in newtonlike_methods
        deriv = xx -> ForwardDiff.jacobian(gradf, xx)
        test_newtonlike(gradf, deriv, XX, method, 25, tol)
    end
end

@testset "Complex roots" begin
    X = -5..5
    Xc = Complex(X, X)
    f(z) = z^3 - 1

    # Default
    rts = roots(f, Xc)
    @test length(rts) == 3

    # Bisection
    rts = roots(f, Xc, Bisection, 1e-3)
    @test length(rts) == 7

    for method in newtonlike_methods
        deriv = z -> 3*z^2
        test_newtonlike(f, deriv, Xc, method, 3, 1e-10)
    end
end


text:
module test_lupin_mergewith

using Test
using Lupin

a = (foo = 0.0, bar = 42.0)
b = (baz = 17, bar = 4711)
@test mergewith(+, a, b) == mergewith(+)(a, b) == (foo = 0.0, bar = 4753.0, baz = 17)

A = [(a=1, b=4), (a=2, b=5), (a=3, b=6)]
@test mergewith(+, A...) == mergewith(+)(A...) == (a=6, b=15)

getfirstarg(x, y) = x
getsecondarg(x, y) = y
@test mergewith(getfirstarg, (a=1,), (a=2,))          == (a = 1,)
@test mergewith(getsecondarg, (a=1,), (a=2,))         == (a = 2,)
@test mergewith(getsecondarg, (a=1,), (a=2,), (a=3,)) == (a = 3,)

@test mergewith(vcat, (a=[1],), (a=[2],)) == (a = [1, 2],)
@test mergewith(vcat, (a=["a"],), (a=[2],)) == (a = ["a", 2],)

end # module test_lupin_mergewith


text:
using ODEInterfaceDiffEq, DiffEqProblemLibrary, DiffEqBase
using Test

using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_mm_linear

prob = prob_ode_mm_linear

@test_throws ErrorException solve(prob,dopri5(),dt=1//2^4)

@test_throws ErrorException solve(prob,dop853();dt=1//2^(4))

@test_throws ErrorException solve(prob,odex();dt=1//2^(4))

sol =solve(prob,seulex();dt=1//2^(4))

sol =solve(prob,radau();dt=1//2^(4))

sol =solve(prob,radau5();dt=1//2^(4))

sol =solve(prob,rodas();dt=1//2^(4))

@test_throws ErrorException solve(prob,ddeabm();dt=1//2^(4))

sol =solve(prob,ddebdf();dt=1//2^(4))


text: Dict{String, Any}("true" => Dict{String, Any}(), "false" => Dict{String, Any}(), "nan" => Dict{String, Any}(), "inf" => Dict{String, Any}())


text:
# Parameters for Liquid Density
DataLiquidDensity = Dict("1,1,1-Trichloroethane"=>[[0.47476,0.27258,0.29333],[545.00,242.75,545.00]],
            "1,1,2-Trichloroethane"=>[[0.47455,0.25475,0.31000],[602.00,236.50,602.00]],
            "1,1-Dichloroethane"=>[[0.41231,0.26533,0.28700],[523.00,176.19,523.00]],
            "1,2-Dichloroethane"=>[[0.46501,0.28742,0.31041],[561.00,237.49,561.00]],
            "1,3-Butadiene"=>[[0.24597,0.27227,0.29074],[425.37,164.25,425.37]],
            "1,4-Dioxane"=>[[0.37018,0.28130,0.30470],[587.00,284.95,587.00]],
            "1-Butanol"=>[[0.26891,0.26674,0.24570],[562.93,183.85,562.93]],
            "1-Butene"=>[[0.23224,0.26630,0.28530],[419.59,87.80,419.59]],
            "1-Decene"=>[[0.23981,0.25776,0.28562],[617.05,206.89,617.05]],
            "1-Nonane"=>[[0.23364,0.25556,0.28571],[595.65,219.63,595.65]],
            "1-Octene"=>[[0.23682,0.25649,0.28571],[566.60,171.45,566.60]],
            "1-Propanol"=>[[0.27684,0.27200,0.24940],[536.71,146.95,536.71]],
            "Acetaldehyde"=>[[0.28207,0.26004,0.27760],[461.00,150.15,461.00]],
            "Acetic acid"=>[[0.35182,0.26954,0.26843],[592.71,289.81,592.71]],
            "Acetic anhydride"=>[[0.33578,0.24080,0.26990],[569.15,200.15,569.15]],
            "Acetone"=>[[0.27728,0.25760,0.29903],[508.20,178.45,508.20]],
            "Acrylic acid"=>[[0.34645,0.25822,0.30701],[615.00,286.65,615.00]],
            "Ammonia"=>[[0.23689,0.25471,0.28870],[405.65,195.41,405.65]],
            "Aniline"=>[[0.31190,0.25000,0.28571],[699.00,267.13,699.00]],
            "Benzene"=>[[0.30090,0.26770,0.28180],[562.16,278.68,562.16]],
            "Butyric acid"=>[[0.31132,0.26192,0.27997],[628.00,267.95,628.00]],
            "Carbon disulfide"=>[[0.47589,0.28749,0.32260],[552.00,161.58,552.00]],
            "Carbon dioxide"=>[[0.46382,0.26160,0.29030],[304.19,216.58,304.19]],
            "Carbon monoxide"=>[[0.29818,0.27655,0.29053],[132.92,68.15,132.92]],
            "Carbon tetrachloride"=>[[0.56607,0.27663,0.29000],[556.35,250.33,556.35]],
            "Chlorine"=>[[0.56600,0.27315,0.28830],[417.15,172.12,417.15]],
            "Chloroform"=>[[0.49807,0.25274,0.28766],[536.40,209.63,536.40,1.480]],
            "Cyclohexane"=>[[0.27376,0.27408,0.28511],[553.54,279.69,553.54]],
            "Cyclohexanol"=>[[0.29681,0.24340,0.28570],[625.15,296.60,625.15]],
            "Cyclopropane"=>[[0.25880,0.27400,0.28571],[397.91,145.73,397.91]],
            "Dichloromethane"=>[[0.45965,0.25678,0.29020],[510.00,178.01,510.00]],
            "Diethyl ether"=>[[0.27267,0.27608,0.29358],[466.70,156.85,466.70]],
            "Diethyl ketone"=>[[0.25635,0.24291,0.27364],[560.95,234.18,560.95]],
            "Dimethylamine"=>[[0.24110,0.26785,0.24800],[437.65,180.96,437.65]],
            "Ethane"=>[[0.20087,0.27330,0.28330],[305.42,90.35,305.42]],
            "Ethanol"=>[[0.26570,0.26395,0.23670],[516.25,159.05,516.25]],
            "Ethyl acetate"=>[[0.30654,0.25856,0.27800],[523.30,189.60,523.30]],
            "Ethyl chloride"=>[[0.32259,0.27464,0.23140],[460.35,136.75,460.35]],
            "Ethylbenzene"=>[[0.28889,0.26438,0.29210],[617.17,178.20,617.17]],
            "Ethylene"=>[[0.21428,0.28061,0.28571],[282.36,104.01,282.36]],
            "Ethylene glycol"=>[[0.32503,0.25499,0.17200],[645.00,260.15,645.00]],
            "Ethylene oxide"=>[[0.31402,0.26089,0.28253],[469.15,161.45,469.15]],
            "Fluorine"=>[[0.57092,0.28518,0.29000],[144.31,53.48,144.31]],
            "Glycerol"=>[[0.34908,0.24902,0.15410],[723.00,291.33,723.00]],
            "Hydrogen"=>[[0.03125,0.34730,0.27560],[33.18,13.95,33.18]],
            "Hydrogen chloride"=>[[0.44134,0.26957,0.31870],[324.65,159.97,324.65]],
            "Hydrogen cyanide"=>[[0.19501,0.18589,0.2820],[456.65,259.91,465.65]],
            "Hydrogen peroxide"=>[[0.43776,0.24982,0.28770],[730.15,272.74,730.15]],
            "i-Butane"=>[[0.22281,0.27294,0.27301],[408.14,113.54,408.14]],
            "Methane"=>[[0.15998,0.28810,0.27700],[190.58,90.67,190.58]],
            "Methanol"=>[[0.27197,0.27192,0.23310],[512.58,175.47,512.58]],
            "Methyl bromide"=>[[0.60859,0.26292,0.28030],[467.00,179.55,467.00]],
            "Methyl chloride"=>[[0.35821,0.26109,0.28690],[416.25,175.45,416.25]],
            "Methyl isobutyl ketone"=>[[0.26654,0.25887,0.28571],[571.40,189.15,571.40]],
            "Methylamine"=>[[0.20168,0.21405,0.22750],[430.05,179.69,430.05]],
            "m-Xylene"=>[[0.27866,0.25925,0.27243],[617.05,225.30,617.05]],
            "Naphthalene"=>[[0.30619,0.25037,0.28300],[748.35,353.43,748.35]],
            "n-Butane"=>[[0.22827,0.27240,0.28630],[425.18,134.86,425.18]],
            "n-Decane"=>[[0.23276,0.25240,0.28570],[618.45,243.49,618.45]],
            "n-Heptane"=>[[0.23237,0.26020,0.27910],[540.26,182.57,540.26]],
            "n-Hexane"=>[[0.23242,0.26500,0.27810],[507.43,177.84,507.43]],
            "Nitrobenzene"=>[[0.36140,0.24731,0.28570],[719.00,278.91,719.00]],
            "Nitrogen"=>[[0.31205,0.28479,0.29250],[126.10,63.15,126.10]],
            "n-Nonane"=>[[0.23364,0.25556,0.28571],[595.65,219.63,595.65]],
            "n-Octane"=>[[0.22807,0.25476,0.26940],[568.83,216.38,568.83]],
            "n-Pentane"=>[[0.23143,0.26923,0.28215],[469.65,143.42,469.65]],
            "Oxygen"=>[[0.43533,0.28772,0.29240],[154.58,54.35,154.58]],
            "o-Xylene"=>[[0.28381,0.26083,0.27410,630.37,247.98,630.37]],
            "Phenol"=>[[0.41476,0.32162,0.32120],[694.25,314.06,694.25]],
            "p-Xylene"=>[[0.27984,0.26003,0.27900],[616.26,286.41,616.26]],
            "Pyridine"=>[[0.30752,0.24333,0.30450],[619.95,231.53,619.95]],
            "Styrene"=>[[0.29383,0.26315,0.28570],[648.00,242.54,648.00]],
            "Toluene"=>[[0.29999,0.27108,0.29889],[591.79,178.18,591.79]],
            "Vinyl acetate"=>[[0.31843,0.25803,0.28270],[524.00,180.35,524.00]],
            "Water"=>[[0.34710,0.27400,0.28571],[647.13,273.16,647.13]])

# Parameters for Gas Viscosity
DataGasViscosity=Dict("1,1,1-Trichloroethane"=>[-19.216 4.0308E-01 -6.8618E-05 347.0 997.0],
            "1,1,2-Trichloroethane"=>[-8.293 3.3989E-01 -5.3678E-05 387.0 987.0],
            "1,1-Dichloroethane"=>[-12.991 4.0085E-01 -1.1779E-04 320.0 472.0],
            "1,2-Dichloroethane"=>[1.025 3.1792E-01 -4.1853E-05 322.0 561.0],
            "1,3-Butadiene"=>[10.256 2.6833E-01 -4.1148E-05 250.0 650.0],
            "1,4-Dioxane"=>[-16701 3.4988E-01 -5.3736E-05 374.0 994.0],
            "1-Butanol"=>[-11.144 2.8790E-01 -5.6275E-05 391.0 1000.0],
            "1-Butene"=>[-9.143 3.1562E-01 -8.4164E-05 175.0 800.0],
            "1-Decene"=>[4.015 1.8083E-01 -3.8216E-05 444.0 1000.0],
            "1-Nonane"=>[-6.802 1.8688E-01 3.4929E-07 273.0 773.0],
            "1-Octene"=>[2.722 2.0327E-01 -4.3879E-05 394.0 994.0],
            "1-Propanol"=>[-14.894 3.2171E-01 -5.8021E-05 200.0 1000.0],
            "Acetaldehyde"=>[0.069 3.0246E-01 -4.2372E-05 294.0 1000.0],
            "Acetic acid"=>[-28.660 2.3510E-01 2.2087E-04 366.0 523.0],
            "Acetic anhydride"=>[-1.485 2.8869E-01 -2.3391E-05 295.0 993.0],
            "Acetone"=>[-4.055 2.6655E-01 -5.6936E-06 300.0 650.0],
            "Acrylic acid"=>[-6.532 3.0600E-01 -4.6620E-05 287.0 1000.0],
            "Ammonia"=>[-7.874 3.6700E-01 -4.4700E-06 195.0 1000.0],
            "Aniline"=>[-6.918 2.5935E-01 -3.4348E-05 458.0 1000.0],
            "Benzene"=>[-0.151 2.5706E-01 -8.9797E-06 287.0 628.0],
            "Butyric acid"=>[-5.781 2.6159E-01 -3.4903E-05 268.0 1000.0],
            "Carbon disulfide"=>[-7.700 3.6594E-01 -2.5416E-05 273.0 583.0],
            "Carbon dioxide"=>[11.811 4.9838E-01 -1.0851E-04 195.0 1500.0],
            "Carbon monoxide"=>[23.811 5.3944E-01 -1.5411E-04 68.0 1250.0],
            "Carbon tetrachloride"=>[-7.745 3.9481E-01 -1.1150E-04 280.0 800.0],
            "Chloroform"=>[-4.392 3.7309E-01 -5.1696E-05 250.0 700.0],
            "Chlorine"=>[-3.571 4.8700E-01 -8.5300E-05 200.0 1000.0],
            "Cyclohexane"=>[1.190 2.4542E-01 -3.8334E-05 315.0 600.0],
            "Cyclohexanol"=>[-13.542 2.9086E-01 -3.9472E-05 400.0 1000.0],
            "Cyclopropane"=>[-9.521 3.7037E-01 -1.3138E-04 240.0 446.0],
            "Dichloromethane"=>[-20.372 4.3745E-01 -7.7549E-05 273.0 993.0],
            "Diethyl ether"=>[-7.932 3.0235E-01 -7.3858E-05 200.0 1000.0],
            "Diethyl ketone"=>[-3.593 2.4308E-01 -3.7613E-05 234.0 1000.0],
            "Dimethylamine"=>[-9.275 2.8958E-01 -1.3875E-06 250.0 450.0],
            "Ethane"=>[0.514 3.3449E-01 -7.1071E-05 150.0 1000.0],
            "Ethanol"=>[1.499 3.0741E-01 -4.4479E-05 200.0 1000.0],
            "Ethyl acetate"=>[-9.259 3.0725E-01 -7.1069E-05 190.0 1000.0],
            "Ethyl chloride"=>[0.458 3.2827E-01 -1.2467E-05 213.0 523.0],
            "Ethylbenzene"=>[-4.267 2.4735E-01 -5.4264E-05 409.0 1000.0],
            "Ethylene"=>[-3.985 3.8726E-01 -1.1227E-04 150.0 1000.0],
            "Ethylene glycol"=>[-7.178 3.1246E-01 -4.4028E-05 260.0 1000.0],
            "Ethylene oxide"=>[-12.180 3.7672E-01 -7.7599E-05 250.0 1000.0],
            "Fluorine"=>[-0.811 8.9800E-01 -3.9600E-04 90.0 500.0],
            "Glycerol"=>[-23.119 2.8879E-01 -3.4277E-05 563.0 993.0],
            "Hydrogen"=>[27.758 2.1200E-01 -3.2800E-05 150.0 1500.0],
            "Hydrogen chloride"=>[-9.118 5.5500E-01 -1.1100E-04 200.0 1000.0],
            "Hydrogen cyanide"=>[-8.486 9.0368E-02 7.9146E-05 300.0 425.0],
            "Hydrogen peroxide"=>[8.039 2.7000E-01 8.2900E-05 373.0 600.0],
            "i-Butane"=>[-4.731 2.9131E-01 -8.0995E-05 150.0 1000.0],
            "Methane"=>[3.844 4.0112E-01 -1.4303E-04 91.0 850.0],
            "Methanol"=>[-14.236 3.8935E-01 -6.2762E-05 240.0 1000.0],
            "Methyl bromide"=>[-27.740 5.5901E-01 -7.0942E-05 260.0 440.0],
            "Methyl chloride"=>[-1.374 3.8627E-01 -4.8650E-05 230.0 700.0],
            "Methyl isobutyl ketone"=>[-3.237 2.3310E-01 -3.5612E-05 189 1000],
            "Methylamine"=>[-5.334 3.4181E-01 -7.4297E-05 267.0 1000.0],
            "m-Xylene"=>[-21.620 2.7820E-01 -6.0531E-05 250.0 1000.0],
            "Naphthalene"=>[-16.789 2.5406E-01 -3.5495E-05 353.0 1000.0],
            "n-Butane"=>[-4.946 2.9001E-01 -6.9665E-05 150.0 1200.0],
            "n-Decane"=>[-7.297 1.8506E-01 -4.8008E-06 287.0 783.0],
            "n-Heptane"=>[-10.378 2.4401E-01 -5.4003E-05 338.0 700.0],
            "n-Hexane"=>[-8.222 2.6229E-01 -5.7366E-05 300.0 1000.0],
            "Nitrobenzene"=>[-16.569 2.9184E-01 -2.5523E-05 450.0 1000.0],
            "Nitrogen"=>[42.606 4.7500E-01 -9.8800E-05 150.0 1500.0],
            "n-Nonane"=>[-6.802 1.8688E-01 3.4929E-07 273.0 773.0],
            "n-Octane"=>[3.940 1.6640E-01 1.4470E-05 374.0 670.0],
            "n-Pentane"=>[-3.202 2.6746E-01 -6.6178E-05 303.0 900.0],
            "Oxygen"=>[44.224 5.6200E-01 -1.1300E-04 150.0 1500.0],
            "o-Xylene"=>[-19.763 2.8022E-01 -5.9293E-05 250.0 1000.0],
            "Phenol"=>[-7.185 2.7179E-01 -3.6205E-05 455.0 1000.0],
            "Propane"=>[-5.462 3.2722E-01 -1.0672E-04 193.0 750.0],
            "p-Xylene"=>[-17.226 2.5098E-01 -2.8232E-05 286.0 1000.0],
            "Pyridine"=>[-5.739 2.7135E-01 -1.7202E-05 369.0 998.0],
            "Styrene"=>[-10.035 2.5191E-01 -3.7932E-05 243.0 1000.0],
            "Sulfur dioxide"=>[-11.103 5.0200E-01 -1.0800E-04 200.0 1000.0],
            "Toluene"=>[1.787 2.3566E-01 -9.3508E-06 275.0 600.0],
            "Vinyl acetate"=>[-7.462 3.0466E-01 -5.7544E-05 346.0 1000.0],
            "Water"=>[-36.826 4.2900E-01 -1.6200E-05 280.0 1073.0])

# Parameters for Liquid Viscosity
DataLiquidViscosity=Dict("1,1,1-Trichloroethane"=>[-3.9096 7.0709E+02 7.5847E-03 -9.1662E-06 243 545],
            "1,1,2-Trichloroethane"=>[-3.2716 6.8810E+02 4.8932E-03 -5.4671E-06 237 602],
            "1,1-Dichloroethane"=>[-3.8388 5.9046E+02 8.0953E-03 -9.9210E-06 176 523],
            "1,2-Dichloroethane"=>[-0.1656 2.7576E+02 -3.3493E-03 1.4093E-06 245 561],
            "1,3-Butadiene"=>[0.3772 7.9658E+01 -5.8889E-03 2.9221E-06 250 425],
            "1,4-Dioxane"=>[-7.5724 1.3813E+03 1.3556E-02 -1.1464E-05 288 587],
            "1-Butanol (n-Butanol)"=>[-5.3970 1.3256E+03 6.2223E-03 -5.5062E-06 250 563],
            "1-Butene"=>[-4.9218 4.9503E+02 1.4390E-02 -2.0853E-05 160 420],
            "1-Decene"=>[-6.8845 1.1003E+03 1.4341E-02 -1.3520E-05 273 617],
            "1-Nonane (n-Nonane)"=>[-6.0742 9.6861E+02 1.2677E-02 -1.2675E-05 220 596],
            "1-Octene"=>[-5.6209 8.1305E+02 1.2523E-02 -1.3384E-05 250 567],
            "1-Propanol (n-Propanol)"=>[-3.7702 9.9151E+02 4.0836E-03 -5.4586E-06 220 537],
            "Acetaldehyde"=>[-6.6171 6.8123E+02 1.9979E-02 -2.5563E-05 260 461],
            "Acetic acid"=>[-3.8937 7.8482E+02 6.6650E-03 -7.5606E-06 290 593],
            "Acetic anhydride"=>[-17.3580 2.3611E+03 4.2734E-02 -3.8202E-05 265 569],
            "Acetone"=>[-7.2126 9.0305E+02 1.8385E-02 -2.0353E-05 223 508],
            "Acrylic acid"=>[5.9215 2.4408E+03 3.4383E-02 -2.7677E-05 293 615],
            "Ammonia"=>[-8.5910 8.7640E+02 2.6810E-02 -3.6120E-05 195 406],
            "Aniline"=>[-13.8625 2.5109E+03 2.5681E-02 -1.8281E-05 268 699],
            "Benzene"=>[-7.4005 1.1815E+03 1.4888E-02 -1.3713E-05 285 562],
            "Butyric acid"=>[-7.9846 1.3636E+03 1.6315E-02 -1.4511E-05 268 628],
            "carbon disulfide"=>[-9.1108 1.1216E+03 2.3216E-02 -2.2648E-05 235 552],
            "Carbon dioxide"=>[-19.4921 1.5948E+03 7.9274E-02 -1.2025E-04 219 304],
            "Carbon monoxide"=>[-1.1224 5.7858E+01 -4.9174E-03 8.2233E-06 69 133],
            "Carbon tetrachloride"=>[-6.4564 1.0379E+03 1.4021E-02 -1.4107E-05 265 556],
            "Chlorine"=>[-0.7681 1.5140E+02 -8.0650E-04 4.0750E-07 172 417],
            "Chloroform"=>[-47831 6.9902E+02 1.0929E-02 -1.2244E-05 210 536],
            "Cyclohexane"=>[.7423 -2.5322E+02 -1.6927E-02 1.2472E-05 285 554],
            "Cyclohexanol"=>[-5.3792 1.8793E+03 1.7011E-03 1.0187E-07 303 625],
            "Cyclopropane"=>[-3.2541 3.2192E+02 9.9766E-03 -1.8191E-05 146 398],
            "Dichloromethane"=>[5.1043 6.8653E+02 1.2459E-02 -1.4540E-05 208 510],
            "Diethyl ether"=>[-8.5060 1.0020E+03 2.2753E-02 -2.5780E-05 233 467],
            "Diethyl ketone"=>[-9.2905 1.2716E+03 2.1925E-02 -2.1036E-05 274 561],
            "Dimethylamine"=>[-11.5558 1.2126E+03 3.4999E-02 -4.1253E-05 240 438],
            "Ethane"=>[-4.2694 2.8954E+02 1.7111E-02 -3.6092E-05 98 305],
            "Ethanol"=>[-6.4406 1.1176E+03 1.3721E-02 -1.5465E-05 240 516],
            "Ethyl acetate"=>[-3.6861 5.5228E+02 8.0018E-03 -1.0439E-05 220 523],
            "Ethyl chloride"=>[-4.4279 5.1891E+02 1.2035E-02 -1.6620E-05 150 460],
            "Ethylbenzene"=>[-5.2585 8.3065E+02 1.0784E-02 -1.0618E-05 210 617],
            "Ethylene"=>[-4.5611 3.0811E+02 1.8030E-02 -3.8145E-05 105 282],
            "Ethylene glycol"=>[-16.9728 3.1886E+03 3.2537E-02 -2.4480E-05 261 645],
            "Ethylene oxide"=>[-5.7794 6.7020E+02 1.5686E-02 -1.9462E-05 190 469],
            "Fluorine"=>[-1.5760 8.5630E+01 -4.0730E-04 -2.7250E-06 54 145],
            "Glycerol"=>[-18.2152 4.2305E+03 2.8705E-02 -1.8648E-05 293 723],
            "Hydrogen"=>[-7.0154 4.0791E+01 2.3714E-01 -4.0830E-03 14 33],
            "Hydrogen chloride"=>[-1.5150 1.9460E+02 3.0670E-03 -1.3760E-05 159 325],
            "Hydrogen cyanide"=>[-12.0812 1.3183E+03 3.5234E-02 -4.0185E-05 260 457],
            "Hydrogen peroxide"=>[-1.6150 5.0380E+02 3.5010E-04 -1.1680E-06 273 728],
            "i-Butane (iso-Butane)"=>[-13.4207 1.3131E+03 4.4329E-02 -5.5793E-05 190 408],
            "Methane"=>[-7.3801 3.1925E+02 4.7934E-02 -1.4120E-04 91 191 - 0.020],
            "Methyl bromide"=>[-9.5533 1.0306E+03 2.8322E-02 -3.1920E-05 193 467],
            "Methyl chloride"=>[-7.3473 8.5395E+02 1.9485E-02 -2.3484E-05 249 416],
            "Methyl isobutyl ketone"=>[-3.0570 5.0050E+02 6.5038E-03 -8.8243E-06 246 571],
            "Methylamine"=>[-9.4670 9.8286E+02 2.8918E-02 -3.5672E-05 180 430],
            "m-Xylene"=>[-6.0517 9.2460E+02 1.2583E-02 -1.1850E-05 225 617],
            "Naphthalene"=>[-10.3716 1.8572E+03 1.9320E-02 -1.4012E-05 353 748],
            "n-Butane"=>[-6.8590 6.7393E+02 2.1973E-02 -3.0686E-05 180 425],
            "n-Decane"=>[-6.0716 1.0177E+03 1.2247E-02 -1.1892E-05 243 618],
            "n-Heptane"=>[-5.7782 8.0587E+02 1.3355E-02 -1.4794E-05 183 540],
            "n-Hexane"=>[-5.0715 6.5536E+02 1.2349E-02 -1.5042E-05 178 507],
            "Nitrobenzene"=>[-7.7710 1.4019E+03	1.4653E-02 -1.1512E-05 273 719],
            "Nitrogen"=>[-15.6104 4.6505E+02 1.6259E-01 -6.3353E-04 63 125],
            "n-Nonane"=>[-6.0742 9.6861E+02 1.2677E-02 -1.2675E-05 220 596],
            "n-Octane"=>[-5.9245 8.8809E+02 1.2955E-02 -1.3596E-05 216 569],
            "n-Pentane"=>[-7.1711 7.4736E+02 2.1697E-02 -2.7176E-05 143	470],
            "Oxygen"=>[-5.0957 1.7983E+02 3.9779E-02 -1.4664E-04 54 150],
            "o-Xylene"=>[-7.8805 1.2500E+03 1.6116E-02 -1.3993E-05 268 630],
            "Phenol"=>[.5349 4.2620E+02 -9.1577E-03 6.2322E-06 318 694],
            "p-Xylene"=>[-9.4655 1.4400E+03 1.9910E-02 -1.6994E-05 288 616],
            "Pyridine"=>[-6.8100 1.1496E+03 1.3229E-02 -1.1661E-05 232 620],
            "Styrene"=>[-8.0291 1.2666E+03 1.6127E-02 -1.3475E-05 243 648],
            "Sulfur dioxide"=>[-2.6700 4.0670E+02 6.1440E-03 -1.2540E-05 200 431],
            "Toluene"=>[-5.1649 8.1068E+02 1.0454E-02 -1.0488E-05 200 592],
            "Vinyl acetate"=>[-9.0671 1.1863E+03 2.2663E-02	-2.3208E-05	250 524],
            "Water"=>[-10.2158 1.7925E+03 1.7730E-02 -1.2631E-05 273 643])


text:
"""
First-order MAML (FOMAML) from
Model-Agnostic Meta-Learning for Fast Adapdation of Deep Networks [1].

[1] https://arxiv.org/abs/1703.03400
"""
function fomaml(model; meta_opt=Descent(0.01), inner_opt=Descent(0.02), epochs=30_000,
              n_tasks=3, train_batch_size=10, eval_batch_size=10, eval_interval=1000)

    weights = params(model)
    dist = Uniform(-5, 5)
    testx = Float32.(range(-5, stop=5, length=50))

    for i in 1:epochs
        prev_weights = deepcopy(Flux.data.(weights))

        for _ in 1:n_tasks
            task = SineWave()

            x = Float32.(rand(dist, train_batch_size))
            y = task(x)
            grad = Flux.Tracker.gradient(() -> Flux.mse(model(x'), y'), weights)

            for w in weights
                w.data .-= Flux.Optimise.apply!(inner_opt, w.data, grad[w].data)
            end

            testy = task(testx)
            grad = Flux.Tracker.gradient(() -> Flux.mse(model(testx'), testy'), weights)

            # reset weights and accumulate gradients
            for (w1, w2) in zip(weights, prev_weights)
                w1.data .= w2
                w1.grad .+= grad[w1].data
            end

        end

        Flux.Optimise._update_params!(meta_opt, weights)

        if i % eval_interval == 0
            @printf("Iteration %d, evaluating model on random task...\n", i)
            evalx = Float32.(rand(dist, eval_batch_size))
            eval_model(model, evalx, testx, SineWave())
        end

    end
end

max_stars_repo_name: JuliaTagBot/MetaLearning

max_issues_repo_name: JuliaTagBot/MetaLearning

max_forks_repo_name: JuliaTagBot/MetaLearning


text:
using Documenter
using Literate
using TightBindingLattice

example_directory = joinpath(@__DIR__, "..", "examples")
output_directory = joinpath(@__DIR__, "src/generated")
for filename in [
    "example_group_isomorphism.jl",
    "example_honeycomb_symmetry.jl",
    "example_orthocube.jl",
    "example_kagome_symmetry.jl",
    "example_little_group_2d.jl",
    "example_little_group_3d.jl",
    "example_little_symmetry_kspace_honeycomb.jl",
    "example_point_group_4mm.jl",
    "example_point_group_elements.jl",
]
    example = joinpath(example_directory, filename)
    Literate.markdown(example, output_directory)
end

makedocs(
    modules=[TightBindingLattice],
    doctest=true,
    sitename="TightBindingLattice.jl",
    format=Documenter.HTML(prettyurls=!("local" in ARGS)),
    authors="Kyungmin Lee",
    checkdocs=:all,
    pages = [
        "Home" => "index.md",
        "lattice.md",
        "Symmetry" => [
            "Symmetry/operation-space.md",
            "Symmetry/symmetry-space.md",
            "Symmetry/irrep.md",
        ],
        "Examples" => [
            "Group Isomorphism" => "generated/example_group_isomorphism.md",
            "OrthoCube" => "generated/example_orthocube.md",
            "Point Group Elements" => "generated/example_point_group_elements.md",
            "Point Group 4mm (C₄ᵥ)" => "generated/example_point_group_4mm.md",
            "Little Group in 2D" => "generated/example_little_group_2d.md",
            "Little Group in 3D" => "generated/example_little_group_3d.md",
            "Honeycomb Lattice" => "generated/example_honeycomb_symmetry.md",
            "Honeycomb Lattice in k-space" => "generated/example_little_symmetry_kspace_honeycomb.md",
            "Kagome Lattice" => "generated/example_kagome_symmetry.md",
        ],
        "API Reference" => [
            "Basic" => "API/basic.md",
            "Group" => "API/group.md",
            "Symmetry Operation" => "API/symmetry-operation.md",
            "Symmetry" => "API/symmetry.md",
            "Symmetry Embedding" => "API/symmetry-embedding.md",
            "Symmetry Irrep Component" => "API/symmetry-irrep-component.md",
        ]
    ]
)

deploydocs(
    repo = "github.com/kyungminlee/TightBindingLattice.jl.git",
    devbranch = "dev",
)


text:
# Ост - "Е6" - восток, Норд - "N8" - север, Запад - "W4" - запад, Юг - "S2" - юг.
function main_function!14(r::Robot)
    for side in (Nord, West, Sud, Ost)
        num_steps = put_marker_in_direction_and_get_num_steps(r,side)
        try_move(r,inverse(side), num_steps)
    end
    putmarker!(r)
end

function put_marker_in_direction_and_get_num_steps(r::Robot,side::HorizonSide)
    num_steps=0
    while try_move(r, side)
        putmarker!(r)
        num_steps += 1
    end
    return num_steps
end

try_move(r::Robot, side::HorizonSide, num_steps::Int) =
for _ in 1:num_steps
    try_move(r,side)
end

function try_move(r::Robot, direct_side::HorizonSide)::Bool
    orthogonal_side = clockwise(direct_side)
    reverse_side = inverse(orthogonal_side)
    num_steps=0
    if !isborder(r,direct_side)
        move!(r,direct_side)
        result=true
    else
        while isborder(r,direct_side)
            if !isborder(r, orthogonal_side)
                move!(r, orthogonal_side)
                num_steps += 1
            else
                break
            end
        end
        if !isborder(r,direct_side)
            move!(r,direct_side)
            while isborder(r,reverse_side)
                move!(r,direct_side)
            end
            result = true
        else
            result = false
        end
        while num_steps>0
            num_steps=num_steps-1
            move!(r,reverse_side)
        end
    end
    return result
end

clockwise(side::HorizonSide) = HorizonSide(mod(Int(side)-1,4))
inverse(side::HorizonSide) = HorizonSide(mod(Int(side)+2, 4))


text:
module GPUShowcases

end


text:
module ControlNeuralODE

using Base: @kwdef, ifelse
using Base.Filesystem: mkpath

using Dates: now
using LazyGrids: ndgrid
using ArgCheck: @argcheck
using Formatting: format, sprintf1
using ProgressMeter: Progress, next!
using Infiltrator: @infiltrate
using Statistics: mean, std
using LineSearches: BackTracking
using InfiniteOpt:
    Infinite,
    InfiniteModel,
    OrthogonalCollocation,
    @infinite_parameter,
    @variables,
    @constraint,
    @constraints,
    @objective,
    set_start_value_function,
    optimizer_with_attributes,
    optimizer_model,
    optimize!,
    solution_summary,
    termination_status,
    supports,
    value,
    ∂
using Ipopt: Ipopt
using Optim: LBFGS, BFGS
using ApproxFun: Chebyshev, Fun, Interval
using ForwardDiff: ForwardDiff
using ReverseDiff: ReverseDiff
using Zygote: Zygote
using Flux: glorot_uniform, ADAM, NADAM, sigmoid_fast, tanh_fast
using SciMLBase:
    ODEProblem,
    DECallback,
    remake,
    AbstractODEAlgorithm,
    AbstractSensitivityAlgorithm,
    AbstractODEProblem
using DiffEqCallbacks: FunctionCallingCallback
using OrdinaryDiffEq: solve, AutoTsit5, Rosenbrock23  # Tsit5
using DiffEqSensitivity:
    InterpolatingAdjoint,
    QuadratureAdjoint,
    ZygoteVJP,
    ReverseDiffVJP,
    ForwardDiffSensitivity
using GalacticOptim: GalacticOptim
# using GalacticOptim: AutoForwardDiff, AutoZygote  # does not precompile
using DiffEqFlux: FastChain, FastDense, initial_params, sciml_train
using UnicodePlots: lineplot, lineplot!
using Serialization: serialize, deserialize
using JSON3: JSON3
using CSV: CSV
using Tables: table

# using PyCall: PyObject
using PyPlot: plt, matplotlib, ColorMap

import CommonSolve: solve

export batch_reactor, van_der_pol, reference_tracking, bioreactor, semibatch_reactor

# https://diffeqflux.sciml.ai/stable/ControllingAdjoints/#Choosing-a-sensealg-in-a-Nutshell
const INTEGRATOR = AutoTsit5(Rosenbrock23())
const SENSEALG = ForwardDiffSensitivity()
# const SENSEALG = QuadratureAdjoint(; autojacvec=ReverseDiffVJP())
# const SENSEALG = QuadratureAdjoint(; autojacvec=ZygoteVJP())
# InterpolatingAdjoint(; autojacvec=ZygoteVJP(), checkpointing=true)

include("auxiliaries.jl")
include("simulators.jl")
include("training.jl")
include("plotting.jl")

# case studies
include("systems/batch_reactor.jl")
include("systems/van_der_pol.jl")
include("systems/reference_tracking.jl")
include("systems/bioreactor.jl")
include("systems/semibatch_reactor.jl")

end # module


text:
#!/usr/bin/env julia
#
# Simple DMRG tutorial.  This code integrates the following concepts:
#  - Infinite system algorithm
#  - Finite system algorithm
#
# Copyright 2013-2015 James R. Garrison and Ryan V. Mishmash.
# Open source under the MIT license.  Source code at
# <https://github.com/simple-dmrg/simple-dmrg/>

using Arpack
using LinearAlgebra
using SparseArrays

# Data structures to represent the block and enlarged block objects.
struct Block
    length::Int
    basis_size::Int
    operator_dict::Dict{Symbol,AbstractMatrix{Float64}}
end

struct EnlargedBlock
    length::Int
    basis_size::Int
    operator_dict::Dict{Symbol,AbstractMatrix{Float64}}
end

# For these objects to be valid, the basis size must match the dimension of
# each operator matrix.
isvalid(block::Union{Block,EnlargedBlock}) =
    all(op -> size(op) == (block.basis_size, block.basis_size), values(block.operator_dict))

# Model-specific code for the Heisenberg XXZ chain
model_d = 2  # single-site basis size

Sz1 = [0.5 0.0; 0.0 -0.5]  # single-site S^z
Sp1 = [0.0 1.0; 0.0 0.0]  # single-site S^+

H1 = [0.0 0.0; 0.0 0.0]  # single-site portion of H is zero

function H2(Sz1, Sp1, Sz2, Sp2)  # two-site part of H
    # Given the operators S^z and S^+ on two sites in different Hilbert spaces
    # (e.g. two blocks), returns a Kronecker product representing the
    # corresponding two-site term in the Hamiltonian that joins the two sites.
    J = 1.0
    Jz = 1.0
    return (J / 2) * (kron(Sp1, Sp2') + kron(Sp1', Sp2)) + Jz * kron(Sz1, Sz2)
end

# conn refers to the connection operator, that is, the operator on the edge of
# the block, on the interior of the chain.  We need to be able to represent S^z
# and S^+ on that site in the current basis in order to grow the chain.
initial_block = Block(1, model_d, Dict{Symbol,AbstractMatrix{Float64}}(
    :H => H1,
    :conn_Sz => Sz1,
    :conn_Sp => Sp1,
))

function enlarge_block(block::Block)
    # This function enlarges the provided Block by a single site, returning an
    # EnlargedBlock.
    mblock = block.basis_size
    o = block.operator_dict

    # Create the new operators for the enlarged block.  Our basis becomes a
    # Kronecker product of the Block basis and the single-site basis.  NOTE:
    # `kron` uses the tensor product convention making blocks of the second
    # array scaled by the first.  As such, we adopt this convention for
    # Kronecker products throughout the code.
    I1 = sparse(1.0I, model_d, model_d)
    I_block = sparse(1.0I, mblock, mblock)
    enlarged_operator_dict = Dict{Symbol,AbstractMatrix{Float64}}(
        :H => kron(o[:H], I1) + kron(I_block, H1) + H2(o[:conn_Sz], o[:conn_Sp], Sz1, Sp1),
        :conn_Sz => kron(I_block, Sz1),
        :conn_Sp => kron(I_block, Sp1),
    )

    return EnlargedBlock(block.length + 1,
                         block.basis_size * model_d,
                         enlarged_operator_dict)
end

function rotate_and_truncate(operator, transformation_matrix)
    # Transforms the operator to the new (possibly truncated) basis given by
    # `transformation_matrix`.
    return transformation_matrix' * (operator * transformation_matrix)
end

function single_dmrg_step(sys::Block, env::Block, m::Int)
    # Performs a single DMRG step using `sys` as the system and `env` as the
    # environment, keeping a maximum of `m` states in the new basis.

    @assert isvalid(sys)
    @assert isvalid(env)

    # Enlarge each block by a single site.
    sys_enl = enlarge_block(sys)
    if sys === env  # no need to recalculate a second time
        env_enl = sys_enl
    else
        env_enl = enlarge_block(env)
    end

    @assert isvalid(sys_enl)
    @assert isvalid(env_enl)

    # Construct the full superblock Hamiltonian.
    m_sys_enl = sys_enl.basis_size
    m_env_enl = env_enl.basis_size
    sys_enl_op = sys_enl.operator_dict
    env_enl_op = env_enl.operator_dict
    I_sys_enl = sparse(1.0I, m_sys_enl, m_sys_enl)
    I_env_enl = sparse(1.0I, m_env_enl, m_env_enl)
    superblock_hamiltonian = kron(sys_enl_op[:H], I_env_enl) + kron(I_sys_enl, env_enl_op[:H]) +
                             H2(sys_enl_op[:conn_Sz], sys_enl_op[:conn_Sp], env_enl_op[:conn_Sz], env_enl_op[:conn_Sp])

    # Call ARPACK to find the superblock ground state.  (:SR means find the
    # eigenvalue with the "smallest real" value.)
    #
    # But first, we explicitly modify the matrix so that it will be detected as
    # Hermitian by `eigs`.  (Without this step, the matrix is effectively
    # Hermitian but won't be detected as such due to small roundoff error.)
    superblock_hamiltonian = (superblock_hamiltonian + superblock_hamiltonian') / 2
    (energy,), psi0 = eigs(superblock_hamiltonian, nev=1, which=:SR)

    # Construct the reduced density matrix of the system by tracing out the
    # environment
    #
    # We want to make the (sys, env) indices correspond to (row, column) of a
    # matrix, respectively.  Since the environment (column) index updates most
    # quickly in our Kronecker product structure, psi0 is thus row-major.
    # However, Julia stores matrices in column-major format, so we first
    # construct our matrix in (env, sys) form and then take the transpose.
    psi0 = transpose(reshape(psi0, (env_enl.basis_size, sys_enl.basis_size)))
    rho = Hermitian(psi0 * psi0')

    # Diagonalize the reduced density matrix and sort the eigenvectors by
    # eigenvalue.
    fact = eigen(rho)
    evals, evecs = fact.values, fact.vectors
    permutation = sortperm(evals, rev=true)

    # Build the transformation matrix from the `m` overall most significant
    # eigenvectors.
    my_m = min(length(evals), m)
    indices = permutation[1:my_m]
    transformation_matrix = evecs[:, indices]

    truncation_error = 1 - sum(evals[indices])
    println("truncation error: ", truncation_error)

    # Rotate and truncate each operator.
    new_operator_dict = Dict{Symbol,AbstractMatrix{Float64}}()
    for (name, op) in sys_enl.operator_dict
        new_operator_dict[name] = rotate_and_truncate(op, transformation_matrix)
    end

    newblock = Block(sys_enl.length, my_m, new_operator_dict)

    return newblock, energy
end

function graphic(sys_block::Block, env_block::Block, sys_label::Symbol=:l)
    # Returns a graphical representation of the DMRG step we are about to
    # perform, using '=' to represent the system sites, '-' to represent the
    # environment sites, and '**' to represent the two intermediate sites.
    str = repeat("=", sys_block.length) * "**" * repeat("-", env_block.length)
    if sys_label == :r
        # The system should be on the right and the environment should be on
        # the left, so reverse the graphic.
        str = reverse(str)
    elseif sys_label != :l
        throw(ArgumentError("sys_label must be :l or :r"))
    end
    return str
end

function infinite_system_algorithm(L::Int, m::Int)
    block = initial_block
    # Repeatedly enlarge the system by performing a single DMRG step, using a
    # reflection of the current block as the environment.
    while 2 * block.length < L
        println("L = ", block.length * 2 + 2)
        block, energy = single_dmrg_step(block, block, m)
        println("E/L = ", energy / (block.length * 2))
    end
end

function finite_system_algorithm(L::Int, m_warmup::Int, m_sweep_list::AbstractVector{Int})
    @assert iseven(L)

    # To keep things simple, this dictionary is not actually saved to disk, but
    # we use it to represent persistent storage.
    block_disk = Dict{Tuple{Symbol,Int},Block}()  # "disk" storage for Block objects

    # Use the infinite system algorithm to build up to desired size.  Each time
    # we construct a block, we save it for future reference as both a left
    # (:l) and right (:r) block, as the infinite system algorithm assumes the
    # environment is a mirror image of the system.
    block = initial_block
    block_disk[:l, block.length] = block
    block_disk[:r, block.length] = block
    while 2 * block.length < L
        # Perform a single DMRG step and save the new Block to "disk"
        println(graphic(block, block))
        block, energy = single_dmrg_step(block, block, m_warmup)
        println("E/L = ", energy / (block.length * 2))
        block_disk[:l, block.length] = block
        block_disk[:r, block.length] = block
    end

    # Now that the system is built up to its full size, we perform sweeps using
    # the finite system algorithm.  At first the left block will act as the
    # system, growing at the expense of the right block (the environment), but
    # once we come to the end of the chain these roles will be reversed.
    sys_label, env_label = :l, :r

    # Rename block -> sys_block
    sys_block = block
    block = Block(0, 0, Dict{Symbol,AbstractMatrix{Float64}}())

    for m in m_sweep_list
        while true
            # Load the appropriate environment block from "disk"
            env_block = block_disk[env_label, L - sys_block.length - 2]
            if env_block.length == 1
                # We've come to the end of the chain, so we reverse course.
                sys_block, env_block = env_block, sys_block
                sys_label, env_label = env_label, sys_label
            end

            # Perform a single DMRG step.
            println(graphic(sys_block, env_block, sys_label))
            sys_block, energy = single_dmrg_step(sys_block, env_block, m)

            println("E/L = ", energy / L)

            # Save the block from this step to disk.
            block_disk[sys_label, sys_block.length] = sys_block

            # Check whether we just completed a full sweep.
            if sys_label == :l && 2 * sys_block.length == L
                break  # escape from the "while true" loop
            end
        end
    end
end

#infinite_system_algorithm(100, 20)
finite_system_algorithm(20, 10, [10, 20, 30, 40, 40])

max_stars_repo_name: simple-dmrg/simple-dmrg


text: println.(map(fb, 1:100))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode