text:
#jacobigrowth2.jl

using Distributions
function jacobigrowth2(m1,m2,p,beta)
   CJ=eye(p)
   SJ=zeros(p, p)
   for k=1:m2
       Ju = sqrt(rand(Chisq(beta),1,p))
       Ju *= sqrt( rand(Chisq(beta*p)) / rand(Chisq(beta*(k+m1-p))) ) / norm(Ju)
       _, _, _, CJ, SJ=svd(CJ,[SJ; Ju])
   end
   diag(CJ)
end


text:
##
# This file calcules erfc from the representation as a
# cauchy transform
##

using ApproxFun, SingularIntegralEquations

f=Fun(z->2exp(z^2),PeriodicLine(0.,π/2))

erfc2=z->real(z)>0?-exp(-z^2)*cauchy(f,z):exp(-z^2)*(2-cauchy(f,z))

erfc2(1.)-erfc(1.)

ext: jl


text:
@valueequals mutable struct Decal
    texture::String

    x::Number
    y::Number

    scaleX::Integer
    scaleY::Integer

    Decal(texture::String, x::Number, y::Number, scaleX::Integer=1, scaleY::Integer=1) = new(texture, x, y, scaleX, scaleY)
end

function Base.Dict(d::Decal)
    return Dict{String, Any}(
        "texture" => d.texture,

        "x" => d.x,
        "y" => d.y,

        "scaleX" => d.scaleX,
        "scaleY" => d.scaleY,

        "__name" => "decal",
    )
end

max_stars_repo_name: CelestialCartographers/Maple

max_issues_repo_name: CelestialCartographers/Maple

max_forks_repo_name: CelestialCartographers/Maple


text:
"""
    rewrite_from_left(u::AbstractWord, rewriting)
Rewrites word `u` (from left) using `rewriting` object. The object must implement
`rewrite_from_left!(v::AbstractWord, w::AbstractWord, rewriting)`.
"""
function rewrite_from_left(u::W, rewriting,
    vbuff = BufferWord{T}(0, length(u)),
    wbuff = BufferWord{T}(length(u), 0)
) where {T, W<:AbstractWord{T}}
    isempty(rewriting) && return u
    store!(wbuff, u)
    v = rewrite_from_left!(vbuff, wbuff, rewriting)
    return W(v)
end

"""
    rewrite_from_left!(v::AbstractWord, w::AbstractWord, ::Any)
Trivial rewrite: word `w` is simply stored (copied) to `v`.
"""
rewrite_from_left!(v::AbstractWord, w::AbstractWord, ::Any) = store!(v, w)

"""
    rewrite_from_left!(v::AbstractWord, w::AbstractWord, rule::Rule)
Rewrite word `w` storing the result in `v` by using a single rewriting `rule`.
"""
function rewrite_from_left!(v::AbstractWord, w::AbstractWord, rule::Rule)
    v = resize!(v, 0)
    lhs, rhs = rule
    while !isone(w)
        push!(v, popfirst!(w))
        if issuffix(lhs, v)
            prepend!(w, rhs)
            resize!(v, length(v) - length(lhs))
        end
    end
    return v
end

"""
    rewrite_from_left!(v::AbstractWord, w::AbstractWord, A::Alphabet)
Rewrite word `w` storing the result in `v` by applying free reductions as
defined by the inverses present in alphabet `A`.
"""
function rewrite_from_left!(v::AbstractWord, w::AbstractWord, A::Alphabet)
    v = resize!(v, 0)
    while !isone(w)
        if isone(v)
            push!(v, popfirst!(w))
        else
            # the first check is for monoids only
            if hasinverse(last(v), A) && inv(A, last(v)) == first(w)
                pop!(v)
                popfirst!(w)
            else
                push!(v, popfirst!(w))
            end
        end
    end
    return v
end

abstract type AbstractRewritingSystem{W, O} end

"""
    RewritingSystem{W<:AbstractWord, O<:WordOrdering}
RewritingSystem written as a list of Rules (ordered pairs) of `Word`s together with the ordering.
"""
struct RewritingSystem{W<:AbstractWord, O<:WordOrdering} <: AbstractRewritingSystem{W, O}
    rwrules::Vector{Rule{W}}
    order::O
end

function RewritingSystem(rwrules::Vector{Pair{W,W}}, order::O; bare=false) where
    {W<:AbstractWord, O<:WordOrdering}
    @assert length(alphabet(order)) <= _max_alphabet_length(W) "Type $W can not store words over $(alphabet(order))."

    # add rules from the alphabet
    rls = bare ? Rule{W}[] : rules(W, order)
    # properly orient rwrules
    append!(rls, [Rule{W}(simplifyrule!(deepcopy(a), deepcopy(b),order, balance=true)..., order) for (a, b) in rwrules])

    return RewritingSystem(rls, order)
end

rules(s::RewritingSystem) = Iterators.filter(isactive, s.rwrules)
ordering(s::RewritingSystem) = s.order
alphabet(s::RewritingSystem) = alphabet(ordering(s))

Base.push!(s::RewritingSystem{W}, r::Rule{W}) where W = (push!(s.rwrules, r); s)
Base.empty!(s::RewritingSystem) = (empty!(s.rwrules); s)
Base.empty(s::RewritingSystem{W},o::WordOrdering = ordering(s)) where W =
    RewritingSystem(Rule{W}[], o)
Base.isempty(s::RewritingSystem) = isempty(rules(s))

"""
    rewrite_from_left!(v::AbstractWord, w::AbstractWord, rws::RewritingSystem)
Rewrite word `w` storing the result in `v` by left using rewriting rules of
rewriting system `rws`. See [Sims, p.66]
"""
function rewrite_from_left!(v::AbstractWord, w::AbstractWord, rws::RewritingSystem)
    v = resize!(v, 0)
    while !isone(w)
        push!(v, popfirst!(w))
        for (lhs, rhs) in rules(rws)
            if issuffix(lhs, v)
                prepend!(w, rhs)
                resize!(v, length(v) - length(lhs))
            end
        end
    end
    return v
end

"""
    isirreducible(w::AbstractWord, rws::RewritingSystem)
Returns whether a word is irreducible with respect to a given rewriting system
"""
isirreducible(w::AbstractWord, rws::RewritingSystem) =
    !any(r -> occursin(first(r), w), rules(rws))

"""
    subwords(w::AbstractWord[, minlength=1, maxlength=length(w)])
Return an iterator over all `SubWord`s of `w` of length between `minlength` and `maxlength`.
"""
function subwords(w::AbstractWord, minlength=1, maxlength=length(w))
    n = length(w)
    (@view(w[i:j]) for i in 1:n for j in i:n if
        minlength <= j-i+1 <= maxlength)
end

"""
    irreduciblesubsystem(rws::RewritingSystem)
Return an array of left sides of rules from rewriting system of which all the
proper subwords are irreducible with respect to this rewriting system.
"""
function irreduciblesubsystem(rws::RewritingSystem{W}) where W
    lsides = W[]
    for rule in rws.rwrules
        lhs = first(rule)
        length(lhs) >= 2 || break
        for sw in subwords(lhs, 2, length(lhs)-1)
            if !isirreducible(sw, rws)
                @debug "subword $sw of $lhs is reducible. skipping!"
                break
            end
        end
        if all(sw->isirreducible(sw, rws), subwords(lhs, 2, length(lhs)-1))
            @debug "all subwords are irreducible; pushing $lhs"
            push!(lsides, lhs)
        end
    end
    return unique!(lsides)
end

"""
    simplifyrule!(lhs::AbstractWord, rhs::AbstractWord, A::Alphabet)
Simplifies both sides of the rule if they start/end with the same invertible word.
"""
function simplifyrule!(lhs::AbstractWord, rhs::AbstractWord, A::Alphabet)
    common_suffix=0
    k = min(length(lhs), length(rhs))
    @inbounds for i in 0:k-1
        l,r = lhs[end-i], rhs[end-i]
        l != r && break
        hasinverse(l, A) || break
        common_suffix += 1
    end

    if !iszero(common_suffix)
        resize!(lhs, length(lhs) - common_suffix)
        resize!(rhs, length(rhs) - common_suffix)
    end

    common_prefix=0
    for (l, r) in zip(lhs,rhs)
        l != r && break
        hasinverse(l, A) || break
        common_prefix += 1
    end

    if !iszero(common_prefix)
        copyto!(lhs, 1, lhs, common_prefix + 1, length(lhs) - common_prefix)
        copyto!(rhs, 1, rhs, common_prefix + 1, length(rhs) - common_prefix)
        resize!(lhs, length(lhs) - common_prefix)
        resize!(rhs, length(rhs) - common_prefix)
    end

    return lhs, rhs
end

function balancerule!(lhs::AbstractWord, rhs::AbstractWord, A::Alphabet)
    while length(lhs) > 2 && length(lhs) > length(rhs)
        hasinverse(last(lhs), A) || break
        push!(rhs, inv(A, pop!(lhs)))
    end

    while length(lhs) > 2 && length(lhs) > length(rhs)
        hasinverse(first(lhs), A) || break
        pushfirst!(rhs, inv(A, popfirst!(lhs)))
    end

    return lhs, rhs
end

function simplifyrule!(lhs::AbstractWord, rhs::AbstractWord, o::Ordering; balance=false)

    lhs, rhs = simplifyrule!(lhs, rhs, alphabet(o))
    if balance
        lhs, rhs = balancerule!(lhs, rhs, alphabet(o))
    end

    return lhs, rhs
end

function Base.show(io::IO, rws::RewritingSystem)
    rls = collect(rules(rws))
    println(io, "Rewriting System with $(length(rls)) active rules ordered by $(ordering(rws)):")
    height = first(displaysize(io))
    A = alphabet(rws)
    if height > length(rls)
        for (i, rule) in enumerate(rls)
            _print_rule(io, i, rule, A)
        end
    else
        for i in 1:height-5
            rule = rls[i]
            _print_rule(io, i, rule, A)
        end

        println(io, "⋮")
        for i in (length(rls)-4):length(rls)
            rule = rls[i]
            _print_rule(io, i, rule, A)
        end
    end
end

function _print_rule(io::IO, i, rule, A)
    (lhs, rhs) = rule
    print(io, i, ". ")
    print_repr(io, lhs, A)
    print(io, "\t → \t")
    print_repr(io, rhs, A)
    println(io, "")
end


text:
using Base.Iterators
using Test

include("../src/Gillespie.jl")
include("../src/Model_CRS.jl")
include("../src/Droplets.jl")
include("../src/Data.jl")

@testset "CRS helpers" begin

    @test Set(generate_all_binary_seqs(2)) == Set(["A", "B", "AA","AB","BA", "BB"])
    @test Set(generate_all_binary_seqs(3)) == Set(["A", "B", "AA","AB","BA", "BB", "AAA", "AAB","ABA","BAA","BBA","BAB","ABB","BBB"])

    @test all_splits("12345") == (["1","12", "123","1234"], ["2345","345", "45", "5"])

    @test countall("10", "1") == 1
    @test countall("10", "0") == 1
    @test countall("111", "1") == 3
    @test countall("101", "1") == 2
    @test countall("BA", "1") == 0
    @test countall("AB", "B") == 1

    @test reduced_mass(1,1) == 0.5
    @test isapprox(reduced_mass(10000,1),1, atol = 1E-4)

    @test calculate_reduced_mass("AAA", "AB", 1, 2) == 1.5
    @test isapprox(calculate_reduced_mass("AAA", "BB", 1, 2), 1.71428571429)

    @test bimolecular_coef(1,1,1,1,1) == 1.0
    @test isapprox(bimolecular_coef(0.001, 300, 8.5, 1000, 10), 0.00001596871942267131199907024517698) # k = 0.001, T = 300, R = 8.5, V = 1000, m = 10
end

@testset "Reaction Funcs" begin

    @test PickReactionID([1.0, 0.0, 0.0]) == 1
    @test PickReactionID([0.0, 0.0, 1.0]) == 3
    @test PickReactionID([0.0, 0.0001, 0.0]) == 2
    ### Include som fancier tests here to make sure relative proportions are correct

end

max_stars_repo_name: croningp/clusterautocatalyticsets

max_issues_repo_name: croningp/clusterautocatalyticsets

max_forks_repo_name: croningp/clusterautocatalyticsets


text:
using DashBootstrapComponents

popovers = html_div([
    dbc_button("Click Me", id = "component-target", n_clicks = 0),
    dbc_popover(
        [
            dbc_popoverheader("Popover header"),
            dbc_popoverbody("And here's some amazing content. Cool!"),
        ],
        target = "component-target",
        trigger = "click",
    ),
]);

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
x = 1

for i in 1:3
    global x+=i
end


text:
#%%
using GdMn2O5
using Parameters
using Base.Threads
using DataFrames
using FileIO
using DrWatson
#%%
df = load(datadir("sims.jld2"))["df"]

model_params_to_test = Dict(
    :K_Gd => [range(0, 1, length=3); range(2, 20, length=3)],
    :K_L  => collect(range(0,1,length=3)),
    :γ    => collect(range(0,2, length=6)),
    :J1   => collect(range(0, 100, length=3)),
    :J2   => 1.0,
    :χ_Gd => 25.0,
    :χ_L  => 1/10,
    :g    => 1.0,
    :α    => 2.5,
    :dipole_strength => collect(range(0,10, length=3)))

function add_simulations!(df::DataFrame, model_params_to_test::Dict)
	mp_to_test = dict_list(model_params_to_test)
	ledgers = [Ledger(Model(;mp_to_test[1]...)) for i=1:nthreads()]
    for m in ledgers
        Entity(m, H(ϕ=0.0))
    end
    for params in mp_to_test
        model = Model(;params...)
        modeldict = type2dict(model)
        xs = 0:2:16
        hmags = repeat([range(0.0, 6, length=50); range(6, 0.0, length=50)], 3)
        for x in xs
            Hϕ = deg2rad(x)
            Pbs  = similar(hmags)
            L1ϕs = similar(hmags)
            L2ϕs = similar(hmags)
            ϕ1s  = similar(hmags)
            ϕ2s  = similar(hmags)
            ϕ3s  = similar(hmags)
            ϕ4s  = similar(hmags)
            ϕ5s  = similar(hmags)
            ϕ6s  = similar(hmags)
            ϕ7s  = similar(hmags)
            ϕ8s  = similar(hmags)
            Es   = similar(hmags)
            @threads for i in 1:length(hmags)
                h = hmags[i]
                m = ledgers[threadid()]
                m[H].data[1] = H(ϕ=Hϕ, magnitude=h)
                # p5 = plot_Lheatmap(m, range(0, 2π, length=20))
                m[Spin].data[1:8] = [Spin(ϕ=x) for x in rand(0:0.01:2π,8)]

                # @time push!(energies, optimize(m, Spin, L).minimum)
                E =  GdMn2O5.best_fitness(optimize(m, Spin, L))
                Pbs[i]  = Pb(m)
                L1ϕs[i] = mod(m[L][1].ϕ, 2π)
                L2ϕs[i] = mod(m[L][1].ϕ, 2π)
                ϕ1s[i]  = mod(m[Spin][1].ϕ, 2π)
                ϕ2s[i]  = mod(m[Spin][2].ϕ, 2π)
                ϕ3s[i]  = mod(m[Spin][3].ϕ, 2π)
                ϕ4s[i]  = mod(m[Spin][4].ϕ, 2π)
                ϕ5s[i]  = mod(m[Spin][5].ϕ, 2π)
                ϕ6s[i]  = mod(m[Spin][6].ϕ, 2π)
                ϕ7s[i]  = mod(m[Spin][7].ϕ, 2π)
                ϕ8s[i]  = mod(m[Spin][8].ϕ, 2π)
                Es[i]   = E
            end
            push!(df, merge(modeldict, Dict(:Pbs  =>Pbs,
                                           :L1ϕs =>L1ϕs,
                                           :L2ϕs =>L2ϕs,
                                           :ϕ1s  =>ϕ1s,
                                           :ϕ2s  =>ϕ2s,
                                           :ϕ3s  =>ϕ3s,
                                           :ϕ4s  =>ϕ4s,
                                           :ϕ5s  =>ϕ5s,
                                           :ϕ6s  =>ϕ6s,
                                           :ϕ7s  =>ϕ7s,
                                           :ϕ8s  =>ϕ8s,
                                           :Es   =>Es,
                                           :Hϕ  =>Hϕ,
                                           :Hmags =>hmags)))
        end
    end
    return df
end
add_simulations!(df, model_params_to_test)
save(datadir("sims.jld2"), "df", df)
#%%


text:
#
# Shared Formulation Definitions
#################################
#
# This is the home of functions that are shared across multiple branches
# of the type hierarchy.  Hence all function in this file should be over
# union types.
#
# The types defined in this file should not be exported because they exist
# only to prevent code replication
#
# Note that Union types are discouraged in Julia,
# https://docs.julialang.org/en/release-0.6/manual/style-guide/#Avoid-strange-type-Unions-1
# and should be used with discretion.
#
# If you are about to add a union type,
# first double check if a different type hierarchy can resolve the issue
# instead.
#

AbstractWRForms = Union{AbstractACTForm, AbstractWRForm, AbstractWRMForm}
AbstractPForms = Union{AbstractACPForm, AbstractACTForm, AbstractDCPForm}

"`t[ref_bus] == 0`"
function constraint_theta_ref(pm::GenericPowerModel{T}, n::Int, i::Int) where T <: AbstractPForms
    @constraint(pm.model, pm.var[:nw][n][:va][i] == 0)
end

"""
```
t[f_bus] - t[t_bus] <= angmax
t[f_bus] - t[t_bus] >= angmin
```
"""
function constraint_voltage_angle_difference(pm::GenericPowerModel{T}, n::Int, f_bus, t_bus, angmin, angmax) where T <: AbstractPForms
    va_fr = pm.var[:nw][n][:va][f_bus]
    va_to = pm.var[:nw][n][:va][t_bus]

    @constraint(pm.model, va_fr - va_to <= angmax)
    @constraint(pm.model, va_fr - va_to >= angmin)
end

function constraint_voltage_magnitude_setpoint(pm::GenericPowerModel{T}, n::Int, i, vm) where T <: AbstractWRForms
    w = pm.var[:nw][n][:w][i]

    @constraint(pm.model, w == vm^2)
end

"Do nothing, no way to represent this in these variables"
function constraint_theta_ref(pm::GenericPowerModel{T}, n::Int, ref_bus::Int) where T <: AbstractWRForms
end

"""
```
sum(p[a] for a in bus_arcs) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == sum(pg[g] for g in bus_gens) - pd - gs*w[i]
sum(q[a] for a in bus_arcs) + sum(q_dc[a_dc] for a_dc in bus_arcs_dc) == sum(qg[g] for g in bus_gens) - qd + bs*w[i]
```
"""
function constraint_kcl_shunt(pm::GenericPowerModel{T}, n::Int, i, bus_arcs, bus_arcs_dc, bus_gens, pd, qd, gs, bs) where T <: AbstractWRForms
    w = pm.var[:nw][n][:w][i]
    pg = pm.var[:nw][n][:pg]
    qg = pm.var[:nw][n][:qg]
    p = pm.var[:nw][n][:p]
    q = pm.var[:nw][n][:q]
    p_dc = pm.var[:nw][n][:p_dc]
    q_dc = pm.var[:nw][n][:q_dc]

    @constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == sum(pg[g] for g in bus_gens) - pd - gs*w)
    @constraint(pm.model, sum(q[a] for a in bus_arcs) + sum(q_dc[a_dc] for a_dc in bus_arcs_dc) == sum(qg[g] for g in bus_gens) - qd + bs*w)
end

"""
```
sum(p[a] for a in bus_arcs) + sum(p_ne[a] for a in bus_arcs_ne) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == sum(pg[g] for g in bus_gens) - pd - gs*w[i]
sum(q[a] for a in bus_arcs) + sum(q_ne[a] for a in bus_arcs_ne) + sum(q_dc[a_dc] for a_dc in bus_arcs_dc) == sum(qg[g] for g in bus_gens) - qd + bs*w[i]
```
"""
function constraint_kcl_shunt_ne(pm::GenericPowerModel{T}, n::Int, i, bus_arcs, bus_arcs_dc, bus_arcs_ne, bus_gens, pd, qd, gs, bs) where T <: AbstractWRForms
    w = pm.var[:nw][n][:w][i]
    pg = pm.var[:nw][n][:pg]
    qg = pm.var[:nw][n][:qg]
    p = pm.var[:nw][n][:p]
    q = pm.var[:nw][n][:q]
    p_ne = pm.var[:nw][n][:p_ne]
    q_ne = pm.var[:nw][n][:q_ne]
    p_dc = pm.var[:nw][n][:p_dc]
    q_dc = pm.var[:nw][n][:q_dc]

    @constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_ne[a] for a in bus_arcs_ne) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == sum(pg[g] for g in bus_gens) - pd - gs*w)
    @constraint(pm.model, sum(q[a] for a in bus_arcs) + sum(q_ne[a] for a in bus_arcs_ne) + sum(q_dc[a_dc] for a_dc in bus_arcs_dc) == sum(qg[g] for g in bus_gens) - qd + bs*w)
end

"""
Creates Ohms constraints (yt post fix indicates that Y and T values are in rectangular form)
"""
function constraint_ohms_yt_from(pm::GenericPowerModel{T}, n::Int, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm) where T <: AbstractWRForms
    p_fr = pm.var[:nw][n][:p][f_idx]
    q_fr = pm.var[:nw][n][:q][f_idx]
    w_fr = pm.var[:nw][n][:w][f_bus]
    wr = pm.var[:nw][n][:wr][(f_bus, t_bus)]
    wi = pm.var[:nw][n][:wi][(f_bus, t_bus)]

    @constraint(pm.model, p_fr ==        g/tm^2*w_fr + (-g*tr+b*ti)/tm^2*wr + (-b*tr-g*ti)/tm^2*wi )
    @constraint(pm.model, q_fr == -(b+c/2)/tm^2*w_fr - (-b*tr-g*ti)/tm^2*wr + (-g*tr+b*ti)/tm^2*wi )
end

"""
Creates Ohms constraints (yt post fix indicates that Y and T values are in rectangular form)
"""
function constraint_ohms_yt_to(pm::GenericPowerModel{T}, n::Int, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm) where T <: AbstractWRForms
    q_to = pm.var[:nw][n][:q][t_idx]
    p_to = pm.var[:nw][n][:p][t_idx]
    w_to = pm.var[:nw][n][:w][t_bus]
    wr = pm.var[:nw][n][:wr][(f_bus, t_bus)]
    wi = pm.var[:nw][n][:wi][(f_bus, t_bus)]

    @constraint(pm.model, p_to ==        g*w_to + (-g*tr-b*ti)/tm^2*wr + (-b*tr+g*ti)/tm^2*-wi )
    @constraint(pm.model, q_to == -(b+c/2)*w_to - (-b*tr+g*ti)/tm^2*wr + (-g*tr-b*ti)/tm^2*-wi )
end

""
function constraint_voltage_angle_difference(pm::GenericPowerModel{T}, n::Int, f_bus, t_bus, angmin, angmax) where T <: AbstractWRForms
    w_fr = pm.var[:nw][n][:w][f_bus]
    w_to = pm.var[:nw][n][:w][t_bus]
    wr = pm.var[:nw][n][:wr][(f_bus, t_bus)]
    wi = pm.var[:nw][n][:wi][(f_bus, t_bus)]

    @constraint(pm.model, wi <= tan(angmax)*wr)
    @constraint(pm.model, wi >= tan(angmin)*wr)
    cut_complex_product_and_angle_difference(pm.model, w_fr, w_to, wr, wi, angmin, angmax)
end

""
function add_bus_voltage_setpoint(sol, pm::GenericPowerModel{T}) where T <: AbstractWRForms
    add_setpoint(sol, pm, "bus", "vm", :w; scale = (x,item) -> sqrt(x))
    # What should the default value be?
    #add_setpoint(sol, pm, "bus", "va", :va; default_value = 0)
end


text:
# SP Constructs #
# ========================== #
"""
    WS(stochasticprogram::TwoStageStochasticProgram, scenario::AbstractScenarioaData; solver = JuMP.UnsetSolver())

Generate a **wait-and-see** (`WS`) model of the two-stage `stochasticprogram`, corresponding to `scenario`.

In other words, generate the first stage and the second stage of the `stochasticprogram` as if `scenario` is known to occur. Optionally, a capable `solver` can be supplied to `WS`. Otherwise, any previously set solver will be used.

See also: [`DEP`](@ref), [`EVP`](@ref)
"""
function WS(stochasticprogram::StochasticProgram{2}, scenario::AbstractScenario; solver::SPSolverType = JuMP.UnsetSolver())
    # Use cached solver if available
    supplied_solver = pick_solver(stochasticprogram, solver)
    # Check that the required generators have been defined
    has_generator(stochasticprogram, :stage_1) || error("First-stage problem not defined in stochastic program. Consider @stage 1.")
    has_generator(stochasticprogram, :stage_2) || error("Second-stage problem not defined in stochastic program. Consider @stage 2.")
    # Return WS model
    return _WS(generator(stochasticprogram,:stage_1), generator(stochasticprogram,:stage_2), stage_parameters(stochasticprogram, 1), stage_parameters(stochasticprogram, 2), scenario, internal_solver(supplied_solver))
end
function _WS(stage_one_generator::Function,
             stage_two_generator::Function,
             stage_one_params::Any,
             stage_two_params::Any,
             scenario::AbstractScenario,
             solver::MathProgBase.AbstractMathProgSolver)
    ws_model = Model(solver = solver)
    stage_one_generator(ws_model, stage_one_params)
    ws_obj = copy(ws_model.obj)
    stage_two_generator(ws_model, stage_two_params, scenario, ws_model)
    append!(ws_obj, ws_model.obj)
    ws_model.obj = ws_obj
    return ws_model
end
"""
    WS_decision(stochasticprogram::TwoStageStochasticProgram, scenario::AbstractScenario; solver = JuMP.UnsetSolver())

Calculate the optimizer of the **wait-and-see** (`WS`) model of the two-stage `stochasticprogram`, corresponding to `scenario`.

Optionally, supply a capable `solver` to solve the wait-and-see problem. The default behaviour is to rely on any previously set solver.

See also: [`WS`](@ref)
"""
function WS_decision(stochasticprogram::StochasticProgram{2}, scenario::AbstractScenario; solver::SPSolverType = JuMP.UnsetSolver())
    # Solve WS model for supplied scenario
    ws_model = WS(stochasticprogram, scenario, solver = solver)
    solve(ws_model)
    # Return WS decision
    decision = ws_model.colVal[1:decision_length(stochasticprogram)]
    if any(isnan.(decision))
        warn("Optimal decision not defined. Check that the EVP model was properly solved.")
    end
    return decision
end
"""
    EWS(stochasticprogram::StochasticProgram; solver = JuMP.UnsetSolver())

Calculate the **expected wait-and-see result** (`EWS`) of the `stochasticprogram`.

In other words, calculate the expectated result of all possible wait-and-see models, using the provided scenarios in `stochasticprogram`. Optionally, a capable `solver` can be supplied to solve the intermediate problems. Otherwise, any previously set solver will be used.

See also: [`VRP`](@ref), [`WS`](@ref)
"""
function EWS(stochasticprogram::StochasticProgram{2}; solver::SPSolverType = JuMP.UnsetSolver())
    # Use cached solver if available
    supplied_solver = pick_solver(stochasticprogram, solver)
    # Abort if no solver was given
    if isa(supplied_solver, JuMP.UnsetSolver)
        error("Cannot determine EWS without a solver.")
    end
    # Solve all possible WS models and compute EWS
    return _EWS(stochasticprogram, internal_solver(supplied_solver))
end
function _EWS(stochasticprogram::TwoStageStochasticProgram{S,SP}, solver::MathProgBase.AbstractMathProgSolver) where {S, SP <: ScenarioProblems}
    return sum([begin
                ws = _WS(stochasticprogram.generator[:stage_1],
                         stochasticprogram.generator[:stage_2],
                         stage_parameters(stochasticprogram, 1),
                         stage_parameters(stochasticprogram, 2),
                         scenario,
                         solver)
                solve(ws)
                probability(scenario)*getobjectivevalue(ws)
                end for scenario in scenarios(stochasticprogram.scenarioproblems)])
end
function _EWS(stochasticprogram::TwoStageStochasticProgram{S,SP}, solver::MathProgBase.AbstractMathProgSolver) where {S, SP <: DScenarioProblems}
    partial_ews = Vector{Float64}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_ews[i] = remotecall_fetch((sp,stage_one_generator,stage_two_generator,stage_one_params,stage_two_params,solver)->begin
                scenarioproblems = fetch(sp)
                isempty(scenarioproblems.scenarios) && return 0.0
                return sum([begin
                            ws = _WS(stage_one_generator,
                                     stage_two_generator,
                                     stage_one_params,
                                     stage_two_params,
                                     scenario,
                                     solver)
                            solve(ws)
                            probability(scenario)*getobjectivevalue(ws)
                            end for scenario in scenarioproblems.scenarios])
                end,
                w,
                stochasticprogram.scenarioproblems[w-1],
                stochasticprogram.generator[:stage_1],
                stochasticprogram.generator[:stage_2],
                stage_parameters(stochasticprogram, 1),
                stage_parameters(stochasticprogram, 2),
                solver)
        end
    end
    return sum(partial_ews)
end
"""
    SAA(stochasticmodel::StochasticModel, sampler::AbstractSampler, n::Integer; solver = JuMP.UnsetSolver())

Generate a **sample average approximation** (`SAA`) instance of size `n` using the model stored in the two-stage `stochasticmodel`, and the provided `sampler`.

Optionally, a capable `solver` can be supplied to `SAA`. Otherwise, any previously set solver will be used.

See also: [`sample!`](@ref)
"""
function SAA(sm::StochasticModel{2}, sampler::AbstractSampler{S}, n::Integer; solver::SPSolverType = JuMP.UnsetSolver(), procs = workers(), kw...) where S <: AbstractScenario
    # Create new stochastic program instance
    saa = StochasticProgram(parameters(sm.parameters[1]; kw...),
                            parameters(sm.parameters[2]; kw...),
                            S,
                            solver,
                            procs)
    sm.generator(saa)
    # Sample n scenarios
    add_scenarios!(saa, n) do
        return sample(sampler, 1/n)
    end
    # Return the SAA instance
    return saa
end
function SAA(sm::StochasticModel{2}, sampler::AbstractSampler{S}, n::Integer; solver::SPSolverType = JuMP.UnsetSolver(), procs = workers(), kw...) where S <: Scenario
    # Create new stochastic program instance
    saa = StochasticProgram(parameters(sm.parameters[1]; kw...),
                            parameters(sm.parameters[2]; kw...),
                            typeof(sample(sampler)),
                            solver,
                            procs)
    sm.generator(saa)
    # Sample n scenarios
    add_scenarios!(saa, n) do
        return sample(sampler, 1/n)
    end
    # Return the SAA instance
    return saa
end
function SAA(sm::StochasticModel{2}, sampler::AbstractSampler{S}, solution::StochasticSolution; solver::SPSolverType = JuMP.UnsetSolver(), procs = workers(), kw...) where S <: AbstractScenario
    if isa(solver, JuMP.UnsetSolver)
        error("Cannot generate SAA from stochastic solution without a solver.")
    end
    n = 16
    CI = confidence_interval(solution)
    α = 1 - confidence(CI)
    while !(confidence_interval(sm, sampler; solver = solver, N = n, M = M, confidence = 1-α) ⊆ CI)
        n = n * 2
    end
    return SAA(sm, sampler, n; solver = solver, procs = procs)
end
"""
    DEP(stochasticprogram::TwoStageStochasticProgram; solver = JuMP.UnsetSolver())

Generate the **deterministically equivalent problem** (`DEP`) of the two-stage `stochasticprogram`.

In other words, generate the extended form the `stochasticprogram` as a single JuMP model. Optionally, a capable `solver` can be supplied to `DEP`. Otherwise, any previously set solver will be used.

See also: [`VRP`](@ref), [`WS`](@ref)
"""
function DEP(stochasticprogram::StochasticProgram{2}; solver::SPSolverType = JuMP.UnsetSolver())
    # Use cached solver if available
    supplied_solver = pick_solver(stochasticprogram, solver)
    # Return possibly cached model
    cache = problemcache(stochasticprogram)
    if haskey(cache,:dep)
        dep = cache[:dep]
        setsolver(dep, supplied_solver)
        return dep
    end
    # Check that the required generators have been defined
    has_generator(stochasticprogram, :stage_1) || error("First-stage problem not defined in stochastic program. Consider @stage 1.")
    has_generator(stochasticprogram, :stage_2) || error("Second-stage problem not defined in stochastic program. Consider @stage 2.")
    # Define first-stage problem
    dep_model = Model(solver = internal_solver(supplied_solver))
    generator(stochasticprogram,:stage_1)(dep_model, stage_parameters(stochasticprogram, 1))
    dep_obj = copy(dep_model.obj)
    # Define second-stage problems, renaming variables according to scenario.
    stage_two_params = stage_parameters(stochasticprogram, 2)
    visited_objs = collect(keys(dep_model.objDict))
    for (i, scenario) in enumerate(scenarios(stochasticprogram))
        generator(stochasticprogram,:stage_2)(dep_model, stage_two_params, scenario, dep_model)
        append!(dep_obj,probability(scenario)*dep_model.obj)
        for (objkey,obj) ∈ filter(kv->kv.first ∉ visited_objs, dep_model.objDict)
            newkey = if isa(obj,JuMP.Variable)
                varname = add_subscript(dep_model.colNames[obj.col], i)
                dep_model.colNames[obj.col] = varname
                dep_model.colNamesIJulia[obj.col] = varname
                newkey = Symbol(varname)
            elseif isa(obj,Array{JuMP.Variable})
                JuMP.fill_var_names(JuMP.REPLMode, dep_model.colNames, obj)
                arrayname = add_subscript(objkey, i)
                for var in obj
                    splitname = split(dep_model.colNames[var.col],"[")
                    varname = @sprintf("%s[%s", add_subscript(splitname[1],i), splitname[2])
                    dep_model.colNames[var.col] = varname
                    dep_model.colNamesIJulia[var.col] = varname
                end
                newkey = Symbol(arrayname)
            elseif isa(obj,JuMP.ConstraintRef)
                arrayname = add_subscript(objkey, i)
                newkey = Symbol(arrayname)
            elseif isa(obj,Array{ConstraintRef})
                arrayname = add_subscript(objkey, i)
                newkey = Symbol(arrayname)
            elseif isa(obj,JuMP.JuMPArray)
                newkey = if isa(obj,JuMP.JuMPArray{JuMP.ConstraintRef})
                    arrayname = add_subscript(objkey, i)
                    newkey = Symbol(arrayname)
                else
                    JuMP.fill_var_names(JuMP.REPLMode, dep_model.colNames, obj)
                    arrayname = add_subscript(dep_model.varData[obj].name, i)
                    newkey = Symbol(arrayname)
                    dep_model.varData[obj].name = newkey
                    for var in obj.innerArray
                        splitname = split(dep_model.colNames[var.col],"[")
                        varname = @sprintf("%s[%s", add_subscript(splitname[1],i), splitname[2])
                        dep_model.colNames[var.col] = varname
                        dep_model.colNamesIJulia[var.col] = varname
                    end
                    newkey
                end
            else
                continue
            end
            dep_model.objDict[newkey] = obj
            delete!(dep_model.objDict,objkey)
            push!(visited_objs,newkey)
        end
    end
    dep_model.obj = dep_obj
    # Cache DEP
    cache[:dep] = dep_model
    # Return DEP
    return dep_model
end
"""
    VRP(stochasticprogram::StochasticProgram; solver = JuMP.UnsetSolver())

Calculate the **value of the recouse problem** (`VRP`) in `stochasticprogram`.

In other words, optimize the stochastic program and return the optimal value. Optionally, supply a capable `solver` to optimize the stochastic program. Otherwise, any previously set solver will be used.

See also: [`EVPI`](@ref), [`EWS`](@ref)
"""
function VRP(stochasticprogram::StochasticProgram; solver::SPSolverType = JuMP.UnsetSolver())
    # Solve DEP
    optimize!(stochasticprogram, solver = solver)
    # Return optimal value
    return optimal_value(stochasticprogram)
end
"""
    EVPI(stochasticprogram::TwoStageStochasticProgram; solver = JuMP.UnsetSolver())

Calculate the **expected value of perfect information** (`EVPI`) of the two-stage `stochasticprogram`.

In other words, calculate the gap between `VRP` and `EWS`. Optionally, supply a capable `solver` to solve the intermediate problems. Otherwise, any previously set solver will be used.

See also: [`VRP`](@ref), [`EWS`](@ref), [`VSS`](@ref)
"""
function EVPI(stochasticprogram::StochasticProgram{2}; solver::SPSolverType = JuMP.UnsetSolver())
    # Use cached solver if available
    supplied_solver = pick_solver(stochasticprogram, solver)
    # Abort if no solver was given
    if isa(supplied_solver, JuMP.UnsetSolver)
        error("Cannot determine EVPI without a solver.")
    end
    # Calculate VRP
    vrp = VRP(stochasticprogram, solver=supplied_solver)
    # Solve all possible WS models and calculate EWS
    ews = _EWS(stochasticprogram, internal_solver(supplied_solver))
    # Return EVPI = EWS-VRP
    return abs(ews-vrp)
end
"""
    EVP(stochasticprogram::TwoStageStochasticProgram; solver = JuMP.UnsetSolver())

Generate the **expected value problem** (`EVP`) of the two-stage `stochasticprogram`.

In other words, generate a wait-and-see model corresponding to the expected scenario over all available scenarios in `stochasticprogram`. Optionally, supply a capable `solver` to `EVP`. Otherwise, any previously set solver will be used.

See also: [`EVP_decision`](@ref), [`EEV`](@ref), [`EV`](@ref), [`WS`](@ref)
"""
function EVP(stochasticprogram::StochasticProgram{2}; solver::SPSolverType = JuMP.UnsetSolver())
    # Use cached solver if available
    supplied_solver = pick_solver(stochasticprogram, solver)
    # Return possibly cached model
    cache = problemcache(stochasticprogram)
    if haskey(cache,:evp)
        evp = cache[:evp]
        setsolver(evp, supplied_solver)
        return evp
    end
    # Create EVP as a wait-and-see model of the expected scenario
    ev_model = WS(stochasticprogram, expected(stochasticprogram), solver = solver)
    # Cache EVP
    cache[:evp] = ev_model
    # Return EVP
    return ev_model
end
"""
    EVP_decision(stochasticprogram::TwoStageStochasticProgram; solver = JuMP.UnsetSolver())

Calculate the optimizer of the `EVP` of the two-stage `stochasticprogram`.

Optionally, supply a capable `solver` to solve the expected value problem. The default behaviour is to rely on any previously set solver.

See also: [`EVP`](@ref), [`EV`](@ref), [`EEV`](@ref)
"""
function EVP_decision(stochasticprogram::StochasticProgram{2}; solver::SPSolverType = JuMP.UnsetSolver())
    # Solve EVP
    evp = EVP(stochasticprogram, solver = solver)
    solve(evp)
    # Return EVP decision
    decision = evp.colVal[1:decision_length(stochasticprogram)]
    if any(isnan.(decision))
        warn("Optimal decision not defined. Check that the EVP model was properly solved.")
    end
    return decision
end
"""
    EV(stochasticprogram::TwoStageStochasticProgram; solver = JuMP.UnsetSolver())

Calculate the optimal value of the `EVP` of the two-stage `stochasticprogram`.

Optionally, supply a capable `solver` to solve the expected value problem. The default behaviour is to rely on any previously set solver.

See also: [`EVP`](@ref), [`EVP_decision`](@ref), [`EEV`](@ref)
"""
function EV(stochasticprogram::StochasticProgram{2}; solver::SPSolverType = JuMP.UnsetSolver())
    # Solve EVP model
    evp = EVP(stochasticprogram; solver = solver)
    solve(evp)
    # Return optimal value
    return getobjectivevalue(evp)
end
"""
    EEV(stochasticprogram::TwoStageStochasticProgram; solver = JuMP.UnsetSolver())

Calculate the **expected value of using the expected value solution** (`EEV`) of the two-stage `stochasticprogram`.

In other words, evaluate the `EVP` decision. Optionally, supply a capable `solver` to solve the intermediate problems. The default behaviour is to rely on any previously set solver.

See also: [`EVP`](@ref), [`EV`](@ref)
"""
function EEV(stochasticprogram::StochasticProgram{2}; solver::SPSolverType = JuMP.UnsetSolver())
    # Solve EVP model
    evp_decision = EVP_decision(stochasticprogram; solver = solver)
    # Calculate EEV by evaluating the EVP decision
    eev = evaluate_decision(stochasticprogram, evp_decision; solver = solver)
    # Return EEV
    return eev
end
"""
    VSS(stochasticprogram::TwoStageStochasticProgram; solver = JuMP.UnsetSolver())

Calculate the **value of the stochastic solution** (`VSS`) of the two-stage `stochasticprogram`.

In other words, calculate the gap between `EEV` and `VRP`. Optionally, supply a capable `solver` to solve the intermediate problems. The default behaviour is to rely on any previously set solver.
"""
function VSS(stochasticprogram::StochasticProgram{2}; solver::SPSolverType = JuMP.UnsetSolver())
    # Solve EVP and determine EEV
    eev = EEV(stochasticprogram; solver = solver)
    # Calculate VRP
    vrp = VRP(stochasticprogram; solver = solver)
    # Return VSS = VRP-EEV
    return abs(vrp-eev)
end
# ========================== #


text:
############################################################################
### Expectation Maximization Algorithm
############################################################################

# An EM Algorithm for MVN-distributed Data with missing values
# Adapted from supplementary Material to the book Machine Learning: A Probabilistic Perspective
# Copyright (2010) Kevin Murphy and Matt Dunham
# found at https://github.com/probml/pmtk3/blob/master/toolbox/BasicModels/gauss/sub/gaussMissingFitEm.m
# and at https://github.com/probml/pmtk3/blob/master/toolbox/Algorithms/optimization/emAlgo.m

# what about random restarts?

# outer function -----------------------------------------------------------------

function em_mvn(
    observed::SemObsMissing;
    start_em = start_em_observed,
    max_iter_em = 100,
    rtol_em = 1e-4,
    kwargs...)

    n_obs, n_man = observed.n_obs, Int(observed.n_man)

    # preallocate stuff?
    𝔼x_pre = zeros(n_man)
    𝔼xxᵀ_pre = zeros(n_man, n_man)

    ### precompute for full cases
    if length(observed.patterns[1]) == observed.n_man
        for row ∈ observed.rows[1]
            row = observed.data_rowwise[row]
            𝔼x_pre += row;
            𝔼xxᵀ_pre += row*row';
        end
    end

    # ess = 𝔼x, 𝔼xxᵀ, ismissing, missingRows, n_obs
    # estepFn = (em_model, data) -> estep(em_model, data, EXsum, EXXsum, ismissing, missingRows, n_obs)

    # initialize
    em_model = start_em(observed; kwargs...)
    em_model_prev = EmMVNModel(zeros(n_man, n_man), zeros(n_man), false)
    iter = 1
    done = false
    𝔼x = zeros(n_man)
    𝔼xxᵀ = zeros(n_man, n_man)

    while !done

        em_mvn_Estep!(𝔼x, 𝔼xxᵀ, em_model, observed, 𝔼x_pre, 𝔼xxᵀ_pre)
        em_mvn_Mstep!(em_model, n_obs, 𝔼x, 𝔼xxᵀ)

        if iter > max_iter_em
            done = true
            @warn "EM Algorithm for MVN missing data did not converge. Likelihood for FIML is not interpretable.
            Maybe try passing different starting values via 'start_em = ...' "
        elseif iter > 1
            # done = isapprox(ll, ll_prev; rtol = rtol)
            done = isapprox(em_model_prev.μ, em_model.μ; rtol = rtol_em) & isapprox(em_model_prev.Σ, em_model.Σ; rtol = rtol_em)
        end

        # print("$iter \n")
        iter = iter + 1
        em_model_prev.μ, em_model_prev.Σ = em_model.μ, em_model.Σ

    end

    # update EM Mode in observed
    observed.em_model.Σ .= em_model.Σ
    observed.em_model.μ .= em_model.μ
    observed.em_model.fitted = true

    return nothing

end

# E and M step ------------------------------------------------------------------------------

function em_mvn_Estep!(𝔼x, 𝔼xxᵀ, em_model, observed, 𝔼x_pre, 𝔼xxᵀ_pre)

    𝔼x .= 0.0
    𝔼xxᵀ .= 0.0

    𝔼xᵢ = copy(𝔼x)
    𝔼xxᵀᵢ = copy(𝔼xxᵀ)

    μ = em_model.μ
    Σ = em_model.Σ

    # Compute the expected sufficient statistics
    for i in 2:length(observed.pattern_n_obs)

        # observed and unobserved vars
        u = observed.patterns_not[i]
        o = observed.patterns[i]

        # precompute for pattern
        V = Σ[u, u] - Σ[u, o] * (Σ[o, o]\Σ[o, u])

        # loop trough data
        for row in observed.rows[i]
            m = μ[u] + Σ[u, o] * ( Σ[o, o] \ (observed.data_rowwise[row]-μ[o]) )

            𝔼xᵢ[u] = m
            𝔼xᵢ[o] = observed.data_rowwise[row]
            𝔼xxᵀᵢ[u, u] = 𝔼xᵢ[u] * 𝔼xᵢ[u]' + V
            𝔼xxᵀᵢ[o, o] = 𝔼xᵢ[o] * 𝔼xᵢ[o]'
            𝔼xxᵀᵢ[o, u] = 𝔼xᵢ[o] * 𝔼xᵢ[u]'
            𝔼xxᵀᵢ[u, o] = 𝔼xᵢ[u] * 𝔼xᵢ[o]'

            𝔼x .+= 𝔼xᵢ
            𝔼xxᵀ .+= 𝔼xxᵀᵢ
        end

    end

    𝔼x .+= 𝔼x_pre
    𝔼xxᵀ .+= 𝔼xxᵀ_pre

end

function em_mvn_Mstep!(em_model, n_obs, 𝔼x, 𝔼xxᵀ)

    em_model.μ = 𝔼x/n_obs;
    Σ = Symmetric(𝔼xxᵀ/n_obs - em_model.μ*em_model.μ')

    # ridge Σ
    # while !isposdef(Σ)
    #     Σ += 0.5I
    # end

    em_model.Σ = Σ

    # diagonalization
    #if !isposdef(Σ)
    #    print("Matrix not positive definite")
    #    em_model.Σ .= 0
    #    em_model.Σ[diagind(em_model.Σ)] .= diag(Σ)
    #else
        # em_model.Σ = Σ
    #end

    return nothing
end

# generate starting values --------------------------------------------------------------

# use μ and Σ of full cases
function start_em_observed(observed::SemObsMissing; kwargs...)

    if (length(observed.patterns[1]) == observed.n_man) & (observed.pattern_n_obs[1] > 1)
        μ = copy(observed.obs_mean[1])
        Σ = copy(Symmetric(observed.obs_cov[1]))
        if !isposdef(Σ)
            Σ = Matrix(Diagonal(Σ))
        end
        return EmMVNModel(Σ, μ, false)
    else
        return start_em_simple(observed, kwargs...)
    end

end

# use μ = O and Σ = I
function start_em_simple(observed::SemObsMissing; kwargs...)
    n_man = Int(observed.n_man)
    μ = zeros(n_man)
    Σ = rand(n_man, n_man); Σ = Σ*Σ'
    # Σ = Matrix(1.0I, n_man, n_man)
    return EmMVNModel(Σ, μ, false)
end

# set to passed values
function start_em_set(observed::SemObsMissing; model_em, kwargs...)
    return em_model
end


text:
"""
gives the ordering for nedelec2 (divergence)
"""
function relorientation(face::SArray{Tuple{3},T,1,3}, tet::SArray{Tuple{4},T,1,4}) where {T}
    # v = setdiff(tet,face)
    # length(v) == 1 || return 0

    # a = something(findfirst(isequal(v[1]), tet),0)

    a = findfirst(x -> !(x in face), tet)
    a == nothing && return 0

    @assert 1 <= a <= 4
    q = findnext(x ->!(x in face), tet, a+1)
    q == nothing || return 0

    v = tet[a]

    Face = deleteat(tet,a)
    # @show a tet face Face indexin(face,Face)
    p = parity(Vector{Int}(indexin(face,Face)))
    return a * (-1)^a * (-1)^p

    # face = [v[1],face[1],face[2],face[3]]
    # w = sortperm(face)
    # b = parity(w)
    # return -a*(-1)^b
end

"""
gives the ordering for nedelec1 (curl)
if a tetrahedron is given by simplex(a,b,c,d), the order is:
    a-b
    a-c
    a-d
    b-c
    b-d
    c-d
"""

function relorientation(edge::SArray{Tuple{2},Int64,1,2}, tet::SArray{Tuple{4},Int64,1,4})
    # v = setdiff(tet, edge)
    # length(v) == 2 || return 0

    count(x -> !(x in edge), tet) != 2 && return 0

    w1 = findfirst(isequal(edge[1]), tet)
    w2 = findfirst(isequal(edge[2]), tet)
    t = tetschoice(w1,w2)
    s = sign(w2-w1)

    return s*t
end

function tetschoice(w1::Int64,w2::Int64)
    a = min(w1,w2)
    b = max(w1,w2)
    if a == 1
        return b-1
    elseif a == 2
        return b+1
    elseif a == 3
        return b+2
    end
end

function relorientation(node::SArray{Tuple{1},Int64,1,1}, face::SArray{Tuple{N},Int64,1,N} where {N})
    idx = something(findfirst(isequal(node[1]), face), 0)
    sgn = (-1)^(idx-1)
    return sgn * idx
end

function relorientation(face, simplex)

    v = setdiff(simplex, face)
    length(v) == 1 || return 0

    # find the position of the missing vertex
    v = v[1]
    #i = Base.findfirst(simplex, v)
    i = something(findfirst(isequal(v), simplex),0)
    s = (-1)^(i-1)

    # remove that vertex from the simplex
    face2 = Array{Int}(undef,length(simplex)-1)
    for j in 1 : i-1
        face2[j] = simplex[j]
    end
    for j in i : length(simplex)-1
        face2[j] = simplex[j+1]
    end

    # get the permutation that maps face to face2
    #p = indexin(face, face2)
    #p = [ findfirst(face2,v) for v in face ]
    p = [something(findfirst(isequal(v),face2),0) for v in face]

    return s * levicivita(p) * i
end


text:
valid = 0

open("day02/day02.in") do io
    global valid
    for line in readlines(io)
        range_all, target_char, test_string = split(line, " ")
        target_char = target_char[1]
        range_min = parse(Int32, split(range_all, "-")[1])
        range_max = parse(Int32, split(range_all, "-")[2])

        char_count = 0
        for ch in test_string
            if ch == target_char
                char_count += 1
            end
        end
        if char_count >= range_min && char_count <= range_max
            valid += 1
        end
    end
end

println(valid)

max_stars_repo_name: detrin/advent-of-code

max_issues_repo_name: detrin/advent-of-code

max_forks_repo_name: detrin/advent-of-code


text:
# Each file of the form "bench_$(name).jl" in this directory is `include`d and
# its last statement is assumed to be a `BenchmarkGroup`.  This group is added
# to the top-level group `SUITE` with the `$name` extracted from the file name.

using BenchmarkTools
const SUITE = BenchmarkGroup()
for file in sort(readdir(@__DIR__))
    if startswith(file, "bench_") && endswith(file, ".jl")
        SUITE[chop(file, head = length("bench_"), tail = length(".jl"))] =
            include(file)
    end
end


text:
using DrWatson
@quickactivate "StatReth"

# %%
using StatisticalRethinking
using CSV
# using DataFrames
using StatsBase
using Distributions
using Plots
using StatsPlots

# %% 4.7, 4.11
d = CSV.read(datadir("exp_raw/Howell_1.csv"), copycols = true)
precis(d)

d2 = filter(row -> row.age >= 18, d)
d2 = d[d.age .>= 18, :]

# %%
d.height
d.weight

f = fit(Histogram, d.weight, nbins = 12)

histogram(d.weight)

# %%
f = fit(Histogram, d.weight, nbins = 12)
# f = fit(Histogram, d.weight, range(0, 75, length = 13))
plot(f)
hline!(range(0, maximum(f.weights), length = 9))
savefig(plotsdir("histogram"))

# %%

f.weights

const BARS = collect("▁▂▃▄▅▆▇█")

function unicode_histogram(data)
  f = fit(Histogram, data, nbins = 12)  # nbins: more like a guideline than a rule, really
  # scale weights between 1 and 8 (length(BARS)) to fit the indices in BARS
  # eps is needed so indices are in the interval [0, 8) instead of [0, 8] which could
  # result in indices 0:8 which breaks things
  scaled = f.weights .* (length(BARS) / maximum(f.weights) - eps())
  indices = floor.(Int, scaled) .+ 1
  return join((BARS[i] for i in indices))
end

for col in eachcol(d)
    println(unicode_histogram(col))
end

fit(histogram)

# %%
 # ▁▂▃▄▅▆▇█

 Example usage:
 ```
 d = CSV.read(datadir("exp_raw/Howell_1.csv"), copycols = true)
 for col in eachcol(d)
     println(unicode_histogram(col))
 end
 ```
 gives
 ```
 ▁▁▁▂▂▂▂▂▂██▆▁
 ▁▃▄▄▃▂▃▆██▅▃▁
 █▆▆▆▆▃▃▁▁
 █▁▁▁▁▁▁▁▁▁█
 ```

 I'm not sure how to incorporate this with `precis` since it would probably turn the array type in to a `Union`. Unicode and font coverage is also necessary but that seems to be pretty good in the Julia community (heavy Unicode symbol usage, strings are UTF-8...).


text:
cd(@__DIR__); include("setups/grid23x22.jl")
gr(dpi = 200)

## build SGWT frame
frame = sgwt_frame(W; nf = 6)
SGWT = reshape(frame, (N, :))
SGWT_dual = (SGWT * SGWT') \ SGWT
f = digit_img[:]

important_idx = sortperm((SGWT' * f).^2; rev = true)
plot(layout = Plots.grid(10,10), size = (2300, 2800))
for i = 1:100
    grid_vector_plot!(important_idx[i], i, SGWT)
    if i == 90
        plot_square!(Nx, Ny; subplot = i)
    end
end
plt = current()
savefig(plt, "../figs/Grid23x22_fdigit3_sgwt_top100.png")


text:
using Q._k

"""
    K_new(x)

Construct a new kdb+ object or create a new reference to an existing one. The
caller is reponsible for calling `r0` once the object is no longer needed.
"""
K_new

"""
    K_Ref(x)

Construct a Julia GC managed object that owns a kdb+ reference.
"""
mutable struct K_Ref
    x::K_
    K_Ref(p::K_) = (x = new(p); finalizer(x, (x)->r0(x.x));x)
end

struct TI
    number::C_
    letter::Char
    name::String
    c_type::Type
    jl_type::Type
    class::Symbol
end

const TYPE_INFO = [
    # num ltr name c_type jl_type super
    TI(1,  'b', "boolean",   B_, Bool,    :_Bool),

    TI(2,  'g', "guid",      U_, UInt128, :_Unsigned),
    TI(4,  'x', "byte",      G_, UInt8,   :_Unsigned),

    TI(5,  'h', "short",     H_, Int16,   :_Signed),
    TI(6,  'i', "int",       I_, Int32,   :_Signed),
    TI(7,  'j', "long",      J_, Int64,   :_Signed),

    TI(8,  'e', "real",      E_, Float32, :_Float),
    TI(9,  'f', "float",     F_, Float64, :_Float),

    TI(10, 'c', "char",      G_, UInt8,   :_Text),
    TI(11, 's', "symbol",    S_, Symbol,  :_Text),

    TI(12, 'p', "timestamp", J_, TimeStamp,   :_Time),
    TI(13, 'm', "month",     I_, Month,   :_Time),
    TI(14, 'd', "date",      I_, Date,    :_Time),
    TI(15, 'z', "datetime",  F_, DateTimeF, :_Time),
    TI(16, 'n', "timespan",  J_, TimeSpan,   :_Period),
    TI(17, 'u', "minute",    I_, Minute,   :_Time),
    TI(18, 'v', "second",    I_, Second,   :_Time),
    TI(19, 't', "time",      I_, Time,   :_Time),
]
const TYPE_CLASSES = unique(t.class for t in TYPE_INFO)
const C_TYPE = merge(
    Dict(KK=>K_, EE=>S_, XT=>K_, XD=>K_, (-EE)=>S_,  # XXX: do we need both ±EE?
         100=>K_, 101=>I_, 102=>I_, 103=>I_, 104=>K_, 105=>K_,
         106=>K_, 107=>K_, 108=>K_, 109=>K_, 110=>K_, 111=>K_, 112=>V_),
    Dict(t.number=>t.c_type for t in TYPE_INFO))

const K_TYPE = Dict(Bool=>KB, UInt128=>UU,
                    UInt8=>KG, Int16=>KH, Int32=>KI, Int64=>KJ,
                    Float32=>KE, Float64=>KF,
                    Char=>KC, Symbol=>KS, Cstring=>KS,
                    Date=>KD,)
# returns type, offset and size
function cinfo(x::K_)
    h = unsafe_load(x)
    t = h.t
    if t < 0   # scalar
        return C_TYPE[-t], (t == -UU ? 16 : 8), ()
    elseif t <= KT
        return C_TYPE[t], 16, (xn(x), )
    elseif t == XT
        return K_, 8, ()
    elseif t == XD
        return K_, 16, (2, )
    elseif t < 77
        return I_, 16, (xn(x), )
    elseif t < XT
        return J_, 16, (xn(x), )
    elseif t == 100  # λ
        return K_, 16, (xn(x), )
    elseif t < 104
        return I_, 8, ()
    elseif t < 106  # projection, composition
        return K_, 16, (xn(x), )
    elseif t < 112  # f', f/, f\, ...
        return K_, 8, ()
    else
        return Ptr{V_}, 16, ()
    end
end
## New reference
K_new(x::K_) = r1(x)
## Conversion of simple types
const _none = ktj(101, 0)
K_new(::Void) = r1(_none)
K_new(x::Bool) = kb(x)
K_new(x::UInt128) = ku(x)
K_new(x::UInt8) = kg(x)
K_new(x::Int16) = kh(x)
K_new(x::Int32) = ki(x)
K_new(x::Integer) = kj(x)
K_new(x::Float32) = ke(x)
K_new(x::Real) = kf(x)
K_new(x::Symbol) = ks(String(x))
K_new(x::Date) = kd(DATE_SHIFT + Dates.value(x))
K_new(x::DateTime) = ktj(-KP, 10^6*Dates.toms(x - DateTime(2000)))
K_new(x::Dates.TimePeriod) = ktj(-KN, Dates.tons(x))
K_new(x::Char) = kc(I_(x))
K_new(x::String) = kp(x)
## Vector conversions
function K_new(a::Vector{T}) where {T<:Number}
    t = K_TYPE[T]
    C = C_TYPE[t]
    n = length(a)
    x = ktn(t, n)
    unsafe_copy!(Ptr{C}(x+16), pointer(a), n)
    return x
end
function K_new(a::AbstractVector{T}) where {T<:Number}
    t = K_TYPE[T]
    C = C_TYPE[t]
    n = length(a)
    x = ktn(t, n)
    for i in 1:n
        unsafe_store!(Ptr{C}(x+16), C(a[i]), i)
    end
    return x
end
function K_new(a::AbstractVector{Symbol})
    t = KS
    n = length(a)
    x = ktn(t, n)
    for i in 1:n
        si = ss(a[i])
        unsafe_store!(Ptr{S_}(x+16), si, i)
    end
    return x
end
function K_new(a::AbstractVector{T}) where {T<:AbstractString}
    n = length(a)
    x = ktn(0, n)
    for i in 1:n
        si = kp(a[i])
        unsafe_store!(Ptr{K_}(x+16), si, i)
    end
    return x
end
function K_new(a::Union{Tuple,Vector{Any}})
    x = ktn(0, 0)
    r = Ref{K_}(x)
    for el in a
        jk(r, el isa K_ ? r1(el) : K_new(el))
    end
    r.x
end
# TODO: Add a fast specialization for the concrete Matrix type.
function K_new(m::AbstractMatrix)
    T = eltype(m)
    t = K_TYPE[T]
    C = C_TYPE[t]
    nrows, ncols = size(m)
    x = ktn(0, ncols)
    for j in 1:ncols
        pcol = kK(x)[j] = ktn(t, nrows)
        for i in 1:nrows
            kX(C, pcol)[i] = _cast(C, m[i,j])
        end
    end
    x
end

function K_new(df::AbstractDataFrame)
    x = K_new(names(df))
    y = ktn(0, xn(x))
    for (i, col) in enumerate(DataFrames.columns(df))
        kK(y)[i] = K_new(col)
    end
    xT(xD(x, y))
end

function asarray(x::K_, own::Bool=true)
    T, o, s = cinfo(x)
    a = unsafe_wrap(Array, Ptr{T}(x + o), s)
    if own
        finalizer(a, b->r0(K_(pointer(b)-o)))
    end
    a
end
# The _none pointer guard - make sure _none is cleaned up eventually.
const _none_array = asarray(_none)
# TODO: Consider using Array(x) instead of asarray(x).
# Base.convert(::Type{Array}, x::K_) = asarray(x)


text:
include(joinpath("..", "..", ".slbootstrap.jl"))
using Test

### Your tests here
@test 1 == 1

max_stars_repo_name: tetsugps/julia

max_issues_repo_name: tetsugps/julia

max_forks_repo_name: tetsugps/julia


text:
using Test
using Clustering

@testset "silhouettes()" begin

local D = [0 1 2 3
           1 0 1 2
           2 1 0 1
           3 2 1 0]

@assert size(D) == (4, 4)

local a = [1, 1, 2, 2]
local c = [2, 2]

@testset "Input checks" begin
    @test_skip silhouettes(a, [1, 1, 2], D) # should throw because cluster counts are inconsistent
    @test_throws ArgumentError silhouettes(a, [3, 2], D)
    @test_throws ArgumentError silhouettes([1, 1, 3, 2], [2, 2], D)
    @test_throws DimensionMismatch silhouettes([1, 1, 2, 2, 2], [2, 3], D)
end

@test @inferred(silhouettes(a, c, D)) ≈ [1.5/2.5, 0.5/1.5, 0.5/1.5, 1.5/2.5]
@test @inferred(silhouettes(a, c, convert(Matrix{Float32}, D))) isa AbstractVector{Float32}
@test silhouettes(a, D) == silhouettes(a, c, D) # c is optional

a = [1, 2, 1, 2]
c = [2, 2]

@test silhouettes(a, c, D) ≈ [0.0, -0.5, -0.5, 0.0]

a = [1, 1, 1, 2]
c = [3, 1]

@test silhouettes(a, c, D) ≈ [0.5, 0.5, -1/3, 0.0]

@testset "silhouettes() handle zero cluster distances correctly" begin
    a = [fill(1, 5); fill(2, 5)]
    d = fill(0, (10, 10))

    @test silhouettes(a, d) == fill(0.0, 10)

    d = fill(1, (10, 10))
    d[1, 2] = d[2, 1] = 5

    @test silhouettes(a, d) == [[-0.5, -0.5]; fill(0.0, 8)]
end

@testset "silhouette() throws an error when degenerated clustering is given" begin
    a = fill(1, 10)
    d = fill(1, (10, 10))
    for i in 1:10; d[i, i] = 0; end

    @test_throws ArgumentError silhouettes(a, d)
end

end


text:
###
# Test using simulated data
#
# For increasing K, simulate datasets and performance
# of various methods SAA, Alpha LOO, and Alpha OR
###
using Distributions, Random, DelimitedFiles
include("../src/JS_SAA_main.jl")

#supp_full, ps_full are d x K matrices with true info per problem
#adds new subproblems in order appear in files
#s is the service level, N is the average amount of data per problem
function convInKtest(numRuns, K_grid, supp_full, ps_full, outPath, N_grid, s;
						usePoisson=true, seed=8675309)
	Random.seed!(seed)
	Kmax = maximum(K_grid)
	@assert Kmax <= size(supp_full, 2) "K_grid exceeds available subproblems"
	@assert size(supp_full) == size(ps_full) "supp_full and ps_full have incompatible dimensions"

	d = size(supp_full, 1)

	#For safety, trim inputs to size Kmax
	supp_full = view(supp_full, 1:d, 1:Kmax)
	ps_full = view(ps_full, 1:d, 1:Kmax)

	p0 = ones(d)/d
	alpha_grid = range(0, stop=180, length=120)

	#set up output file
	f = open("$(outPath).csv", "w")
	writedlm(f, ["Run" "K" "d" "N" "Method" "TruePerf" "time" "alpha"], ',')

	#generate all Kmax subproblems upfront and store in memory
	cs_full, xs_full = JS.genNewsvendorsDiffSupp(supp_full, s, Kmax)
	lam_full = ones(Kmax)

	for (iRun, N) = Iterators.product(1:numRuns, N_grid)
		#simulate data for the run
		Nhats_full = usePoisson ? rand(Poisson(N), Kmax) : N * ones(Int, Kmax)

		#for now, cludge the Nhatk == 0 issue
		Nhats_full[Nhats_full .== 0] .= 1
		mhats_full = JS.sim_path(ps_full, Nhats_full)

		for K in K_grid
			#Take views on evrything for simplicity
			lams = view(lam_full, 1:K)
			supp = view(supp_full, 1:d, 1:K)
			ps = view(ps_full, 1:d, 1:K)
			Nhats = view(Nhats_full, 1:K)
			mhats = view(mhats_full, 1:d, 1:K)
			cs = view(cs_full, 1:d, 1:K)
			xs = view(xs_full, 1:K)

			#for data-driven shrinkage anchor
			phat_avg = vec(mean(mhats ./ Nhats', dims=2))

			#Compute the full-info value once for reference
			if iRun == 1
				t =
				  @elapsed full_info = JS.zstar(xs, cs, ps, lams)
				writedlm(f, [1 K d N "FullInfo" full_info t 0.], ',')
			end

			#SAA
			t =
			  @elapsed perf_SAA = JS.zbar(xs, cs, p0, 0., mhats, ps, lams)
			writedlm(f, [iRun K d N "SAA" perf_SAA t 0.0], ',')

			#Gen the Oracle cost with 1/d anchor
			t =
			  @elapsed alphaOR, min_indx, or_alpha_curve = JS.oracle_alpha(xs, cs, mhats, ps, lams, p0, alpha_grid)
			writedlm(f, [iRun K d N "Oracle" or_alpha_curve[min_indx] t alphaOR], ',')

			#Gen the LOO cost with 1/d anchor
			t =
			  @elapsed alphaLOO, min_indx, looUnsc_curve = JS.loo_alpha(xs, cs, mhats, p0, alpha_grid)
			writedlm(f, [iRun K d N "LOO_unif" or_alpha_curve[min_indx] t alphaLOO], ',')

			##MSE version of alpha
			t =
			  @elapsed alphaMSE, min_indx = JS.mse_estimates(mhats, supp, p0, alpha_grid)
			writedlm(f, [iRun K d N "MSE" or_alpha_curve[min_indx] t alphaMSE], ',')

			#Gen the Oracle cost with GM Anchor
			t =
			  @elapsed alphaOR_GM, min_indx, or_alpha_curve_GM = JS.oracle_alpha(xs, cs, mhats, ps, lams, phat_avg, alpha_grid)
			writedlm(f, [iRun K d N "OraclePhat" or_alpha_curve_GM[min_indx] t alphaOR_GM], ',')

			#Gen the LOO cost with the GM Anchor
			t =
			  @elapsed alphaLOO, min_indx, looUnsc_curve = JS.loo_alpha(xs, cs, mhats, phat_avg, alpha_grid)
			writedlm(f, [iRun K d N "LOO_avg" or_alpha_curve_GM[min_indx] t alphaLOO], ',')

			##MSE version of alpha with GM
			t =
			  @elapsed alphaMSe, min_indx = JS.mse_estimates(mhats, supp, phat_avg, alpha_grid)
			writedlm(f, [iRun K d N "MSE_GM" or_alpha_curve_GM[min_indx] t alphaMSE], ',')

		end  #end K Loop
		flush(f)
	end #end run loop
	close(f)
	"$(outPath).csv"
end #end function

max_stars_repo_name: CausalML/DataPooling

max_issues_repo_name: CausalML/DataPooling

max_forks_repo_name: CausalML/DataPooling


text:
using DPC, CairoMakie
include("utils.jl")

files = ["bio_schema.jl", "inhibition.jl", "multineuron_and.jl", "multineuron_or.jl", "place_cells.jl", "probabilistic.jl", "motifs.jl"]
for file in files
    let
        println("Running $(file) ...")
        include(file)
    end
end

println("DONE.")


text:
# Wrapper for level-3 BLIS routines.
#

# Level3 API of common forms are put together.
#
macro blis_group_level3_hemm(funcname)
    return quote
        @blis_ccall_group($funcname,
                          Cvoid,
                          sidea,   BliSide,
                          uploa,   BliUpLo,
                          conja,   BliConj,
                          transb,  BliTrans,
                          m,       BliDim,
                          n,       BliDim,
                          α,       Ptr{xType},
                          a,       Ptr{xType}, rsa, BliInc, csa, BliInc,
                          b,       Ptr{xType}, rsb, BliInc, csb, BliInc,
                          β,       Ptr{xType},
                          c,       Ptr{xType}, rsc, BliInc, csc, BliInc)
    end
end

@blis_group_level3_hemm hemm
@blis_group_level3_hemm symm

macro blis_group_level3_tr(funcname)
    return quote
        @blis_ccall_group($funcname,
                          Cvoid,
                          sidea,   BliSide,
                          uploa,   BliUpLo,
                          transa,  BliTrans,
                          diaga,   BliDiag,
                          m,       BliDim,
                          n,       BliDim,
                          α,       Ptr{xType},
                          a,       Ptr{xType}, rsa, BliInc, csa, BliInc,
                          b,       Ptr{xType}, rsb, BliInc, csb, BliInc)
    end
end

@blis_group_level3_tr trmm
@blis_group_level3_tr trsm

# Forms with only 1 member are defined directly.
#
@blis_ccall_group(gemm,
                  Cvoid,
                  transa,  BliTrans,
                  transb,  BliTrans,
                  m,       BliDim,
                  n,       BliDim,
                  k,       BliDim,
                  α,       Ptr{xType},
                  a,       Ptr{xType}, rsa, BliInc, csa, BliInc,
                  b,       Ptr{xType}, rsb, BliInc, csb, BliInc,
                  β,       Ptr{xType},
                  c,       Ptr{xType}, rsc, BliInc, csc, BliInc)
@blis_ccall_group(gemmt,
                  Cvoid,
                  uplo,    BliUpLo,
                  transa,  BliTrans,
                  transb,  BliTrans,
                  m,       BliDim,
                  k,       BliDim,
                  α,       Ptr{xType},
                  a,       Ptr{xType}, rsa, BliInc, csa, BliInc,
                  b,       Ptr{xType}, rsb, BliInc, csb, BliInc,
                  β,       Ptr{xType},
                  c,       Ptr{xType}, rsc, BliInc, csc, BliInc)
@blis_ccall_group(herk,
                  Cvoid,
                  uploc,   BliUpLo,
                  transa,  BliTrans,
                  m,       BliDim,
                  k,       BliDim,
                  α,       Ptr{rType},
                  a,       Ptr{xType}, rsa, BliInc, csa, BliInc,
                  β,       Ptr{rType},
                  c,       Ptr{xType}, rsc, BliInc, csc, BliInc)
@blis_ccall_group(her2k,
                  Cvoid,
                  uploc,   BliUpLo,
                  transab, BliTrans,
                  m,       BliDim,
                  k,       BliDim,
                  α,       Ptr{xType},
                  a,       Ptr{xType}, rsa, BliInc, csa, BliInc,
                  b,       Ptr{xType}, rsb, BliInc, csb, BliInc,
                  β,       Ptr{rType},
                  c,       Ptr{xType}, rsc, BliInc, csc, BliInc)
@blis_ccall_group(syrk,
                  Cvoid,
                  uploc,   BliUpLo,
                  transa,  BliTrans,
                  m,       BliDim,
                  k,       BliDim,
                  α,       Ptr{xType},
                  a,       Ptr{xType}, rsa, BliInc, csa, BliInc,
                  β,       Ptr{xType},
                  c,       Ptr{xType}, rsc, BliInc, csc, BliInc)
@blis_ccall_group(syr2k,
                  Cvoid,
                  uploc,   BliUpLo,
                  transab, BliTrans,
                  m,       BliDim,
                  k,       BliDim,
                  α,       Ptr{xType},
                  a,       Ptr{xType}, rsa, BliInc, csa, BliInc,
                  b,       Ptr{xType}, rsb, BliInc, csb, BliInc,
                  β,       Ptr{xType},
                  c,       Ptr{xType}, rsc, BliInc, csc, BliInc)
@blis_ccall_group(trmm3,
                  Cvoid,
                  sidea,   BliSide,
                  uploa,   BliUpLo,
                  transa,  BliTrans,
                  diaga,   BliDiag,
                  transb,  BliTrans,
                  m,       BliDim,
                  n,       BliDim,
                  α,       Ptr{xType},
                  a,       Ptr{xType}, rsa, BliInc, csa, BliInc,
                  b,       Ptr{xType}, rsb, BliInc, csb, BliInc,
                  β,       Ptr{xType},
                  c,       Ptr{xType}, rsc, BliInc, csc, BliInc)


text:
function givens(x::Array{N,1}) where {N <: Number}

    if length(x) != 2
        error("'x' should have length 2")
    end

    z = x[1]
    y = x[2]

    if y == zero(N)
        c = one(N)
        s = zero(N)
    elseif abs(y'y) >= abs(z'z)
        tau = -z/y
        s = one(N)/(sqrt(one(N)+tau'tau))
        c = s*tau
    else
        tau = -y/z
        c = one(N)/(sqrt(one(N)+tau'tau))
        s = c*tau
    end

    g = [c s'; -s c']

    return g

end

function givens!(g::Array{N,2},x::AbstractArray{N,1}) where {N <: Number}

    if length(x) != 2
        error("'x' should have length 2")
    end

    z = x[1]
    y = x[2]

    if y == zero(N)
        c = one(N)
        s = zero(N)
    elseif abs(y'y) >= abs(z'z)
        tau = -z/y
        s = one(N)/(sqrt(one(N)+tau'tau))
        c = s*tau
    else
        tau = -y/z
        c = one(N)/(sqrt(one(N)+tau'tau))
        s = c*tau
    end

    g[1] = c
    g[2] = -s
    g[3] = s'
    g[4] = c'

end


text:
using OpenQuantumSystems
import QuantumOpticsBase, LinearAlgebra, OrdinaryDiffEq, QuadGK, DelayDiffEq

function QME_sI_ansatz_test(
    W0::T,
    tspan::Array,
    agg::Aggregate;
    reltol::AbstractFloat = 1.0e-12,
    abstol::AbstractFloat = 1.0e-12,
    int_reltol::AbstractFloat = 1.0e-4,
    int_abstol::AbstractFloat = 1.0e-4,
    alg::Any = DelayDiffEq.MethodOfSteps(DelayDiffEq.Vern6()),
    fout::Union{Function,Nothing} = nothing,
    kwargs...,
) where {B<:Basis,T<:Operator{B,B}}
    history_fun(p, t) = T(rho0.basis_l, rho0.basis_r, zeros(ComplexF64, size(rho0.data)))
    rho0 = trace_bath(W0, agg.core, agg.tools; vib_basis=agg.operators.vib_basis)
    W0_bath = get_rho_bath(W0, agg.core, agg.tools; vib_basis=agg.operators.vib_basis)
    p = (agg.core, agg.tools, agg.operators, W0, eltype(W0))

    tmp1 = copy(W0.data)
    tmp2 = copy(W0.data)
    dmaster_(t, rho, drho, history_fun, p) = dQME_sI_ansatz_test(
        t,
        rho,
        drho,
        history_fun,
        tmp1,
        tmp2,
        p,
        int_reltol,
        int_abstol,
    )
    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
    x0 = rho0.data
    state = T(rho0.basis_l, rho0.basis_r, rho0.data)
    dstate = T(rho0.basis_l, rho0.basis_r, rho0.data)
    OpenQuantumSystems.integrate_delayed(
        tspan_,
        dmaster_,
        history_fun,
        x0,
        state,
        dstate,
        fout;
        p = p,
        reltol = reltol,
        abstol = abstol,
        alg = alg,
        kwargs...,
    )
end

function dQME_sI_ansatz_test(
    t::AbstractFloat,
    rho::T,
    drho::T,
    history_fun,
    tmp1::Array,
    tmp2::Array,
    p,
    int_reltol::AbstractFloat,
    int_abstol::AbstractFloat,
) where {B<:Basis,T<:Operator{B,B}}
    aggCore, aggTools, aggOperators, W0, elementtype = p

    Ham_II_t = getInteractionHamIPicture(aggOperators.Ham_0, aggOperators.Ham_I, t)
    K = Ham_II_t.data * W0.data - W0.data * Ham_II_t.data
    K_traced = trace_bath(K, aggCore, aggTools; vib_basis=aggOperators.vib_basis)

    kernel_integrated_traced, err = QuadGK.quadgk(
        s -> kernel_sI_ansatz_test(t, s, history_fun, p, tmp1, tmp2, Ham_II_t),
        0,
        t,
        rtol = int_reltol,
        atol = int_abstol,
    )
    drho.data[:, :] = -elementtype(im) * K_traced - kernel_integrated_traced

    return drho
end

function kernel_sI_ansatz_test(t, s, h, p, tmp1, tmp2, Ham_II_t)
    aggCore, aggTools, aggOperators, W0, _ = p

    rho = h(p, s)

    if (typeof(rho) <: Operator)
        rho = rho.data
    end

    Ham_0 = aggOperators.Ham_0
    Ham_I = aggOperators.Ham_I
    Ham = aggOperators.Ham
    Ham_II_s = getInteractionHamIPicture(Ham_0, Ham_I, s)

    U_0_s = evolutionOperator(Ham, s)
    W_s = U_0_s * W0 * U_0_s'
    U_0_op = evolutionOperator(Ham_0, s)
    W_int_s = U_0_op' * W_s * U_0_op

    W_bath_s = get_rho_bath(W_int_s, aggCore, aggTools; vib_basis=aggOperators.vib_basis)
    rho_s = trace_bath(W_int_s, aggCore, aggTools; vib_basis=aggOperators.vib_basis)
    tmp1[:, :] = ad(rho_s.data, W_bath_s.data, aggCore, aggTools)

    tmp2[:, :] = Ham_II_s.data * tmp1 - tmp1 * Ham_II_s.data
    tmp1[:, :] = Ham_II_t.data * tmp2 - tmp2 * Ham_II_t.data

    return trace_bath(tmp1, aggCore, aggTools; vib_basis=aggOperators.vib_basis)
end

function QME_sI_ansatz_const_int(
    W0::T,
    tspan::Array,
    agg::Aggregate;
    reltol::AbstractFloat = 1.0e-12,
    abstol::AbstractFloat = 1.0e-12,
    int_reltol::AbstractFloat = 1.0e-4,
    int_abstol::AbstractFloat = 1.0e-4,
    alg::Any = DelayDiffEq.MethodOfSteps(DelayDiffEq.Vern6()),
    fout::Union{Function,Nothing} = nothing,
    kwargs...,
) where {B<:Basis,T<:Operator{B,B}}
    history_fun(p, t) = T(rho0.basis_l, rho0.basis_r, zeros(ComplexF64, size(rho0.data)))
    rho0 = trace_bath(W0, agg.core, agg.tools; vib_basis=agg.operators.vib_basis)
    W0_bath = get_rho_bath(W0, agg.core, agg.tools; vib_basis=agg.operators.vib_basis)
    p = (agg.core, agg.tools, agg.operators, W0, W0_bath, eltype(W0))

    tmp1 = copy(W0.data)
    tmp2 = copy(W0.data)
    dmaster_(t, rho, drho, history_fun, p) = dQME_sI_ansatz_const_int(
        t,
        rho,
        drho,
        history_fun,
        tmp1,
        tmp2,
        p,
        int_reltol,
        int_abstol,
    )
    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
    x0 = rho0.data
    state = T(rho0.basis_l, rho0.basis_r, rho0.data)
    dstate = T(rho0.basis_l, rho0.basis_r, rho0.data)
    OpenQuantumSystems.integrate_delayed(
        tspan_,
        dmaster_,
        history_fun,
        x0,
        state,
        dstate,
        fout;
        p = p,
        reltol = reltol,
        abstol = abstol,
        alg = alg,
        kwargs...,
    )
end

function dQME_sI_ansatz_const_int(
    t::AbstractFloat,
    rho::T,
    drho::T,
    history_fun,
    tmp1::Array,
    tmp2::Array,
    p,
    int_reltol::AbstractFloat,
    int_abstol::AbstractFloat,
) where {B<:Basis,T<:Operator{B,B}}
    aggCore, aggTools, aggOperators, W0, _, elementtype = p

    Ham_II_t = getInteractionHamIPicture(aggOperators.Ham_0, aggOperators.Ham_I, t)
    K = Ham_II_t.data * W0.data - W0.data * Ham_II_t.data
    K_traced = trace_bath(K, aggCore, aggTools; vib_basis=aggOperators.vib_basis)

    kernel_integrated_traced, err = QuadGK.quadgk(
        s -> kernel_sI_ansatz_const_int(t, s, history_fun, p, tmp1, tmp2, Ham_II_t),
        0,
        t,
        rtol = int_reltol,
        atol = int_abstol,
    )
    drho.data[:, :] = -elementtype(im) * K_traced - kernel_integrated_traced

    return drho
end

function kernel_sI_ansatz_const_int(t, s, h, p, tmp1, tmp2, Ham_II_t)
    aggCore, aggTools, aggOperators, W0, W0_bath, _ = p

    rho = h(p, s)

    if (typeof(rho) <: Operator)
        rho = rho.data
    end

    Ham_0 = aggOperators.Ham_0
    Ham_I = aggOperators.Ham_I
    Ham_II_s = getInteractionHamIPicture(Ham_0, Ham_I, s)

    tmp1[:, :] = ad(rho, W0_bath.data, aggCore, aggTools)

    tmp2[:, :] = Ham_II_s.data * tmp1 - tmp1 * Ham_II_s.data
    tmp1[:, :] = Ham_II_t.data * tmp2 - tmp2 * Ham_II_t.data

    return trace_bath(tmp1, aggCore, aggTools; vib_basis=aggOperators.vib_basis)
end

function QME_sI_ansatz_const_sch(
    W0::T,
    tspan::Array,
    agg::Aggregate;
    reltol::AbstractFloat = 1.0e-12,
    abstol::AbstractFloat = 1.0e-12,
    int_reltol::AbstractFloat = 1.0e-4,
    int_abstol::AbstractFloat = 1.0e-4,
    alg::Any = DelayDiffEq.MethodOfSteps(DelayDiffEq.Vern6()),
    fout::Union{Function,Nothing} = nothing,
    kwargs...,
) where {B<:Basis,T<:Operator{B,B}}
    history_fun(p, t) = T(rho0.basis_l, rho0.basis_r, zeros(ComplexF64, size(rho0.data)))
    rho0 = trace_bath(W0, agg.core, agg.tools; vib_basis=agg.operators.vib_basis)
    W0_bath = get_rho_bath(W0, agg.core, agg.tools; vib_basis=agg.operators.vib_basis)
    p = (agg.core, agg.tools, agg.operators, W0, W0_bath, eltype(W0))

    tmp1 = copy(W0.data)
    tmp2 = copy(W0.data)
    dmaster_(t, rho, drho, history_fun, p) = dQME_sI_ansatz_const_sch(
        t,
        rho,
        drho,
        history_fun,
        tmp1,
        tmp2,
        p,
        int_reltol,
        int_abstol,
    )
    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
    x0 = rho0.data
    state = T(rho0.basis_l, rho0.basis_r, rho0.data)
    dstate = T(rho0.basis_l, rho0.basis_r, rho0.data)
    OpenQuantumSystems.integrate_delayed(
        tspan_,
        dmaster_,
        history_fun,
        x0,
        state,
        dstate,
        fout;
        p = p,
        reltol = reltol,
        abstol = abstol,
        alg = alg,
        kwargs...,
    )
end

function dQME_sI_ansatz_const_sch(
    t::AbstractFloat,
    rho::T,
    drho::T,
    history_fun,
    tmp1::Array,
    tmp2::Array,
    p,
    int_reltol::AbstractFloat,
    int_abstol::AbstractFloat,
) where {B<:Basis,T<:Operator{B,B}}
    aggCore, aggTools, aggOperators, W0, _, elementtype = p

    Ham_II_t = getInteractionHamIPicture(aggOperators.Ham_0, aggOperators.Ham_I, t)
    K = Ham_II_t.data * W0.data - W0.data * Ham_II_t.data
    K_traced = trace_bath(K, aggCore, aggTools; vib_basis=aggOperators.vib_basis)

    kernel_integrated_traced, err = QuadGK.quadgk(
        s -> kernel_sI_ansatz_const_sch(t, s, history_fun, p, tmp1, tmp2, Ham_II_t),
        0,
        t,
        rtol = int_reltol,
        atol = int_abstol,
    )
    drho.data[:, :] = -elementtype(im) * K_traced - kernel_integrated_traced

    return drho
end

function kernel_sI_ansatz_const_sch(t, s, h, p, tmp1, tmp2, Ham_II_t)
    aggCore, aggTools, aggOperators, W0, W0_bath, _ = p

    rho = h(p, s)

    if (typeof(rho) <: Operator)
        rho = rho.data
    end

    Ham_0 = aggOperators.Ham_0
    Ham_I = aggOperators.Ham_I
    Ham_II_s = getInteractionHamIPicture(Ham_0, Ham_I, s)

    U_0_op = evolutionOperator(Ham_0, s)
    W0_int_s = U_0_op' * W0_bath * U_0_op
    tmp1[:, :] = ad(rho, W0_int_s.data, aggCore, aggTools)

    tmp2[:, :] = Ham_II_s.data * tmp1 - tmp1 * Ham_II_s.data
    tmp1[:, :] = Ham_II_t.data * tmp2 - tmp2 * Ham_II_t.data

    return trace_bath(tmp1, aggCore, aggTools; vib_basis=aggOperators.vib_basis)
end

function QME_sI_ansatz_linear_sch(
    W0::T,
    tspan::Array,
    agg::Aggregate;
    reltol::AbstractFloat = 1.0e-12,
    abstol::AbstractFloat = 1.0e-12,
    int_reltol::AbstractFloat = 1.0e-4,
    int_abstol::AbstractFloat = 1.0e-4,
    alg::Any = DelayDiffEq.MethodOfSteps(DelayDiffEq.Vern6()),
    fout::Union{Function,Nothing} = nothing,
    kwargs...,
) where {B<:Basis,T<:Operator{B,B}}
    history_fun(p, t) = T(rho0.basis_l, rho0.basis_r, zeros(ComplexF64, size(rho0.data)))
    rho0 = trace_bath(W0, agg.core, agg.tools; vib_basis=agg.operators.vib_basis)
    W0_bath = get_rho_bath(W0, agg.core, agg.tools; vib_basis=agg.operators.vib_basis)
    p = (agg.core, agg.tools, agg.operators, W0, W0_bath, eltype(W0))

    tmp1 = copy(W0.data)
    tmp2 = copy(W0.data)
    dmaster_(t, rho, drho, history_fun, p) = dQME_sI_ansatz_linear_sch(
        t,
        rho,
        drho,
        history_fun,
        tmp1,
        tmp2,
        p,
        int_reltol,
        int_abstol,
    )
    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
    x0 = rho0.data
    state = T(rho0.basis_l, rho0.basis_r, rho0.data)
    dstate = T(rho0.basis_l, rho0.basis_r, rho0.data)
    OpenQuantumSystems.integrate_delayed(
        tspan_,
        dmaster_,
        history_fun,
        x0,
        state,
        dstate,
        fout;
        p = p,
        reltol = reltol,
        abstol = abstol,
        alg = alg,
        kwargs...,
    )
end

function dQME_sI_ansatz_linear_sch(
    t::AbstractFloat,
    rho::T,
    drho::T,
    history_fun,
    tmp1::Array,
    tmp2::Array,
    p,
    int_reltol::AbstractFloat,
    int_abstol::AbstractFloat,
) where {B<:Basis,T<:Operator{B,B}}
    aggCore, aggTools, aggOperators, W0, _, elementtype = p

    Ham_II_t = getInteractionHamIPicture(aggOperators.Ham_0, aggOperators.Ham_I, t)
    K = Ham_II_t.data * W0.data - W0.data * Ham_II_t.data
    K_traced = trace_bath(K, aggCore, aggTools; vib_basis=aggOperators.vib_basis)

    kernel_integrated_traced, err = QuadGK.quadgk(
        s -> kernel_sI_ansatz_linear_sch(t, s, history_fun, p, tmp1, tmp2, Ham_II_t),
        0,
        t,
        rtol = int_reltol,
        atol = int_abstol,
    )
    drho.data[:, :] = -elementtype(im) * K_traced - kernel_integrated_traced

    return drho
end

function kernel_sI_ansatz_linear_sch(t, s, h, p, tmp1, tmp2, Ham_II_t)
    aggCore, aggTools, aggOperators, _, W0_bath, elementtype = p

    rho = h(p, s)

    if (typeof(rho) <: Operator)
        rho = rho.data
    end

    Ham_0 = aggOperators.Ham_0
    Ham_I = aggOperators.Ham_I
    Ham = aggOperators.Ham
    Ham_II_s = getInteractionHamIPicture(Ham_0, Ham_I, s)

    U_0_op = evolutionOperator(Ham_0, s)
    # linear element
    W_bath_s = W0_bath - elementtype(im) * (Ham * W0_bath - W0_bath * Ham) * s
    W0_int_s = U_0_op' * W_bath_s * U_0_op
    tmp1[:, :] = ad(rho, W0_int_s.data, aggCore, aggTools)

    tmp2[:, :] = Ham_II_s.data * tmp1 - tmp1 * Ham_II_s.data
    tmp1[:, :] = Ham_II_t.data * tmp2 - tmp2 * Ham_II_t.data

    return trace_bath(tmp1, aggCore, aggTools; vib_basis=aggOperators.vib_basis)
end

function QME_sI_ansatz_linear2_sch(
    W0::T,
    tspan::Array,
    agg::Aggregate;
    reltol::AbstractFloat = 1.0e-12,
    abstol::AbstractFloat = 1.0e-12,
    int_reltol::AbstractFloat = 1.0e-4,
    int_abstol::AbstractFloat = 1.0e-4,
    t_mk_bath_count::Integer = 100,
    alg::Any = DelayDiffEq.MethodOfSteps(DelayDiffEq.Vern6()),
    fout::Union{Function,Nothing} = nothing,
    kwargs...,
) where {B<:Basis,T<:Operator{B,B}}
    history_fun(p, t) = T(rho0.basis_l, rho0.basis_r, zeros(ComplexF64, size(rho0.data)))
    rho0 = trace_bath(W0, agg.core, agg.tools; vib_basis=agg.operators.vib_basis)
    W0_bath = get_rho_bath(W0, agg.core, agg.tools; vib_basis=agg.operators.vib_basis)
    t_mk_bath_step = (tspan[end] - tspan[1]) / t_mk_bath_count
    p = (agg.core, agg.tools, agg.operators, W0, W0_bath, t_mk_bath_step, eltype(W0))

    tmp1 = copy(W0.data)
    tmp2 = copy(W0.data)
    dmaster_(t, rho, drho, history_fun, p) = dQME_sI_ansatz_linear2_sch(
        t,
        rho,
        drho,
        history_fun,
        tmp1,
        tmp2,
        p,
        int_reltol,
        int_abstol,
    )
    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
    x0 = rho0.data
    state = T(rho0.basis_l, rho0.basis_r, rho0.data)
    dstate = T(rho0.basis_l, rho0.basis_r, rho0.data)
    OpenQuantumSystems.integrate_delayed(
        tspan_,
        dmaster_,
        history_fun,
        x0,
        state,
        dstate,
        fout;
        p = p,
        reltol = reltol,
        abstol = abstol,
        alg = alg,
        kwargs...,
    )
end

function dQME_sI_ansatz_linear2_sch(
    t::AbstractFloat,
    rho::T,
    drho::T,
    history_fun,
    tmp1::Array,
    tmp2::Array,
    p,
    int_reltol::AbstractFloat,
    int_abstol::AbstractFloat,
) where {B<:Basis,T<:Operator{B,B}}
    aggCore, aggTools, aggOperators, W0, _, _, elementtype = p

    Ham_II_t = getInteractionHamIPicture(aggOperators.Ham_0, aggOperators.Ham_I, t)
    K = Ham_II_t.data * W0.data - W0.data * Ham_II_t.data
    K_traced = trace_bath(K, aggCore, aggTools; vib_basis=aggOperators.vib_basis)

    kernel_integrated_traced, err = QuadGK.quadgk(
        s -> kernel_sI_ansatz_linear2_sch(t, s, history_fun, p, tmp1, tmp2, Ham_II_t),
        0,
        t,
        rtol = int_reltol,
        atol = int_abstol,
    )
    drho.data[:, :] = -elementtype(im) * K_traced - kernel_integrated_traced

    return drho
end

function kernel_sI_ansatz_linear2_sch(t, s, h, p, tmp1, tmp2, Ham_II_t)
    aggCore, aggTools, aggOperators, _, W0_bath, t_mk_bath_step, elementtype = p

    rho = h(p, s)

    if (typeof(rho) <: Operator)
        rho = rho.data
    end

    Ham_0 = aggOperators.Ham_0
    Ham_I = aggOperators.Ham_I
    Ham = aggOperators.Ham
    Ham_II_s = getInteractionHamIPicture(Ham_0, Ham_I, s)

    U_0_op = evolutionOperator(Ham_0, s)
    # linear element on intervals
    W_bath_s = deepcopy(W0_bath)
    s_rest = s
    while s_rest > t_mk_bath_step
        W_bath_s = W_bath_s - elementtype(im) * (Ham * W_bath_s - W_bath_s * Ham) * t_mk_bath_step
        s_rest -= t_mk_bath_step
    end
    if s_rest > 0
        W_bath_s = W_bath_s - elementtype(im) * (Ham * W_bath_s - W_bath_s * Ham) * s_rest
    end

    W0_int_s = U_0_op' * W_bath_s * U_0_op
    tmp1[:, :] = ad(rho, W0_int_s.data, aggCore, aggTools)

    tmp2[:, :] = Ham_II_s.data * tmp1 - tmp1 * Ham_II_s.data
    tmp1[:, :] = Ham_II_t.data * tmp2 - tmp2 * Ham_II_t.data

    return trace_bath(tmp1, aggCore, aggTools; vib_basis=aggOperators.vib_basis)
end

function QME_sI_ansatz_upart1_sch(
    W0::T,
    tspan::Array,
    agg::Aggregate;
    reltol::AbstractFloat = 1.0e-12,
    abstol::AbstractFloat = 1.0e-12,
    int_reltol::AbstractFloat = 1.0e-4,
    int_abstol::AbstractFloat = 1.0e-4,
    alg::Any = DelayDiffEq.MethodOfSteps(DelayDiffEq.Vern6()),
    fout::Union{Function,Nothing} = nothing,
    kwargs...,
) where {B<:Basis,T<:Operator{B,B}}
    history_fun(p, t) = T(rho0.basis_l, rho0.basis_r, zeros(ComplexF64, size(rho0.data)))
    rho0 = trace_bath(W0, agg.core, agg.tools; vib_basis=agg.operators.vib_basis)
    W0_bath = get_rho_bath(W0, agg.core, agg.tools; vib_basis=agg.operators.vib_basis)
    p = (agg.core, agg.tools, agg.operators, W0, W0_bath, eltype(W0))

    tmp1 = copy(W0.data)
    tmp2 = copy(W0.data)
    dmaster_(t, rho, drho, history_fun, p) = dQME_sI_ansatz_upart1_sch(
        t,
        rho,
        drho,
        history_fun,
        tmp1,
        tmp2,
        p,
        int_reltol,
        int_abstol,
    )
    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
    x0 = rho0.data
    state = T(rho0.basis_l, rho0.basis_r, rho0.data)
    dstate = T(rho0.basis_l, rho0.basis_r, rho0.data)
    OpenQuantumSystems.integrate_delayed(
        tspan_,
        dmaster_,
        history_fun,
        x0,
        state,
        dstate,
        fout;
        p = p,
        reltol = reltol,
        abstol = abstol,
        alg = alg,
        kwargs...,
    )
end

function dQME_sI_ansatz_upart1_sch(
    t::AbstractFloat,
    rho::T,
    drho::T,
    history_fun,
    tmp1::Array,
    tmp2::Array,
    p,
    int_reltol::AbstractFloat,
    int_abstol::AbstractFloat,
) where {B<:Basis,T<:Operator{B,B}}
    aggCore, aggTools, aggOperators, W0, _, elementtype = p

    Ham_II_t = getInteractionHamIPicture(aggOperators.Ham_0, aggOperators.Ham_I, t)
    K = Ham_II_t.data * W0.data - W0.data * Ham_II_t.data
    K_traced = trace_bath(K, aggCore, aggTools; vib_basis=aggOperators.vib_basis)

    kernel_integrated_traced, err = QuadGK.quadgk(
        s -> kernel_sI_ansatz_upart1_sch(t, s, history_fun, p, tmp1, tmp2, Ham_II_t),
        0,
        t,
        rtol = int_reltol,
        atol = int_abstol,
    )
    drho.data[:, :] = -elementtype(im) * K_traced - kernel_integrated_traced

    return drho
end

function kernel_sI_ansatz_upart1_sch(t, s, h, p, tmp1, tmp2, Ham_II_t)
    aggCore, aggTools, aggOperators, _, W0_bath, elementtype = p

    rho = h(p, s)

    if (typeof(rho) <: Operator)
        rho = rho.data
    end

    indicesMap = aggTools.indicesMap
    Ham_0 = aggOperators.Ham_0
    Ham_I = aggOperators.Ham_I
    Ham = aggOperators.Ham
    Ham_II_s = getInteractionHamIPicture(Ham_0, Ham_I, s)

    U_0_op = evolutionOperator(Ham_0, s)
    # linear element on intervals
    W_bath_s = deepcopy(W0_bath)
    for el_st in 1:aggCore.molCount+1
        a1 = indicesMap[el_st][1]
        a2 = indicesMap[el_st][end]
        U_el_part = evolution_el_part(Ham.data, s, el_st, el_st, indicesMap)
        W_bath_s.data[a1:a2, a1:a2] = U_el_part * W_bath_s.data[a1:a2, a1:a2] * adjoint(U_el_part)
    end

    W0_int_s = U_0_op' * W_bath_s * U_0_op
    tmp1[:, :] = ad(rho, W0_int_s.data, aggCore, aggTools)

    tmp2[:, :] = Ham_II_s.data * tmp1 - tmp1 * Ham_II_s.data
    tmp1[:, :] = Ham_II_t.data * tmp2 - tmp2 * Ham_II_t.data

    return trace_bath(tmp1, aggCore, aggTools; vib_basis=aggOperators.vib_basis)
end

function QME_sI_ansatz_upart1_int(
    W0::T,
    tspan::Array,
    agg::Aggregate;
    reltol::AbstractFloat = 1.0e-12,
    abstol::AbstractFloat = 1.0e-12,
    int_reltol::AbstractFloat = 1.0e-4,
    int_abstol::AbstractFloat = 1.0e-4,
    alg::Any = DelayDiffEq.MethodOfSteps(DelayDiffEq.Vern6()),
    fout::Union{Function,Nothing} = nothing,
    kwargs...,
) where {B<:Basis,T<:Operator{B,B}}
    history_fun(p, t) = T(rho0.basis_l, rho0.basis_r, zeros(ComplexF64, size(rho0.data)))
    rho0 = trace_bath(W0, agg.core, agg.tools; vib_basis=agg.operators.vib_basis)
    W0_bath = get_rho_bath(W0, agg.core, agg.tools; vib_basis=agg.operators.vib_basis)
    p = (agg.core, agg.tools, agg.operators, W0, W0_bath, eltype(W0))

    tmp1 = copy(W0.data)
    tmp2 = copy(W0.data)
    dmaster_(t, rho, drho, history_fun, p) = dQME_sI_ansatz_upart1_int(
        t,
        rho,
        drho,
        history_fun,
        tmp1,
        tmp2,
        p,
        int_reltol,
        int_abstol,
    )
    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
    x0 = rho0.data
    state = T(rho0.basis_l, rho0.basis_r, rho0.data)
    dstate = T(rho0.basis_l, rho0.basis_r, rho0.data)
    OpenQuantumSystems.integrate_delayed(
        tspan_,
        dmaster_,
        history_fun,
        x0,
        state,
        dstate,
        fout;
        p = p,
        reltol = reltol,
        abstol = abstol,
        alg = alg,
        kwargs...,
    )
end

function dQME_sI_ansatz_upart1_int(
    t::AbstractFloat,
    rho::T,
    drho::T,
    history_fun,
    tmp1::Array,
    tmp2::Array,
    p,
    int_reltol::AbstractFloat,
    int_abstol::AbstractFloat,
) where {B<:Basis,T<:Operator{B,B}}
    aggCore, aggTools, aggOperators, W0, _, elementtype = p

    Ham_II_t = getInteractionHamIPicture(aggOperators.Ham_0, aggOperators.Ham_I, t)
    K = Ham_II_t.data * W0.data - W0.data * Ham_II_t.data
    K_traced = trace_bath(K, aggCore, aggTools; vib_basis=aggOperators.vib_basis)

    kernel_integrated_traced, err = QuadGK.quadgk(
        s -> kernel_sI_ansatz_upart1_int(t, s, history_fun, p, tmp1, tmp2, Ham_II_t),
        0,
        t,
        rtol = int_reltol,
        atol = int_abstol,
    )
    drho.data[:, :] = -elementtype(im) * K_traced - kernel_integrated_traced

    return drho
end

function kernel_sI_ansatz_upart1_int(t, s, h, p, tmp1, tmp2, Ham_II_t)
    aggCore, aggTools, aggOperators, _, W0_bath, elementtype = p

    rho = h(p, s)

    if (typeof(rho) <: Operator)
        rho = rho.data
    end

    indicesMap = aggTools.indicesMap
    Ham_0 = aggOperators.Ham_0
    Ham_I = aggOperators.Ham_I
    Ham = aggOperators.Ham
    Ham_II_s = getInteractionHamIPicture(Ham_0, Ham_I, s)

    # linear element on intervals
    W_bath_s = deepcopy(W0_bath)
    for el_st in 1:aggCore.molCount+1
        a1 = indicesMap[el_st][1]
        a2 = indicesMap[el_st][end]
        U_el_part = evolution_el_part(Ham.data, s, el_st, el_st, indicesMap)
        W_bath_s.data[a1:a2, a1:a2] = U_el_part * W_bath_s.data[a1:a2, a1:a2] * adjoint(U_el_part)
    end

    tmp1[:, :] = ad(rho, W_bath_s.data, aggCore, aggTools)

    tmp2[:, :] = Ham_II_s.data * tmp1 - tmp1 * Ham_II_s.data
    tmp1[:, :] = Ham_II_t.data * tmp2 - tmp2 * Ham_II_t.data

    return trace_bath(tmp1, aggCore, aggTools; vib_basis=aggOperators.vib_basis)
end

function QME_sI_ansatz_upart2_sch(
    W0::T,
    tspan::Array,
    agg::Aggregate;
    reltol::AbstractFloat = 1.0e-12,
    abstol::AbstractFloat = 1.0e-12,
    int_reltol::AbstractFloat = 1.0e-4,
    int_abstol::AbstractFloat = 1.0e-4,
    alg::Any = DelayDiffEq.MethodOfSteps(DelayDiffEq.Vern6()),
    fout::Union{Function,Nothing} = nothing,
    kwargs...,
) where {B<:Basis,T<:Operator{B,B}}
    history_fun(p, t) = T(rho0.basis_l, rho0.basis_r, zeros(ComplexF64, size(rho0.data)))
    rho0 = trace_bath(W0, agg.core, agg.tools; vib_basis=agg.operators.vib_basis)
    W0_bath = get_rho_bath(W0, agg.core, agg.tools; vib_basis=agg.operators.vib_basis)
    p = (agg.core, agg.tools, agg.operators, W0, W0_bath, eltype(W0))

    tmp1 = copy(W0.data)
    tmp2 = copy(W0.data)
    dmaster_(t, rho, drho, history_fun, p) = dQME_sI_ansatz_upart2_sch(
        t,
        rho,
        drho,
        history_fun,
        tmp1,
        tmp2,
        p,
        int_reltol,
        int_abstol,
    )
    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
    x0 = rho0.data
    state = T(rho0.basis_l, rho0.basis_r, rho0.data)
    dstate = T(rho0.basis_l, rho0.basis_r, rho0.data)
    OpenQuantumSystems.integrate_delayed(
        tspan_,
        dmaster_,
        history_fun,
        x0,
        state,
        dstate,
        fout;
        p = p,
        reltol = reltol,
        abstol = abstol,
        alg = alg,
        kwargs...,
    )
end

function dQME_sI_ansatz_upart2_sch(
    t::AbstractFloat,
    rho::T,
    drho::T,
    history_fun,
    tmp1::Array,
    tmp2::Array,
    p,
    int_reltol::AbstractFloat,
    int_abstol::AbstractFloat,
) where {B<:Basis,T<:Operator{B,B}}
    aggCore, aggTools, aggOperators, W0, _, elementtype = p

    Ham_II_t = getInteractionHamIPicture(aggOperators.Ham_0, aggOperators.Ham_I, t)
    K = Ham_II_t.data * W0.data - W0.data * Ham_II_t.data
    K_traced = trace_bath(K, aggCore, aggTools; vib_basis=aggOperators.vib_basis)

    kernel_integrated_traced, err = QuadGK.quadgk(
        s -> kernel_sI_ansatz_upart2_sch(t, s, history_fun, p, tmp1, tmp2, Ham_II_t),
        0,
        t,
        rtol = int_reltol,
        atol = int_abstol,
    )
    drho.data[:, :] = -elementtype(im) * K_traced - kernel_integrated_traced

    return drho
end

function kernel_sI_ansatz_upart2_sch(t, s, h, p, tmp1, tmp2, Ham_II_t)
    aggCore, aggTools, aggOperators, _, W0_bath, elementtype = p

    rho = h(p, s)

    if (typeof(rho) <: Operator)
        rho = rho.data
    end

    indicesMap = aggTools.indicesMap
    Ham_0 = aggOperators.Ham_0
    Ham_I = aggOperators.Ham_I
    Ham = aggOperators.Ham
    Ham_II_s = getInteractionHamIPicture(Ham_0, Ham_I, s)

    U_0_op = evolutionOperator(Ham_0, s)
    # linear element on intervals
    W_bath_s = deepcopy(W0_bath)
    for el_st1 in 1:aggCore.molCount+1
        a1 = indicesMap[el_st1][1]
        a2 = indicesMap[el_st1][end]
        for el_st2 in 1:aggCore.molCount+1
            b1 = indicesMap[el_st2][1]
            b2 = indicesMap[el_st2][end]
            U_el_part = evolution_el_part(Ham.data, s, el_st1, el_st2, indicesMap)
            W_bath_s.data[a1:a2, b1:b2] = U_el_part * W_bath_s.data[a1:a2, b1:b2] * adjoint(U_el_part)
        end
    end

    W0_int_s = U_0_op' * W_bath_s * U_0_op
    tmp1[:, :] = ad(rho, W0_int_s.data, aggCore, aggTools)

    tmp2[:, :] = Ham_II_s.data * tmp1 - tmp1 * Ham_II_s.data
    tmp1[:, :] = Ham_II_t.data * tmp2 - tmp2 * Ham_II_t.data

    return trace_bath(tmp1, aggCore, aggTools; vib_basis=aggOperators.vib_basis)
end

function QME_sI_ansatz_upart2_int(
    W0::T,
    tspan::Array,
    agg::Aggregate;
    reltol::AbstractFloat = 1.0e-12,
    abstol::AbstractFloat = 1.0e-12,
    int_reltol::AbstractFloat = 1.0e-4,
    int_abstol::AbstractFloat = 1.0e-4,
    alg::Any = DelayDiffEq.MethodOfSteps(DelayDiffEq.Vern6()),
    fout::Union{Function,Nothing} = nothing,
    kwargs...,
) where {B<:Basis,T<:Operator{B,B}}
    history_fun(p, t) = T(rho0.basis_l, rho0.basis_r, zeros(ComplexF64, size(rho0.data)))
    rho0 = trace_bath(W0, agg.core, agg.tools; vib_basis=agg.operators.vib_basis)
    W0_bath = get_rho_bath(W0, agg.core, agg.tools; vib_basis=agg.operators.vib_basis)
    p = (agg.core, agg.tools, agg.operators, W0, W0_bath, eltype(W0))

    tmp1 = copy(W0.data)
    tmp2 = copy(W0.data)
    dmaster_(t, rho, drho, history_fun, p) = dQME_sI_ansatz_upart2_int(
        t,
        rho,
        drho,
        history_fun,
        tmp1,
        tmp2,
        p,
        int_reltol,
        int_abstol,
    )
    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
    x0 = rho0.data
    state = T(rho0.basis_l, rho0.basis_r, rho0.data)
    dstate = T(rho0.basis_l, rho0.basis_r, rho0.data)
    OpenQuantumSystems.integrate_delayed(
        tspan_,
        dmaster_,
        history_fun,
        x0,
        state,
        dstate,
        fout;
        p = p,
        reltol = reltol,
        abstol = abstol,
        alg = alg,
        kwargs...,
    )
end

function dQME_sI_ansatz_upart2_int(
    t::AbstractFloat,
    rho::T,
    drho::T,
    history_fun,
    tmp1::Array,
    tmp2::Array,
    p,
    int_reltol::AbstractFloat,
    int_abstol::AbstractFloat,
) where {B<:Basis,T<:Operator{B,B}}
    aggCore, aggTools, aggOperators, W0, _, elementtype = p

    Ham_II_t = getInteractionHamIPicture(aggOperators.Ham_0, aggOperators.Ham_I, t)
    K = Ham_II_t.data * W0.data - W0.data * Ham_II_t.data
    K_traced = trace_bath(K, aggCore, aggTools; vib_basis=aggOperators.vib_basis)

    kernel_integrated_traced, err = QuadGK.quadgk(
        s -> kernel_sI_ansatz_upart2_int(t, s, history_fun, p, tmp1, tmp2, Ham_II_t),
        0,
        t,
        rtol = int_reltol,
        atol = int_abstol,
    )
    drho.data[:, :] = -elementtype(im) * K_traced - kernel_integrated_traced

    return drho
end

function kernel_sI_ansatz_upart2_int(t, s, h, p, tmp1, tmp2, Ham_II_t)
    aggCore, aggTools, aggOperators, _, W0_bath, elementtype = p

    rho = h(p, s)

    if (typeof(rho) <: Operator)
        rho = rho.data
    end

    indicesMap = aggTools.indicesMap
    Ham_0 = aggOperators.Ham_0
    Ham_I = aggOperators.Ham_I
    Ham = aggOperators.Ham
    Ham_II_s = getInteractionHamIPicture(Ham_0, Ham_I, s)

    # linear element on intervals
    W_bath_s = deepcopy(W0_bath)
    for el_st1 in 1:aggCore.molCount+1
        a1 = indicesMap[el_st1][1]
        a2 = indicesMap[el_st1][end]
        for el_st2 in 1:aggCore.molCount+1
            b1 = indicesMap[el_st2][1]
            b2 = indicesMap[el_st2][end]
            U_el_part = evolution_el_part(Ham.data, s, el_st1, el_st2, indicesMap)
            W_bath_s.data[a1:a2, b1:b2] = U_el_part * W_bath_s.data[a1:a2, b1:b2] * adjoint(U_el_part)
        end
    end

    tmp1[:, :] = ad(rho, W_bath_s.data, aggCore, aggTools)

    tmp2[:, :] = Ham_II_s.data * tmp1 - tmp1 * Ham_II_s.data
    tmp1[:, :] = Ham_II_t.data * tmp2 - tmp2 * Ham_II_t.data

    return trace_bath(tmp1, aggCore, aggTools; vib_basis=aggOperators.vib_basis)
end

max_stars_repo_name: detrin/OpenQuantumSystems

max_issues_repo_name: detrin/OpenQuantumSystems

max_forks_repo_name: detrin/OpenQuantumSystems


text:
#### naive sampling

struct CategoricalDirectSampler{T<:Real,Ts<:AbstractVector{T}} <: Sampleable{Univariate,Discrete}
    prob::Ts

    function CategoricalDirectSampler{T,Ts}(p::Ts) where {
        T<:Real,Ts<:AbstractVector{T}}
        isempty(p) && throw(ArgumentError("p is empty."))
        new{T,Ts}(p)
    end
end

CategoricalDirectSampler(p::Ts) where {T<:Real,Ts<:AbstractVector{T}} =
    CategoricalDirectSampler{T,Ts}(p)

ncategories(s::CategoricalDirectSampler) = length(s.prob)

function rand(rng::AbstractRNG, s::CategoricalDirectSampler)
    p = s.prob
    n = length(p)
    i = 1
    c = p[1]
    u = rand(rng)
    while c < u && i < n
        c += p[i += 1]
    end
    return i
end


text:
function _parse_typestring(s::AbstractString)
    s = strip(s)

    (isempty(s)||s[1]=='(') && error("typename parse error in $S")
    brace_pos = findfirst("(", s)
    if isnothing(brace_pos)
        type_name = Symbol(s)
        return is_ch_type(type_name) ? TypeAst(type_name) : s

    end
    brace_pos = brace_pos[1]
    s[end] != ')' && error("typename parse error in $s")
    type_name = Symbol(strip(s[1:brace_pos-1]))
    !is_ch_type(type_name) && error("typename parse error in $s")

    ast = TypeAst(type_name)
    inner = SubString(s, brace_pos + 1, length(s) - 1)

    cursor = 1
    elem_pos = 1
    opened_braces = 0
    while true
        range = findnext(r"\(|\)|,", inner, cursor)

        if isnothing(range) || isempty(range)

            if elem_pos <= length(inner)
                push!(ast, _parse_typestring(inner[elem_pos:end]))
            end
            break
        end
        pos = first(range)
        inner[pos] == '(' && (opened_braces += 1)
        inner[pos] == ')' && (opened_braces -= 1)
        if inner[pos] == ',' && opened_braces == 0
            push!(ast, _parse_typestring(inner[elem_pos:pos-1]))
            elem_pos = last(range) + 1
        end
        cursor = last(range) + 1
    end
    return ast
end

function parse_typestring(s::AbstractString)
    result = _parse_typestring(s)
    !(result isa TypeAst) && error("typename parse error in $s")
    return result
end