text:
using DifferentialEquations
using LSODA
using JLD2
using ArgParse
using StatsBase
using Distributions
using ThreadTools

parser = ArgParseSettings(allow_ambiguous_opts=false)
@add_arg_table! parser begin
    "output"
        help = "output JLD2"
        required = true
    "pset"
        help = "parameter set: snic, sn, hopf, sl"
        required = true
    "form"
        help = "noise form: add or mul"
        required = true
    "initstate"
        help = "initial state: hi or lo"
        required = true
    "signals"
        help = "signal values"
        required = true
    "--noises"
        help = "noise values"
        default = "0.5:0.5:1.5"
    "--stochtime"
        help = "stochastic simulation endpoint"
        arg_type = Float64
        default = 50.0
    "--initcells"
        help = "number of cells selected"
        arg_type = Int
        default = 4
end
args = parse_args(parser)

cellspecs = Dict("add" => 9, "mul" => 12)
const cellspec = cellspecs[args["form"]]
const divlimit = 75.0

function model_mmi2_add(du, u, p, t)
    sR, a1, a2, b1, b2, kR, kr, K1, K2, koff, sr, sR0, kG, σr, σR, σG, τ, V = p
    r, R, c1A, c1B, c2 = u[1:cellspec:end], u[2:cellspec:end], u[3:cellspec:end], u[4:cellspec:end], u[5:cellspec:end]
    sig = sR

    @. du[1:cellspec:end] = sr - kr * r + koff * (c1A+c1B) - K1 * koff * 2R * r +
                koff * 2c2 - K2 * koff * (c1A+c1B) * r +
                (c1A+c1B) * kR * a1 + 2c2 * kR * a2 + u[7:cellspec:end]

    @. du[2:cellspec:end] = sig + sR0 - kR * R + koff * (c1A+c1B) - K1 * koff * 2R * r +
                (c1A+c1B) * kr * b1 + u[8:cellspec:end]

    @. du[3:cellspec:end] = K1 * koff * R * r - koff * c1A +
                koff * c2 - K2 * koff * c1A * r +
                c2 * 1kr * b2 - c1A * kR * a1 - c1A * kr * b1

    @. du[4:cellspec:end] = K1 * koff * R * r - koff * c1B +
                koff * c2 - K2 * koff * c1B * r +
                c2 * 1kr * b2 - c1B * kR * a1 - c1B * kr * b1

    @. du[5:cellspec:end] = K2 * koff * (c1A+c1B) * r - koff * 2c2 -
                c2 * 2kr * b2 - c2 * kR * a2

    @. du[6:cellspec:end] = kG + u[9:cellspec:end]

    for i in 7:9
        du[i:cellspec:end] = -u[i:cellspec:end]/τ
    end
end

function stoch_mmi2_add(du, u, p, t)
    sR, a1, a2, b1, b2, kR, kr, K1, K2, koff, sr, sR0, kG, σr, σR, σG, τ, V = p
    for i in 1:6
        du[i:cellspec:end] .= 0
    end
    du[7:cellspec:end] .= sqrt(2.0*σr^2.0/τ)
    du[8:cellspec:end] .= sqrt(2.0*σR^2.0/τ)
    du[9:cellspec:end] .= sqrt(2.0*σG^2.0/τ)
end

function model_mmi2_mul(du, u, p, t)
    sR, a1, a2, b1, b2, kR, kr, K1, K2, koff, sr, sR0, kG, σr, σR, σG, τ, V = p
    r, R, c1A, c1B, c2 = u[1:cellspec:end], u[2:cellspec:end], u[3:cellspec:end], u[4:cellspec:end], u[5:cellspec:end]
    sig = sR

    @. du[1:cellspec:end] = sr - kr * r + koff * (c1A+c1B) - K1 * koff * 2R * r +
                koff * 2c2 - K2 * koff * (c1A+c1B) * r +
                (c1A+c1B) * kR * a1 + 2c2 * kR * a2 + u[7:cellspec:end]

    @. du[2:cellspec:end] = sig + sR0 - kR * R + koff * (c1A+c1B) - K1 * koff * 2R * r +
                (c1A+c1B) * kr * b1 + u[8:cellspec:end]

    @. du[3:cellspec:end] = K1 * koff * R * r - koff * c1A +
                koff * c2 - K2 * koff * c1A * r +
                c2 * 1kr * b2 - c1A * kR * a1 - c1A * kr * b1 + u[9:cellspec:end]

    @. du[4:cellspec:end] = K1 * koff * R * r - koff * c1B +
                koff * c2 - K2 * koff * c1B * r +
                c2 * 1kr * b2 - c1B * kR * a1 - c1B * kr * b1 + u[10:cellspec:end]

    @. du[5:cellspec:end] = K2 * koff * (c1A+c1B) * r - koff * 2c2 -
                c2 * 2kr * b2 - c2 * kR * a2 + u[11:cellspec:end]

    @. du[6:cellspec:end] = kG + u[12:cellspec:end]

    for i in 7:12
        du[i:cellspec:end] = -u[i:cellspec:end]/τ
    end
end

function stoch_mmi2_mul(du, u, p, t)
    sR, a1, a2, b1, b2, kR, kr, K1, K2, koff, sr, sR0, kG, σr, σR, σG, τ, V = p
    for i in 1:6
        du[i:cellspec:end] .= 0
    end
    for i in 7:11
        du[i:cellspec:end] = sqrt(2.0*σR^2.0/τ) * u[(i - 6):cellspec:end]
    end
    du[12:cellspec:end] .= sqrt(2.0*σG^2.0/τ)
end

function condition_mmi2_split(u, t, integrator)
    for i in 1:5
        minconc = minimum(u[i:cellspec:end])
        if minconc < -0.01 # Adjust if needed
            println("Possibly diverging conc of $minconc for species $i at t = $t")
            terminate!(integrator)
        end
    end
    return maximum(u[6:cellspec:end]) >= divlimit
end

function affect_mmi2_split(integrator)
    #println("Affecting at t = $(integrator.t)")
    length(integrator.u) > 5000 * cellspec && error("Runaway division")
    sR, a1, a2, b1, b2, kR, kr, K1, K2, koff, sr, sR0, kG, σr, σR, σG, τ, V = integrator.p
    while maximum(integrator.u[6:cellspec:end]) >= divlimit
        prevend = length(integrator.u)
        maxcell = argmax(integrator.u[6:cellspec:end])
        maxbase = (maxcell - 1) * cellspec
        resize!(integrator, prevend + cellspec)
        u = integrator.u
        u[maxbase + 6] >= divlimit || error("Index $(maxbase + 6) should be above division threshold")
        u[maxbase + 6] = 0
        for spec in 1:5
            conc = u[maxbase + spec]
            #print("Species $spec conc was $conc, ")
            if conc > 0
                dist = truncated(Normal(0.5, 0.5 / sqrt(2 * conc * V)), 0., 1.)
                newconc = 2 * rand(dist) * conc
                u[end - cellspec + spec] = newconc
                u[maxbase + spec] = 2 * conc - newconc
                #println("now $newconc and $(2 * conc - newconc)")
            else
                #println("unchanged (< 0)")
                u[end - cellspec + spec] = conc
            end
        end
        u[(end - cellspec + 7):end] = u[(maxbase + 7):(maxbase + cellspec)]
        #println("Split cell $maxcell, now $(length(u) / cellspec |> Int) cells")
    end
    return nothing
end

callback_mmi2_split = DiscreteCallback(condition_mmi2_split, affect_mmi2_split)

p_snic = [ # SNIC, like Figure 5B
    0.3, #1 sR
    1.0, #2 a1
    12.0, #3 a2
    1.0, #4 b1
    4.0, #5 b2
    1.0, #6 kR
    0.25, #7 kr
    1000.0, #8 K1
    1000.0, #9 K2
    100.0, #10 koff
    1.0, #11 sr
    5.7, #12 sR0    total 6
    divlimit / 6.0, #13 kG
    0.25, #14 σr    replaced
    1.00, #15 σR    replaced
    1.00, #16 σG    replaced
    0.1, #17 τ
    2000.0 #18 V
   ];

p_hopf = [ # Hopf only, like Figure 3B
    0.5, #1 sR
    1.0, #2 a1
    12.0, #3 a2
    1.0, #4 b1
    7.0, #5 b2
    1.0, #6 kR
    0.25, #7 kr
    1000.0, #8 K1
    1000.0, #9 K2
    100.0, #10 koff
    1.0, #11 sr
    3.1, #12 sR0    total 3.6
    divlimit / 6.0, #13 kG
    0.25, #14 σr    replaced
    1.00, #15 σR    replaced
    1.00, #16 σG    replaced
    0.1, #17 τ
    2000.0 #18 V
    ];

p_sn = [ # SN, like Figure 3F
    3.0, #1 sR
    1.0, #2 a1
    4.0, #3 a2
    0.5, #4 b1
    0.1, #5 b2
    1.0, #6 kR
    2.0, #7 kr
    1000.0, #8 K1
    1000.0, #9 K2
    100.0, #10 koff
    1.0, #11 sr
    0.0, #12 sR0    total 3
    divlimit / 6.0, #13 kG
    0.25, #14 σr    replaced
    1.00, #15 σR    replaced
    1.00, #16 σG    replaced
    0.1, #17 τ
    2000.0 #18 V
   ];

p_sl = [ # Saddle-loop, like Figure 5G
    1.0, #1 sR
    1.0, #2 a1
    12.0, #3 a2
    1.0, #4 b1
    3.0, #5 b2
    1.0, #6 kR
    0.25, #7 kr
    1000.0, #8 K1
    1000.0, #9 K2
    100.0, #10 koff
    1.0, #11 sr
    6.9, #12 sR0    total 7.9
    divlimit / 6.0, #13 kG
    0.25, #14 σr    replaced
    1.00, #15 σR    replaced
    1.00, #16 σG    replaced
    0.1, #17 τ
    2000.0 #18 V
  ];

function timecourse_array(sols, savevars, savedat=0.1)
    timepoints = 0:savedat:sols[end].t[end]
    maxcells = map(s -> length(s[end]) / cellspec, sols) |> maximum |> Int
    timecourse = fill!(Array{Float64, 3}(undef, length(savevars), length(timepoints), maxcells), NaN)
    for (tindex, time) in enumerate(timepoints)
        for (sindex, sol) in enumerate(sols)
            (time < sol.t[begin] || time > sol.t[end]) && continue
            state = sol(time)
            curcells = length(state) / cellspec |> Int
            timecourse[:, tindex, :] .= NaN
            for (target, source) in enumerate(savevars)
                timecourse[target, tindex, 1:curcells] = state[source:cellspec:end]
            end
        end
    end
    return timecourse
end

function sim_pop(model, model_stoch, u0, pars, initcells; selectR::String="lo", saveat=0.1)
    selectR == "lo" || selectR == "hi" || error("invalid selectR")
    pars_nogrowth = copy(pars)
    pars_nogrowth[13] = 0 # kG = 0
    prob_ode = ODEProblem(model, u0, (0., 200.), pars_nogrowth);
    r0s = range(0.001, 10.0, length=20)
    R0s = range(0.001, 10.0, length=20)
    function prob_func_cells_init(prob,i,repeat)
        m, n = i ÷ length(r0s), i % length(r0s)
        if n != 0
            r0, R0 = r0s[m+1], R0s[n]
        else
            r0, R0 = r0s[m], R0s[end]
        end
        remake(prob, u0=[r0; R0; u0[3:end]])
    end
    prob_cells_init = EnsembleProblem(prob_ode, prob_func=prob_func_cells_init)
    @time sim_init = solve(prob_cells_init, lsoda(), saveat=saveat, dt=0.01,
        reltol=1e-6, abstol=1e-6, trajectories=length(r0s)*length(R0s));
    sol_init = sim_init[:, 400:end, :] # The selection of state is sensitive to the start index
    sol_reduced = dropdims(sum(sol_init[2:5, :, :]; dims=1), dims=1)
    if selectR == "lo"
        sol_reduced[isnan.(sol_reduced)] .= Inf
        idsel = argmin(sol_reduced)
    else
        sol_reduced[isnan.(sol_reduced)] .= -Inf
        idsel = argmax(sol_reduced)
    end
    u0_single = sol_init[:, idsel[1], idsel[2]]
    println("Starting stochastic simulation with initial condition ", u0_single)
    u0_cells = repeat(u0_single, initcells)
    u0_cells[6:cellspec:end] = rand(Uniform(0, divlimit), initcells)
    prob = SDEProblem(model, model_stoch, u0_cells, tspan, pars)
    sol = solve(prob, callback=callback_mmi2_split, saveat=saveat, dt=0.001, maxiter=1e12);
    sols = [sol]
    while sols[end].t[end] < tspan[2]
        println("Interrupted at t = $(sols[end].t[end]) with status $(sol.retcode)")
        if sol.retcode == :Terminated
            fullcourse = timecourse_array(sols, 1:cellspec, saveat)
            stopindex = size(fullcourse)[2] # Adjust savepoint selection if needed
            while minimum(filter(!isnan, fullcourse[1:5, stopindex, :])) < -0.0005
                stopindex -= 1
            end
            stopindex -= 1
            trest = ((stopindex - 1) * saveat, tspan[2])
            println("Using savepoint t = $(trest[1])")
            ulast = filter(!isnan, fullcourse[:, stopindex, :])
        else
            stopindex = length(sols[end])
            while !isinteger(round(sols[end].t[stopindex] / saveat; digits=5))
                stopindex = stopindex - 1
            end
            println("Using endpoint t = $(sols[end].t[stopindex])")
            ulast = sols[end].u[stopindex]
            trest = (sols[end].t[stopindex], tspan[2])
        end
        prob_resume = remake(prob, u0=ulast, tspan=trest)
        sol = solve(prob_resume, callback=callback_mmi2_split, saveat=saveat, dt=0.001, maxiter=1e12)
        push!(sols, sol)
    end
    return timecourse_array(sols, 1:5, saveat)
end

u0s = Dict("add" => [0.8, 0.1, 0., 0., 0., 0., 0., 0., 0.], "mul" => [0.8, 0.1, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
models = Dict("add" => model_mmi2_add, "mul" => model_mmi2_mul)
stochs = Dict("add" => stoch_mmi2_add, "mul" => stoch_mmi2_mul)
psets = Dict("snic" => p_snic, "sn" => p_sn, "hopf" => p_hopf, "sl" => p_sl)

model = models[args["form"]]
stoch = stochs[args["form"]]
u0 = u0s[args["form"]]
p = psets[args["pset"]]
tspan = (0., args["stochtime"])
noises = eval(Meta.parse(args["noises"]))
signals = eval(Meta.parse(args["signals"]))

all_results = map(noises) do noise
    results = map(signals) do signal
        signal => sim_pop(model, stoch, u0, [signal; p[2:end-5]; noise*0.25; noise; noise; p[(end-1):end]], args["initcells"]; selectR=args["initstate"])
    end
    noise => Dict(results)
end
save(args["output"], "results", Dict(all_results))


text:
struct Preference{T}
  possibles::Vector{T}
  min::T
  max::T
end

Preference(possibles::Vector{T}) where T <: Real = Preference(sort(unique(possibles)), minimum(possibles), maximum(possibles))
Preference(possibles::AbstractRange) = Preference(collect(possibles), minimum(possibles), maximum(possibles))

Base.string(x::Preference) = string(sort(x.possibles))
Base.print(io::IO, x::Preference) = print(io, string(x))
Base.show(io::IO, x::Preference) = print(io, "Ratings Preference: ", x)

"""
    minimum(preference::Preference)
Return lower value of rating preferences.
"""
Base.minimum(preference::Preference) = preference.min

"""
    maximum(preference::Preference)
Return highest value of rating preferences.
"""
Base.maximum(preference::Preference) = preference.max

"""
    size(preference::Preference)
Return number of elements contained in the rating preferences.
"""
Base.size(preference::Preference) = length(preference.possibles)

"""
    eltype(preference::Preference)
Return rating type.
"""
Base.eltype(preference::Preference) = eltype(preference.possibles)

"""
    round(rating::T, preference::Preference{T}) where T
Returns the nearest integral value of the rating preferences.
"""
function Base.round(rating::T, preference::Preference{T}) where T
  ratings = sort(preference.possibles)

  m = abs.(rating .- ratings)

  return ratings[findall(r->r == minimum(m), m)[end]]
end

Base.round(rating::Float64, preference::Preference{Int}) = round(Base.convert(Int, round(rating)), preference)
Base.round(rating, preference::Preference{T}) where T = round(Base.convert(T, rating), preference)

"""
    correct(rating::Real, preference::Preference)
Returns the correct value preserving the interval of rating
preferences.
"""
function correct(rating::Real, preference::Preference)
  if rating > maximum(preference)
    return maximum(preference)
  elseif rating < minimum(preference)
    return minimum(preference)
  end
  return rating
end

recommendation(preference::Preference) = maximum(preference) - round(1/3 * (maximum(preference) -  minimum(preference)))

ext: jl


text:
using Dates

# Widgets built on top of more basic widgets

"""
    frame(w) -> f

Return the GtkFrame `f` associated with widget `w`.
"""
frame(f::GtkFrame) = f

################# A movie-player widget ##################

struct Player{P} <: Widget
    signal::Signal{Int}
    widget::P
    preserved::Vector

    function Player{P}(signal::Signal{Int}, widget, preserved) where P
        obj = new{P}(signal, widget, preserved)
        gc_preserve(frame(widget), obj)
        obj
    end
end
Player(signal::Signal{Int}, widget::P, preserved) where {P} =
    Player{P}(signal, widget, preserved)

frame(p::Player) = frame(p.widget)

struct PlayerWithTextbox
    range::UnitRange{Int}     # valid values for index
    direction::Signal{Int8}   # +1 = forward, -1 = backward, 0 = not playing
    # GUI elements
    frame::GtkFrame
    scale::Slider{Int}
    entry::Textbox
    play_back::Button
    step_back::Button
    stop::Button
    step_forward::Button
    play_forward::Button
end

frame(p::PlayerWithTextbox) = p.frame

function PlayerWithTextbox(builder, index::Signal, range::AbstractUnitRange, id::Integer=1)
    1 <= id <= 2 || error("only 2 player widgets are defined in player.glade")
    direction = Signal(Int8(0))
    frame = Gtk.G_.object(builder,"player_frame$id")
    scale = slider(range; widget=Gtk.G_.object(builder,"index_scale$id"), signal=index)
    entry = textbox(first(range); widget=Gtk.G_.object(builder,"index_entry$id"), signal=index, range=range)
    play_back = button(; widget=Gtk.G_.object(builder,"play_back$id"))
    step_back = button(; widget=Gtk.G_.object(builder,"step_back$id"))
    stop = button(; widget=Gtk.G_.object(builder,"stop$id"))
    step_forward = button(; widget=Gtk.G_.object(builder,"step_forward$id"))
    play_forward = button(; widget=Gtk.G_.object(builder,"play_forward$id"))

    # Fix up widget properties
    set_gtk_property!(scale.widget, "round-digits", 0)  # glade/gtkbuilder bug that I have to set this here?

    # Link the buttons
    clampindex(i) = clamp(i, minimum(range), maximum(range))
    preserved = [map(x->push!(direction, -1), signal(play_back); init=nothing),
                 map(x->(push!(direction, 0); push!(index, clampindex(value(index)-1))),
                     signal(step_back); init=nothing),
                 map(x->push!(direction, 0), signal(stop); init=nothing),
                 map(x->(push!(direction, 0); push!(index, clampindex(value(index)+1))),
                     signal(step_forward); init=nothing),
                 map(x->push!(direction, +1), signal(play_forward); init=nothing)]
    function advance(widget)
        i = value(index) + value(direction)
        if !(i ∈ range)
            push!(direction, 0)
            i = clampindex(i)
        end
        push!(index, i)
        nothing
    end
    # Stop playing if the widget is destroyed
    signal_connect(frame, :destroy) do widget
        push!(direction, 0)
    end
    # Start the timer
    push!(preserved, map(advance, fpswhen(map(x->x!=0, direction), 30)))
    # Create the player object
    PlayerWithTextbox(range, direction, frame, scale, entry, play_back, step_back, stop, step_forward, play_forward), preserved
end
function PlayerWithTextbox(index::Signal, range::AbstractUnitRange, id::Integer=1)
    builder = GtkBuilder(filename=joinpath(splitdir(@__FILE__)[1], "player.glade"))
    PlayerWithTextbox(builder, index, range, id)
end

player(range::AbstractRange{Int}; style="with-textbox", id::Int=1) =
    player(Signal(first(range)), range; style=style, id=id)

"""
    player(range; style="with-textbox", id=1)
    player(slice::Signal{Int}, range; style="with-textbox", id=1)

Create a movie-player widget. This includes the standard play and stop
buttons and a slider; style "with-textbox" also includes play
backwards, step forward/backward, and a textbox for entering a
slice by keyboard.

You can create up to two player widgets for the same GUI, as long as
you pass `id=1` and `id=2`, respectively.
"""
function player(cs::Signal, range::AbstractUnitRange; style="with-textbox", id::Int=1)
    if style == "with-textbox"
        widget, preserved = PlayerWithTextbox(cs, range, id)
        return Player(cs, widget, preserved)
    end
    error("style $style not recognized")
end

Base.unsafe_convert(::Type{Ptr{Gtk.GLib.GObject}}, p::PlayerWithTextbox) =
    Base.unsafe_convert(Ptr{Gtk.GLib.GObject}, frame(p))

################# A time widget ##########################

struct TimeWidget{T <: Dates.TimeType} <: InputWidget{T}
    signal::Signal{T}
    widget::GtkFrame
end

"""
    timewidget(time)

Return a time widget that includes the `Time` and a `GtkFrame` with the hour, minute, and
second widgets in it. You can specify the specific `GtkFrame` widget (useful when using the `Gtk.Builder` and `glade`). Time is guaranteed to be positive.
"""
function timewidget(t1::Dates.Time; widget=nothing, signal=nothing)
    zerotime = Dates.Time(0,0,0) # convenient since we'll use it frequently
    b = Gtk.GtkBuilder(filename=joinpath(@__DIR__, "time.glade"))
    if signal == nothing
        signal = Signal(t1) # this is the input signal, we can push! into it to update the widget
    end
    S = map(signal) do x
        (Dates.Second(x), x) # crop the seconds from the Time signal, but keep the time for the next (minutes) crop
    end
    M = map(S) do x
        x = last(x) # this is the time
        (Dates.Minute(x), x) # crop the minutes out of this tuple signal, and again, keep hold of the time for the next (hour) crop
    end
    H = map(M) do x
        x = last(x)
        (Dates.Hour(x), x) # last crop, we have the hours now, and the time is kept as well
    end
    t2 = map(last, H) # here is the final time
    bind!(signal, t2) # we connect the input and output times so that any update to the resulting time will go into the input signal and actually show on the widgets
    Sint = Signal(Dates.value(first(value(S)))) # necessary for now, until range-like GtkReactive.widgets can accept other ranges.
    Ssb = spinbutton(-1:60, widget=b["second"], signal=Sint) # allow for values outside the actual range of seconds so that we'll be able to increase and decrease minutes.
    foreach(Sint) do x
        Δ = Dates.Second(x) - first(value(S)) # how much did we change by, this should always be ±1
        new_t = value(signal) + Δ # new time
        new_t = new_t < zerotime ? zerotime : new_t # julia Time is allowed negative values, here we correct for that
        new_x = Dates.Second(new_t) # new seconds
        push!(S, (new_x, new_t)) # update that specific widget, here the magic begins, this update will cascade down the widget-line...
    end
    Sint2 = map(src -> Dates.value(Dates.Second(src)), t2) # Any change in the value of the seconds, namely 60 -> 0, needs to loop back into the beginning of this last chain of events.
    Sint3 = droprepeats(Sint2) # important, otherwise we get an endless update loop
    bind!(Sint, Sint3, false) # final step of connecting the two
    # everything is the same for minutes:
    Mint = Signal(Dates.value(first(value(M))))
    Msb = spinbutton(-1:60, widget=b["minute"], signal=Mint)
    foreach(Mint) do x
        Δ = Dates.Minute(x) - first(value(M))
        new_t = value(signal) + Δ
        new_t = new_t < zerotime ? zerotime : new_t
        new_x = Dates.Minute(new_t)
        push!(M, (new_x, new_t))
    end
    Mint2 = map(src -> Dates.value(Dates.Minute(src)), t2)
    Mint3 = droprepeats(Mint2)
    bind!(Mint, Mint3, false)
    # while I think this next part is not entirely necessary for Hours, my brain hurts and I want this to be over. It works.
    Hint = Signal(Dates.value(first(value(H))))
    Hsb = spinbutton(0:23, widget=b["hour"], signal=Hint)
    foreach(Hint) do x
        Δ = Dates.Hour(x) - first(value(H))
        new_t = value(signal) + Δ
        new_t = new_t < zerotime ? zerotime : new_t
        new_x = Dates.Hour(new_t)
        push!(H, (new_x, new_t))
    end
    Hint2 = map(src -> Dates.value(Dates.Hour(src)), t2)
    Hint3 = droprepeats(Hint2)
    bind!(Hint, Hint3, false)

    if widget == nothing
        return TimeWidget(signal, b["frame"])
    else
        push!(widget, b["frame"])
        return TimeWidget(signal, widget)
    end
end

"""
    datetimewidget(datetime)

Return a datetime widget that includes the `DateTime` and a `GtkBox` with the
year, month, day, hour, minute, and second widgets in it. You can specify the
specific `SpinButton` widgets for the hour, minute, and second (useful when using
`Gtk.Builder` and `glade`). Date and time are guaranteed to be positive.
"""
function datetimewidget(t1::DateTime; widget=nothing, signal=nothing)
    zerotime = DateTime(0,1,1,0,0,0)
    b = Gtk.GtkBuilder(filename=joinpath(@__DIR__, "datetime.glade"))
    # the same logic is applied here as for `timewidget`
    if signal == nothing
        signal = Signal(t1)
    end
    S = map(signal) do x
        (Dates.Second(x), x)
    end
    M = map(S) do x
        x = last(x)
        (Dates.Minute(x), x)
    end
    H = map(M) do x
        x = last(x)
        (Dates.Hour(x), x)
    end
    d = map(H) do x
        x = last(x)
        (Dates.Day(x), x)
    end
    m = map(d) do x
        x = last(x)
        (Dates.Month(x), x)
    end
    y = map(m) do x
        x = last(x)
        (Dates.Year(x), x)
    end
    t2 = map(last, y)
    bind!(signal, t2)
    Sint = Signal(Dates.value(first(value(S))))
    Ssb = spinbutton(-1:60, widget=b["second"], signal=Sint)
    foreach(Sint) do x
        Δ = Dates.Second(x) - first(value(S))
        new_t = value(signal) + Δ
        new_t = new_t < zerotime ? zerotime : new_t
        new_x = Dates.Second(new_t)
        push!(S, (new_x, new_t))
    end
    Sint2 = map(src -> Dates.value(Dates.Second(src)), t2)
    Sint3 = droprepeats(Sint2)
    bind!(Sint, Sint3, false)
    Mint = Signal(Dates.value(first(value(M))))
    Msb = spinbutton(-1:60, widget=b["minute"], signal=Mint)
    foreach(Mint) do x
        Δ = Dates.Minute(x) - first(value(M))
        new_t = value(signal) + Δ
        new_t = new_t < zerotime ? zerotime : new_t
        new_x = Dates.Minute(new_t)
        push!(M, (new_x, new_t))
    end
    Mint2 = map(src -> Dates.value(Dates.Minute(src)), t2)
    Mint3 = droprepeats(Mint2)
    bind!(Mint, Mint3, false)
    Hint = Signal(Dates.value(first(value(H))))
    Hsb = spinbutton(-1:24, widget=b["hour"], signal=Hint)
    foreach(Hint) do x
        Δ = Dates.Hour(x) - first(value(H))
        new_t = value(signal) + Δ
        new_t = new_t < zerotime ? zerotime : new_t
        new_x = Dates.Hour(new_t)
        push!(H, (new_x, new_t))
    end
    Hint2 = map(src -> Dates.value(Dates.Hour(src)), t2)
    Hint3 = droprepeats(Hint2)
    bind!(Hint, Hint3, false)
    dint = Signal(Dates.value(first(value(d))))
    dsb = spinbutton(-1:32, widget=b["day"], signal=dint)
    foreach(dint) do x
        Δ = Dates.Day(x) - first(value(d))
        new_t = value(signal) + Δ
        new_t = new_t < zerotime ? zerotime : new_t
        new_x = Dates.Day(new_t)
        push!(d, (new_x, new_t))
    end
    dint2 = map(src -> Dates.value(Dates.Day(src)), t2)
    dint3 = droprepeats(dint2)
    bind!(dint, dint3, false)
    mint = Signal(Dates.value(first(value(m))))
    msb = spinbutton(-1:13, widget=b["month"], signal=mint)
    foreach(mint) do x
        Δ = Dates.Month(x) - first(value(m))
        new_t = value(signal) + Δ
        new_t = new_t < zerotime ? zerotime : new_t
        new_x = Dates.Month(new_t)
        push!(m, (new_x, new_t))
    end
    mint2 = map(src -> Dates.value(Dates.Month(src)), t2)
    mint3 = droprepeats(mint2)
    bind!(mint, mint3, false)
    yint = Signal(Dates.value(first(value(y))))
    ysb = spinbutton(-1:10000, widget=b["year"], signal=yint)
    foreach(yint) do x
        Δ = Dates.Year(x) - first(value(y))
        new_t = value(signal) + Δ
        new_t = new_t < zerotime ? zerotime : new_t
        new_x = Dates.Year(new_t)
        push!(y, (new_x, new_t))
    end
    yint2 = map(src -> Dates.value(Dates.Year(src)), t2)
    yint3 = droprepeats(yint2)
    bind!(yint, yint3, false)

    if widget == nothing
        return TimeWidget(signal, b["frame"])
    else
        push!(widget, b["frame"])
        return TimeWidget(signal, widget)
    end
end


text:
struct ZeroPolynomialSet{DT <: AbstractSemialgebraicSet,
                         BT <: MB.AbstractPolynomialBasis,
                         MT <: AbstractMonomial,
                         MVT <: AbstractVector{MT}} <: MOI.AbstractVectorSet
    domain::DT
    basis::Type{BT}
    monomials::MVT
end
MOI.dimension(set::ZeroPolynomialSet) = length(set.basis)
Base.copy(set::ZeroPolynomialSet) = set

# `x`-in-`PlusMinusSet(set)` iff `x`-in-`set` and `-x`-in-`set`.
struct PlusMinusSet{ST <: MOI.AbstractVectorSet} <: MOI.AbstractVectorSet
    set::ST
end
MOI.dimension(set::PlusMinusSet) = MOI.dimension(set.set)
Base.copy(set::PlusMinusSet) = set


text:
using Base: splat

test = """
    .#.
    ..#
    ###
    """

puzzle = """
    #.#####.
    #..##...
    .##..#..
    #.##.###
    .#.#.#..
    #.##..#.
    #####..#
    ..#.#.##
    """

Base.@kwdef mutable struct PocketDim{d}
    space::Dict{CartesianIndex{d}, Bool} = Dict()
    min::CartesianIndex{d} = zero(CartesianIndex{d})
    max::CartesianIndex{d} = zero(CartesianIndex{d})
end

Base.keys(pd::PocketDim) = pd.min:pd.max
Base.values(pd::PocketDim) = pd[keys(pd)]
Base.getindex(pd::PocketDim, i::CartesianIndices) = getindex.(Ref(pd), i)
Base.getindex(pd::PocketDim, i::AbstractVector) = getindex.(Ref(pd), i)
Base.getindex(pd::PocketDim, i::CartesianIndex) = get(pd.space, i, false)

function Base.getindex(pd::PocketDim, i...)
    vs = pd[CartesianIndices(Tuple(minimum(j):maximum(j) for j in i))]
    return any(j -> j isa AbstractRange, i) ? vs : only(vs)
end

function Base.setindex!(pd::PocketDim, v, i...)
    pd.space[CartesianIndex(i...)] = v
    pd.min = min(pd.min, CartesianIndex(i...))
    pd.max = max(pd.max, CartesianIndex(i...))
end

function Base.show(io::IO, pd::PocketDim)
    for z = CartesianIndex(pd.min.I[3:end]):CartesianIndex(pd.max.I[3:end])
        println(io, join(map((n, v) -> "$n=$v", "zw", z.I), ','))
        for r in pd.min[1]:pd.max[1]
            row = pd[r, pd.min[2]:pd.max[2], z.I...]
            println(io, join(map(v -> v ? '#' : '.', row), ' '))
        end
        println(io)
    end
end

function parse_input(input::AbstractString; d::Int = 3)
    lines = split(input, r"\s+"; keepempty = false)
    plane = map.(==("#"), split.(lines, "")) |> splat(hcat) |> permutedims
    pd = PocketDim{d}()
    for k in keys(plane)
        pd[k - one(k), fill(0, d - 2)...] = plane[k]
    end
    return pd
end

function neighbors(ci::CartesianIndex)
    return setdiff((ci - one(ci)):(ci + one(ci)), [ci])
end

function step!(pd::PocketDim)::PocketDim
    change = []
    for i in (pd.min - one(pd.min)):(pd.max + one(pd.max))
        cube = pd[i]
        nbs = pd[neighbors(i)] |> sum
        if cube && !(2 <= nbs <= 3)
            push!(change, i => false)
        elseif !cube && nbs == 3
            push!(change, i => true)
        end
    end
    for (i, v) in change
        pd[i] = v
    end
    return pd
end

function run!(pd::PocketDim, cycles::Int = 6)::Int
    for _ in 1:cycles
        step!(pd)
    end
    return sum(values(pd))
end

q1(input) = run!(parse_input(input; d = 3))
q2(input) = run!(parse_input(input; d = 4))

@assert q1(test) == 112
@assert q2(test) == 848

println("Q1: ", q1(puzzle))
println("Q2: ", q2(puzzle))

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
module TestSQS
using AWS
using AWS.SQS
using Base.Test

function check_resp(resp, resp_obj_type=Any)
    info("response: ", resp.obj)
    @test 200 <= resp.http_code <= 206
    @test isa(resp, AWS.SQS.SQSResponse)
    @test isa(resp.obj, resp_obj_type)
    @test !isa(resp.obj, SQSError)
end

function runtests(env, config)
    awsAccountID = get(config, "awsAccountID", "")

    if isempty(awsAccountID)
        warn("skipping SQS tests, awsAccountID not set")
        return
    end

    queueName = "MyTest1"

    println("list all queues")
    resp = ListQueues(env)
    check_resp(resp, ListQueuesResult)
    @test isa(resp.obj.queueUrl, Union{Void,Vector{String}})

    println("create queue")
    attributes = Dict(
        "DelaySeconds" => "300",
        "VisibilityTimeout" => "120"
    )
    resp = CreateQueue(env; queueName=queueName, attribute=attributes)
    check_resp(resp, CreateQueueResult)
    @test endswith(resp.obj.queueUrl, queueName)

    qurl = resp.obj.queueUrl
    @show qurl
    @test qurl == GetQueueUrl(env; queueName=queueName).obj.queueUrl

    println("get queue attributes")
    resp = GetQueueAttributes(env; queueUrl=qurl, attributeName=["All"])
    check_resp(resp, AWS.SQS.GetQueueAttributesResult)
    @test isa(resp.obj.attribute, Dict{String,String})
    @test !isempty(resp.obj.attribute)

    println("send message")
    msgAttributes = Dict(
        "some-attribute"  => MessageAttributeValue(; dataType="Number.integer-ish", stringValue="0"),
        "other-attribute" => MessageAttributeValue(; dataType="String.yyy",         stringValue="My yyy string"),
        "bin-attribute"   => MessageAttributeValue(; dataType="Binary.jpg",         binaryValue=[0x0,0x1,0x2,0x3,0x4,0x5,0x6])
    )
    resp = SendMessage(env; queueUrl=qurl, delaySeconds=0, messageBody="test", messageAttribute=msgAttributes)
    check_resp(resp, SendMessageResult)
    msgid = resp.obj.messageId

    println("receive messages")
    resp = ReceiveMessage(env; queueUrl=qurl, attributeName=["All"], messageAttributeName=["All"])
    check_resp(resp, ReceiveMessageResult)
    @test isa(resp.obj.message, Vector{Message})
    info("received ", length(resp.obj.message), " messages")

    msgids = [x.messageId for x in resp.obj.message]
    @test msgid in msgids
    msgbodies = [x.body for x in resp.obj.message]
    @test "test" in msgbodies

    println("delete messages")
    for hndl in [x.receiptHandle for x in resp.obj.message]
        resp = DeleteMessage(env; queueUrl=qurl, receiptHandle=hndl)
        check_resp(resp)
        info("deleted ", hndl)
    end

    println("send message batch")
    id1 = "SendMessageBatchMessage1"
    delaySeconds1 = 10
    msgAttributes1 = Dict(
        "some-attribute"  => MessageAttributeValue(; dataType="Number.integer-ish", stringValue="111"),
        "other-attribute" => MessageAttributeValue(; dataType="String.yyy",         stringValue="My 111 string"),
        "bin-attribute"   => MessageAttributeValue(; dataType="Binary.jpg",         binaryValue=[0x0,0x1,0x2,0x3,0x4,0x5,0x6])
    )
    messageBody1 = "Message Body 1"
    sendMessageBatchRequestEntry1 = SendMessageBatchRequestEntry(; delaySeconds=delaySeconds1, id=id1, messageAttribute=msgAttributes1, messageBody=messageBody1)

    id2 = "SendMessageBatchMessage2"
    delaySeconds2 = 10
    msgAttributes2 = Dict(
        "some-attribute"  => MessageAttributeValue(; dataType="Number.integer-ish", stringValue="222"),
        "other-attribute" => MessageAttributeValue(; dataType="String.yyy",         stringValue="My 222 string"),
        "bin-attribute"   => MessageAttributeValue(; dataType="Binary.jpg",         binaryValue=[0x0,0x1,0x2,0x3,0x4,0x5,0x6])
    )
    messageBody2 = "Message Body 2"
    sendMessageBatchRequestEntry2 = SendMessageBatchRequestEntry(; delaySeconds=delaySeconds2, id=id2, messageAttribute=msgAttributes2, messageBody=messageBody2)

    id3 = "SendMessageBatchMessage3"
    delaySeconds3 = 10
    messageBody3 = "Message Body 3"
    sendMessageBatchRequestEntry3 = SendMessageBatchRequestEntry(; delaySeconds=delaySeconds3, id=id3, messageBody=messageBody3)

    id4 = "SendMessageBatchMessage4"
    delaySeconds4 = 10
    messageBody4 = "Message Body 4"
    sendMessageBatchRequestEntry4 = SendMessageBatchRequestEntry(; delaySeconds=delaySeconds4, id=id4, messageBody=messageBody4)

    sendMessageBatchRequestEntry = [
        sendMessageBatchRequestEntry1,
        sendMessageBatchRequestEntry2,
        sendMessageBatchRequestEntry3,
        sendMessageBatchRequestEntry4
    ]

    req = SendMessageBatchRequest(; sendMessageBatchRequestEntry=sendMessageBatchRequestEntry, queueUrl=qurl)
    resp = SendMessageBatch(env, req)
    check_resp(resp)

    println("Sleeping for 20 secs !!!")
    sleep(20)

    println("change message visibility")
    resp = ReceiveMessage(env; queueUrl=qurl)
    check_resp(resp, ReceiveMessageResult)
    msg = resp.obj.message[1]

    req = ChangeMessageVisibilityRequest(; queueUrl=qurl, receiptHandle=msg.receiptHandle, visibilityTimeout=30)
    resp = ChangeMessageVisibility(env, req)
    check_resp(resp)

    println("change message visibility batch")
    resp1 = ReceiveMessage(env; queueUrl=qurl)
    msg1 = resp1.obj.message[1]
    id1 = msg1.messageId

    resp2 = ReceiveMessage(env; queueUrl=qurl)
    msg2 = resp2.obj.message[1]
    id2 = msg2.messageId

    changeMessageVisibilityBatchRequestEntry1 = ChangeMessageVisibilityBatchRequestEntry(; id=id1, receiptHandle=msg1.receiptHandle, visibilityTimeout=60)
    changeMessageVisibilityBatchRequestEntry2 = ChangeMessageVisibilityBatchRequestEntry(; id=id2, receiptHandle=msg2.receiptHandle, visibilityTimeout=90)
    changeMessageVisibilityBatchRequestEntry = [
        changeMessageVisibilityBatchRequestEntry1,
        changeMessageVisibilityBatchRequestEntry2
    ]

    req = ChangeMessageVisibilityBatchRequest(; queueUrl=qurl, changeMessageVisibilityBatchRequestEntry=changeMessageVisibilityBatchRequestEntry)
    resp = ChangeMessageVisibilityBatch(env, req)
    check_resp(resp)

    println("delete message batch")
    deleteMessageBatchRequestEntry1 = DeleteMessageBatchRequestEntry(; id=id1, receiptHandle=msg1.receiptHandle)
    deleteMessageBatchRequestEntry2 = DeleteMessageBatchRequestEntry(; id=id2, receiptHandle=msg2.receiptHandle)
    deleteMessageBatchRequestEntry = [
        deleteMessageBatchRequestEntry1,
        deleteMessageBatchRequestEntry2
    ]

    req = DeleteMessageBatchRequest(; deleteMessageBatchRequestEntry=deleteMessageBatchRequestEntry, queueUrl=qurl)
    resp = DeleteMessageBatch(env, req)
    check_resp(resp)

    println("list dead letter source queues")
    resp = ListDeadLetterSourceQueues(env; queueUrl=qurl)
    check_resp(resp)

    println("purge queue")
    resp = PurgeQueue(env; queueUrl=qurl)
    check_resp(resp)

    println("add permission")
    awsAccountIds = [awsAccountID, awsAccountID]
    actionNames = ["SendMessage", "ReceiveMessage"]

    resp = AddPermission(env; queueUrl=qurl, label="My Permission 1", aWSAccountId=awsAccountIds, actionName=actionNames)
    check_resp(resp)

    println("remove  permission")
    resp = RemovePermission(env; queueUrl=qurl, label="My Permission 1")
    check_resp(resp)

    println("testing delete queue")
    resp = DeleteQueue(env; queueUrl=qurl)
    check_resp(resp)
end

end # module TestSQS


text:
using PowerDynamics
using ForwardDiff
using Distributions:Uniform
import PowerDynamics:rhs
using DifferentialEquations
using LinearAlgebra
using Optim
#using DiffEqCallbacks
include("../PDpatches.jl")
include("system.jl")
include("PowerFlow.jl")

"""
# neuer text?
A method to find random pertubations in the case of systems with algebraic
constraints. There is then a set of valid initial conditions which lies on a so
called constrained manifold.
The system must still lie on this manifold after the random pertubation.
This could be applied to Basin Stability Methods in the case of complexer
Generator Methods which contain algebraic constraints.
"""
function RandPertWithConstrains(pg::PowerGrid, operationpoint, node, interval_V, interval_θ,Q; ProjectionMethod = "Optim")
    #Pi = pg.nodes[node].P # Power of the specific node
    g(θ,i,u) = Q - 3 * i * exp(-1im * θ) * u * exp(-1im * θ) # Constraint, Power Flow Analysis
    z = [operationpoint[node,:θ], operationpoint[node,:iabs],operationpoint[node,:v]]

    #Jacg_imag = ForwardDiff.jacobian(x -> imag(g(x[1],x[2],x[3])), z)
    #Jacg_real = ForwardDiff.jacobian(x -> real(g(x[1],x[2],x[3])), z)

    #Jacg = Jacg_real + 1im .* Jacg_imag  # I could not find a package which is able to handle complex valued function
    Jacg = ForwardDiff.gradient(x -> real(g(x[1],x[2],x[3])), z) .+ 1im .* ForwardDiff.gradient(x -> imag(g(x[1],x[2],x[3])), z)

    kerJacg = nullspace(Matrix(Jacg')) # kernel of the Jacobian, returns orthogonal vectors

    num_basis_vec = size(kerJacg,2)

    Frand = [rand(Uniform(interval_θ[1],interval_θ[2])), rand(Uniform(interval_V[1],interval_V[2]))]

    dz = kerJacg * Frand

    if ProjectionMethod == "Optim"
        p_dash = ProjOptim(dz, g) # Calculates the projection onto M using Optim

    elseif ProjectionMethod == "ODE"
        p_dash = ProjODE(g, node, operationpoint,dz, Q)

    else
        ArgumentError("Please use either Optim or ODE as the ProjectionMethod.")
    end

    PerturbedState = copy(operationpoint) # otherwise operationpoint is overwritten
    PerturbedState[2,:θ] = abs(p_dash[1])
    PerturbedState[2,:u] = p_dash[3]
    return PerturbedState, p_dash
end

"""
Calculates and returns the orthogonal Projection of a vector w on the vector v.
"""
function Projvw(v,w)
    dot(w,v) * v * (1 / norm(v)^2)
end

"""
Calculates the Projection from the Tangetial Room TzM to the the Manifold M by
minimizing the function:
    f(x) = a * |dz-p''| + b * |g(p')|.
Where g(p') is the Constraint at point p'.
Inputs:
        dz: The calculated pertubation inside of TzM
        g: The constraint function
        Q: Reactive Power at the perturbed node
Outputs:
        p_dash: Solution of the ODE problem, a Pertubation on M
"""
function ProjOptim(dz, g::Function; a = 1.0, b = 100.0)
    Proj(x) = a * norm(dz - [x[1] + 1im * x[2], x[3] + 1im * x[4], x[5] + 1im * x[6]]) +
              b * abs(g(x[1] + 1im * x[2], x[3] + 1im * x[4], x[5] + 1im * x[6]))

    x0 = [real(dz[1]), imag(dz[1]),
          real(dz[2]), imag(dz[2]),
          real(dz[3]), imag(dz[3])]

    proj = optimize(Proj,x0)

    minimum = proj.minimizer
    p_dash = [minimum[1] + 1im * minimum[2],
              minimum[3] + 1im * minimum[4],
              minimum[5] + 1im * minimum[6]]

    return p_dash
end

"""
Calculates the Projection from the Tangetial Room TzM to the the Manifold M by
solving a ODE of the form:
    f(x) = y' = - g(dz, y).
Where dz is the result of the pertubation and a point in TzM.
Which could result in a fixed point g(x_v, y*) = 0.
Inputs:
        dz: The calculated pertubation inside of TzM
        g: The constraint function
        Q: Reactive Power at the perturbed node
        endtime: Integration time of the problem
Outputs:
        p_dash: The result of optimization, a Pertubation on M
"""
function ProjODE(g::Function,node,op,dz,Q;endtime = 10.0)
    tspan = (0.0,endtime)

    u0 = [dz[1], dz[2], dz[3], g(dz[1], dz[2], dz[3]), op[node,:ω]]
    vars = powergrid.nodes[node]
    #function constraint(du,u,p,t)
    #    du[4] = - (Q - 3 * u[3] * u[2] * exp(-2im * u[1]))
    #end
    function constraint(du,u,p,t)
        Ω_H = vars.Ω / (2 * vars.H)
        i_c = 1im * u[2] * exp(-1im * u[1])
        e_c = 1im * u[3] * exp(-1im * u[1])
        p = real(u[3] * conj(u[2]))
        e_d = 0                         # simplification from K.Schmietendorf paper
        e_q = imag(e_c)
        i_d = real(i_c)
        i_q = imag(i_c)

        du[1] = u[5]
        du[5] = (vars.P - vars.D * u[5] - p - (vars.X_q_dash - vars.X_d_dash) * i_d * i_q) * Ω_H
        de_q = (1 / vars.T_d_dash) * (vars.E_f - e_q + i_d * (vars.X_d - vars.X_d_dash))
        # -> u = e_q * exp(1im * θ)
        du[3] = de_q * exp(1im * u[1]) + u[3] * 1im * u[5]
        du[4] = u[4] - (Q - 3 * u[3] * u[2] * exp(-2im * u[1]))
    end

    prob = ODEProblem(constraint,u0,tspan)

    sol = solve(prob, Rosenbrock23(autodiff=false), abstol=1e-8,reltol=1e-6,tspan=Inf, force_dtmin = true)
    p_dash = sol[1:3,end]
end

function rhs(pg::PowerGrid)
    network_dynamics(map(construct_vertex, pg.nodes), map(construct_edge, pg.lines), pg.graph)
end

"""
A method to assures that the conditions given from g stay on the contrained
manifold. A DifferentialEquations Callback is used.
This DOES NOT work because of a type missmatch regarding complex numbers
Inputs:
    pg: Power grid, a graph containg nodes and lines
    endtime: Endtime of the simulation
    u0: current state of the powergrid pg
Outputs:
    solution: a PowerGrid solution containg the timeseries and the Powergrid pg
"""
function SimManifoldCd(pg::PowerGrid, Q, u0::State, endtime)
    ode = ODEProblem(rhs(pg), u0.vec, (0.0, endtime), nothing)

    function g(resid, x, e_s, e_d, p, t)
            i = total_current(e_s, e_d) + Y_n * (x[1] + x[2] * im)
            u = x[1] + x[2] * im
            θ = x[3]
            ω = x[4]
            resid[1] = abs(Q) - 3 * abs(i * u * exp((-2im) * θ))
    end

    function g1(resid, x, e_s, e_d, p, t)
        resid[1] = x[2]^2 + x[1]^2 - 2
        resid[2] = 0
    end

    cb = ManifoldProjection(g1)

    dqsol = solve(ode, Rodas5(autodiff = false), save_everystep=false, callback=cb)
    solution = PowerGridSolution(dqsol, pg)
    return solution
end

max_stars_repo_name: luap-pik/PowerDynIC

max_issues_repo_name: luap-pik/PowerDynIC

max_forks_repo_name: luap-pik/PowerDynIC


text:
abstract PersistentMap{K, V} <: Associative{K, V}

type NotFound end

immutable KVPair{K, V}
    key::K
    value::V
end

Base.convert(::Type{Tuple}, kv::KVPair) = (kv.key, kv.value)
Base.convert{K, V}(::Type{KVPair{K, V}}, kv::KVPair) =
    KVPair{K, V}(convert(K, kv.key), convert(V, kv.value))

Base.isequal(kv1::KVPair, kv2::KVPair) =
    isequal(kv1.key, kv2.key) && isequal(kv1.value, kv2.value)
Base.isequal(kv::KVPair, tup::Tuple) = isequal(convert(Tuple, kv), tup)
Base.isequal(tup::Tuple, kv::KVPair) = isequal(kv, tup)
==(kv1::KVPair, kv2::KVPair) = kv1.key == kv2.key && kv1.value == kv2.value
==(kv::KVPair, tup::Tuple) = convert(Tuple, kv) == tup
==(tup::Tuple, kv::KVPair) = kv == tup

Base.hash(kv::KVPair) = @compat UInt(Base.hash(kv.key, hash(kv.value)))

@compat Base.show(io::IO, ::MIME"text/plain", kv::KVPair) = print(io, "$(kv.key) => $(kv.value)")


text:
module EntropyHub

export
    ApEn,
    AttnEn,
    BubbEn,
    CoSiEn,
    CondEn,
    DispEn,
    DistEn,
    EnofEn,
    FuzzEn,
    GridEn,
    IncrEn,
    K2En,
    PermEn,
    PhasEn,
    SampEn,
    SlopEn,
    SyDyEn,
    SpecEn,

    XApEn,
    XCondEn,
    XDistEn,
    XFuzzEn,
    XK2En,
    XPermEn,
    XSpecEn,
    XSampEn,

    SampEn2D,
    FuzzEn2D,
    DistEn2D,
    DispEn2D,
    PermEn2D,
    EspEn2D,

    MSobject,
    MSEn,
    cMSEn,
    rMSEn,
    hMSEn,
    XMSEn,
    cXMSEn,
    rXMSEn,
    hXMSEn,

    EMD,
    ExampleData

# using Reexport

include("./_ApEn.jl"),
include("./_AttnEn.jl"),
include("./_BubbEn.jl"),
include("./_CondEn.jl"),
include("./_CoSiEn.jl"),
include("./_DispEn.jl"),
include("./_DistEn.jl"),
include("./_FuzzEn.jl"),
include("./_EnofEn.jl"),
include("./_GridEn.jl"),
include("./_IncrEn.jl"),
include("./_K2En.jl"),
include("./_PermEn.jl"),
include("./_PhasEn.jl"),
include("./_SampEn.jl"),
include("./_SpecEn.jl"),
include("./_SyDyEn.jl"),
include("./_SlopEn.jl"),

include("./_XApEn.jl"),
include("./_XFuzzEn.jl"),
include("./_XDistEn.jl"),
include("./_XCondEn.jl"),
include("./_XSpecEn.jl"),
include("./_XSampEn.jl"),
include("./_XK2En.jl"),
include("./_XPermEn.jl"),

include("./_MSobject.jl"),
include("./_ExampleData.jl"),

include("./_MSEn.jl"),
include("./_XMSEn.jl"),
include("./_cMSEn.jl"),
include("./_cXMSEn.jl"),
include("./_rMSEn.jl"),
include("./_rXMSEn.jl"),
include("./_hMSEn.jl"),
include("./_hXMSEn.jl"),

include("./_SampEn2D.jl"),
include("./_FuzzEn2D.jl"),
include("./_DistEn2D.jl"),
include("./_DispEn2D.jl"),
include("./_PermEn2D.jl"),
include("./_EspEn2D.jl"),

#=
@reexport using ._ApEn
@reexport using ._MSobject

@reexport using ._AttnEn
@reexport using ._BubbEn
@reexport using ._CondEn
@reexport using ._CoSiEn
@reexport using ._DistEn
@reexport using ._DispEn
@reexport using ._EnofEn
@reexport using ._FuzzEn =#

# Base Entropies:
using ._ApEn: ApEn
using ._AttnEn: AttnEn
using ._BubbEn: BubbEn
using ._CoSiEn: CoSiEn
using ._CondEn: CondEn
using ._DispEn: DispEn
using ._DistEn: DistEn
using ._EnofEn: EnofEn
using ._FuzzEn: FuzzEn
using ._GridEn: GridEn
using ._IncrEn: IncrEn
using ._K2En: K2En
using ._PermEn: PermEn
using ._PhasEn: PhasEn
using ._SampEn: SampEn
using ._SlopEn: SlopEn
using ._SpecEn: SpecEn
using ._SyDyEn: SyDyEn

# Cross Entropies:
using ._XApEn: XApEn
using ._XCondEn: XCondEn
using ._XDistEn: XDistEn
using ._XFuzzEn: XFuzzEn
using ._XK2En: XK2En
using ._XPermEn: XPermEn
using ._XSampEn: XSampEn
using ._XSpecEn: XSpecEn

# Bidimensional Entropies
using ._SampEn2D: SampEn2D
using ._DistEn2D: DistEn2D
using ._FuzzEn2D: FuzzEn2D
using ._DispEn2D: DispEn2D
using ._PermEn2D: PermEn2D
using ._EspEn2D: EspEn2D

# (cross) Multiscale Entropies
using ._MSobject: MSobject

using ._MSEn: MSEn, EMD
using ._cMSEn: cMSEn
using ._rMSEn: rMSEn
using ._hMSEn: hMSEn

using ._XMSEn: XMSEn
using ._cXMSEn: cXMSEn
using ._rXMSEn: rXMSEn
using ._hXMSEn: hXMSEn

using ._ExampleData: ExampleData

greet() = print(raw"""
     ___  _   _  _____  _____  ____  ____  _     _
    |  _|| \ | ||_   _||     \|    ||    || \   / |   ___________
    | \_ |  \| |  | |  |   __/|    ||  __| \ \_/ /   /  _______  \
    |  _|| \ \ |  | |  |   \  |    || |     \   /   |  /  ___  \  |
    | \_ | |\  |  | |  | |\ \ |    || |      | |    | |  /   \  | |
    |___||_| \_|  |_|  |_| \_||____||_|      |_|   _|_|__\___/  | |
     _   _  _   _  ____                           / |__\______\/  |
    | | | || | | ||    \     An open-source      |  /\______\__|_/
    | |_| || | | ||    |     toolkit for         | |  /   \  | |
    |  _  || | | ||    \     entropic time-      | |  \___/  | |
    | | | || |_| ||     \    series analysis     |  \_______/  |
    |_| |_|\_____/|_____/                         \___________/
    """)

raw"""
 ___  _   _  _____  _____  ____  ____  _     _
|  _|| \ | ||_   _||     \|    ||    || \   / |   ___________
| \_ |  \| |  | |  |   __/|    ||  __| \ \_/ /   /  _______  \
|  _|| \ \ |  | |  |   \  |    || |     \   /   |  /  ___  \  |
| \_ | |\  |  | |  | |\ \ |    || |      | |    | |  /   \  | |
|___||_| \_|  |_|  |_| \_||____||_|      |_|   _|_|__\___/  | |
 _   _  _   _  ____                           / |__\______\/  |
| | | || | | ||    \     An open-source      |  /\______\__|_/
| |_| || | | ||    |     toolkit for         | |  /   \  | |
|  _  || | | ||    \     entropic time-      | |  \___/  | |
| | | || |_| ||     \    series analysis     |  \_______/  |
|_| |_|\_____/|_____/                         \___________/

EntropyHub functions belong to one of five main classes/categories:
    Base Entropies             >>  e.g. Approximate Entropy (ApEn),
                                        Sample Entropy (SampEn)
    Cross Entropies            >>  e.g. Cross-Approximate Entropy (XApEn)
                                        Cross-Sample Entropy (XSampEn)
    Bidimensional Entropies    >>  e.g. Bidimensional Sample Entropy (SampEn2D)
                                        Bidimensional Fuzzy Entropy (FuzzEn2D)
    Multiscale Entropies       >>  e.g. Multiscale Sample Entropy (MSEn)
                                        Refined Multiscale Sample Entropy (rMSEn)
                                        Composite Multiscale Sample Entropy (cMSEn)
    Multiscale Cross Entropies >>  e.g. Multiscale Cross-Sample Entropy (XMSEn)
                                        Refined Multiscale Cross-Sample Entropy (rXMSEn)

_________________________________________________________________________
Base Entropies                                        |	Function Name
______________________________________________________|__________________
Approximate Entropy                               	  |	ApEn
Sample Entropy                                		    |	SampEn
Fuzzy Entropy                                 		    |	FuzzEn
Kolmogorov Entropy                            		    | K2En
Permutation Entropy                           		    |	PermEn
Conditional Entropy                           		    |	CondEn
Distribution Entropy                          		    |	DistEn
Spectral Entropy                              		    |	SpecEn
Dispersion Entropy                            		    |	DispEn
Symbolic Dynamic Entropy                          	  |	SyDyEn
Increment Entropy                                 	  |	IncrEn
Cosine Similarity Entropy                         	  |	CoSiEn
Phase Entropy                                         |	PhasEn
Slope Entropy                                      	  |	SlopEn
Bubble Entropy                                		    |	BubbEn
Gridded Distribution Entropy                          |	GridEn
Entropy of Entropy                            		    |	EnofEn
Attention Entropy                                     |	AttnEn

_________________________________________________________________________
Cross Entropies                                       |	Function Name
______________________________________________________|__________________
Cross Sample Entropy                                  |	XSampEn
Cross Approximate Entropy                             |	XApEn
Cross Fuzzy Entropy                                   |	XFuzzEn
Cross Permutation Entropy                             |	XPermEn
Cross Conditional Entropy                             |	XCondEn
Cross Distribution Entropy                            |	XDistEn
Cross Spectral Entropy                                |	XSpecEn
Cross Kolmogorov Entropy                              |	XK2En

_________________________________________________________________________
Bidimensional Entropies                              |	Function Name
_____________________________________________________|__________________
Bidimensional Sample Entropy                         |	SampEn2D
Bidimensional Fuzzy Entropy                          |	FuzzEn2D
Bidimensional Distribution Entropy                   |	DistEn2D
Bidimensional Dispersion Entropy                     |	DispEn2D
Bidimensional Permutation Entropy                    |	PermEn2D
Bidimensional Espinosa Entropy                       |	EspEn2D

_________________________________________________________________________
Multiscale Entropy Functions                          | Function Name
______________________________________________________|__________________
Multiscale Entropy Object                             |   MSobject
                                                      |
Multiscale Entropy                                    |   MSEn
Composite/Refined-Composite Multiscale Entropy        |   cMSEn
Refined Multiscale Entropy                            |   rMSEn
Hierarchical Multiscale Entropy Object                |   hMSEn

_________________________________________________________________________
Multiscale Entropies	MSEn                            |	Function Name
_________________________________________________________________________
Multiscale Sample Entropy                             |
Multiscale Approximate Entropy                        |
Multiscale Fuzzy Entropy                              |
Multiscale Permutation Entropy                        |
Multiscale Dispersion Entropy                         |
Multiscale Cosine Similarity Entropy                  |
Multiscale Symblic Dynamic Entropy                    |	  MSobject
Multiscale Conditional Entropy                        |	     +
Multiscale Entropy of Entropy                         |  MSEn / cMSEn
Multiscale Gridded Distribution Entropy               |	rMSEn / hMSEn
Multiscale Slope Entropy                              |
Multiscale Phase Entropy                              |
Multiscale Kolmogorov Entropy                         |
Multiscale Distribution Entropy                       |
Multiscale Bubble Entropy                             |
Multiscale Increment Entropy			                    |
Multiscale Attention Entropy                          |

_________________________________________________________________________
Multiscale Cross-Entropy Functions                    |   Function Name
______________________________________________________|__________________
Multiscale Cross-Entropy Object                       |   MSobject
                                                      |
Multiscale Cross-Entropy                              |   XMSEn
Composite/Refined-Composite Multiscale Cross-Entropy  |   cXMSEn
Refined Multiscale Entropy                            |   rXMSEn
Hierarchical Multiscale Entropy Object                |   hXMSEn

_________________________________________________________________________
Multiscale Cross-Entropies                            |	Function Name
_________________________________________________________________________
Multiscale Cross-Sample Entropy                       |
Multiscale Cross-Approximate Entropy                  |
Multiscale Cross-Fuzzy Entropy                        |	MSobject
Multiscale Cross-Permutation Entropy                  |	    +
Multiscale Cross-Distribution Entropy                 |	XMSEn / cXMSEn
Multiscale Cross-Kolmogorov Entropy                   |	rXMSEn / hXMSEn
Multiscale Cross-Conditional Entropy                  |

  We kindly ask that if you use EntropyHub in your research, to please
  include the following citation with the appropriate version number,
  as well as original articles upon which functions are derived:

  Matthew W. Flood and Bernd Grimm (2021),
  "EntropyHub - An open source toolkit for entropic time series analysis"
  PLoS ONE 16(11):e0295448,
  DOI:  10.1371/journal.pone.0259448
  https://www.EntropyHub.xyz

  © Copyright 2021 Matthew W. Flood, EntropyHub

  Licensed under the Apache License, Version 2.0 (the "License");
  you may not use this file except in compliance with the License.
  You may obtain a copy of the License at

       http://www.apache.org/licenses/LICENSE-2.0

  Unless required by applicable law or agreed to in writing, software
  distributed under the License is distributed on an "AS IS" BASIS,
  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  See the License for the specific language governing permissions and
  limitations under the License.

  For Terms of Use see https://github.com/MattWillFlood/EntropyHub
"""

  end


text:
"""
	RemoteWebDriver(capabilities::Capabilities;
                    host::AbstractString = "localhost",
                    port::Integer = 4444,
                    path::AbstractString = "",
					kwargs...)::RemoteWebDriver

Specifies a remote web driver according to [W3C](https://www.w3.org/TR/webdriver).

```jldoctest
julia> capabilities = Capabilities("chrome")
Remote WebDriver Capabilities
browserName: chrome
julia> wd = RemoteWebDriver(capabilities, host = ENV["WEBDRIVER_HOST"], port = parse(Int, ENV["WEBDRIVER_PORT"]))
Remote WebDriver
julia> status(wd) # Ready to accept new sessions?
true
```
"""
struct RemoteWebDriver{C<:Capabilities}
    addr::String
    capabilities::C
    kw::Dict{String,Union{Int,String}}
    function RemoteWebDriver(
        capabilities::Capabilities;
        host::AbstractString = "localhost",
        port::Integer = 4444,
        path::AbstractString = "",
        kwargs...,
    )::RemoteWebDriver
        addr = URI(scheme = "http", host = host, port = port, path = path)
        isvalid(addr) || throw(ArgumentError("$addr is invalid"))
        addr = string(addr)
        kw = Dict{String,Union{Int,String}}(
            string(k) => isa(v, Bool) ? string(v) : v for (k, v) ∈ kwargs if !isnothing(v)
        )
        new{typeof(capabilities)}(addr, capabilities, kw)
    end
end
summary(io::IO, obj::RemoteWebDriver) = println(io, "Remote WebDriver")
function show(io::IO, obj::RemoteWebDriver)
    print(io, summary(obj))
end


text:
module TestAtBackslash

using AtBackslash: @\, AtBackslash
using Documenter: doctest
using Test: @test, @testset

@testset "All bound" begin
    @test (x = 1, y = 2) |> @\(_..., a = :x + :y, y = 10) === (x = 1, y = 10, a = 3)
    @test (x = 1, y = 2) |> @\(:y, :x) === (y = 2, x = 1)
    @test (x = 1, y = 2) |> @\(:x < :y) === true
    @test (x = 1, y = 2) |> @\(:x < :y < 3) === true
    @test (x = 1, y = 2) |> @\(:x + 10) == 11
    @test (x = 1, y = 2) |> @\(; :x) === (x = 1,)
    @test (x = 1, y = 2) |> @\[:x, :y] == [1, 2]
    @test (x = 1, y = 2) |> @\[:x :y] == [1 2]
    @test (x = 1, y = 2) |> @\[:x :y; :y :x] == [1 2; 2 1]
    @test (x = 1, y = 2) |> @\[_] == [(x = 1, y = 2)]
    @test (x = 1, y = 2) |> @\[_ _] == [(x = 1, y = 2) (x = 1, y = 2)]
    @test (x = 1, y = 2) |> @\[(; :x, :y)] == [(x = 1, y = 2)]
    @test (x = 1, y = 2) |> @\[(; x = :x, y = :y)] == [(x = 1, y = 2)]
    @test (x = 1, y = 2, name = :x) |> @\(_[:name]) == 1
    @test (x = 1, y = 2, name = :x) |> @\(; :name => _[:name] + 10) === (x = 11,)
    @test (x = 1, y = 2, name = :x) |> @\((; :name => _[:name] + 10)) === (x = 11,)
    @test (x = 1, y = 2, name = :x) |> @\(; :name => _[:name] + 10, :y) === (x = 11, y = 2)
    @test (x = 1, y = 2, name = :x) |> @\((; :name => _[:name] + 10, :y)) === (
        x = 11,
        y = 2,
    )
    @test (x = 1, y = 2, name = :x) |> @\(; :name => _[:name] + 10, y = :y) === (
        x = 11,
        y = 2,
    )
    @test (x = 1, y = 2, name = :x) |> @\((; :name => _[:name] + 10, y = :y)) === (
        x = 11,
        y = 2,
    )
end

struct ObjectProxy
    dict::AbstractDict{Symbol}
end

Base.getproperty(obj::ObjectProxy, name::Symbol) = getfield(obj, :dict)[name]

@testset "NonNamedTuple" begin
    @test 1 + 2im |> @\(:re, :im) == (re = 1, im = 2)
    @test ObjectProxy(Dict(:x => 1)) |> @\(:x, y = :x + 1) == (x = 1, y = 2)
end

@testset "Shadowing" begin
    x = :nonlocal
    @test (x = 1, y = 2) |> @\[:x, :y] == [1, 2]
    @test (x = 1, y = 2) |> @\(_..., a = :x + :y, y = 10) === (x = 1, y = 10, a = 3)
    @test (x = 1, y = 2) |> @\(:x, :y) === (x = 1, y = 2)
end

@testset "Free variable" begin
    z = 3
    @test (x = 1, y = 2) |> @\[:x, :y, z] == [1, 2, 3]
    @test (x = 1, y = 2) |> @\(_..., a = :x + :y, y = z) === (x = 1, y = 3, a = 3)
    @test (x = 1, y = 2) |> @\(z, :x) === (z = 3, x = 1)
end

@testset "macrocall" begin
    z = 3
    result = (x = 1, y = 2) |> @\(:x, z = @timed(z))
    @test result.x === 1
    @test result.z[1] === z
    @test result.z[5].total_time isa Number
end

@testset "quote" begin
    @test identity |> @\(_(:(x = 1, y = 2))) == :(x = 1, y = 2)
end

@testset "Non-Symbol QuoteNode" begin
    ex = :(@\(; x=$(QuoteNode(:(some(expr))))))
    @test nothing |> @eval($ex) == (x = :(some(expr)),)
end

@testset "Explicit nonlocal" begin
    x = :nonlocal
    @test (x = 1, y = 2) |> @\(a = :x, b = x) == (a = 1, b = :nonlocal)
    @test (x = 1, y = 2) |> @\(x, :y) == (x = :nonlocal, y = 2)
    @test (x = 1, y = 2) |> @\(x = string(x, :y)) == (x = "nonlocal2",)
    @test (x = 1, y = 2) |> @\(string(x, :y)) == "nonlocal2"
    @test (x = 1, y = 2) |> @\(; x) == (x = :nonlocal,)
    @test (x = 1, y = 2) |> @\[x, :y] == [:nonlocal, 2]
end

@testset "Interpolation" begin
    x = :nonlocal
    @test (x = 1, y = 2) |> @\(a = :x, b = x, c = $:x) == (a = 1, b = :nonlocal, c = :x)
    @test (x = 1, y = 2) |> @\(identity(:x)) == 1
    @test (x = 1, y = 2) |> @\(identity(x)) == :nonlocal
    @test (x = 1, y = 2) |> @\(identity($:x)) == :x
    @test (x = 1, y = 2) |> @\($:x) == :x
    @test (x = 1, y = 2) |>
          @\(merge((; a = $:x), (; :x, y = x))) == (a = :x, x = 1, y = :nonlocal)
end

macro test_error(ex)
    quote
        let err = nothing
            @test try
                $(esc(ex))
                false
            catch err
                true
            end
            err
        end
    end
end

@testset "Error handling" begin
    ⊏ = occursin

    @testset "@\\ :arg" begin
        err = @test_error @eval @\ :arg
        @test raw"""
        Ambiguous expression: `@\(:arg)`
        Use `@\(_.:arg)` if you mean `x -> x.:arg`.
        Use `@\(; :arg)` if you mean `x -> (:arg = x.:arg,)`.
        """ ⊏ sprint(showerror, err)
    end

    @testset "@\\ arg" begin
        err = @test_error @eval @\ arg
        @test raw"""
        Ambiguous expression: `@\(arg)`
        Use `@\(_.arg)` if you mean `x -> x.arg`.
        Use `@\(; arg)` if you mean `x -> (arg = x.arg,)`.
        """ ⊏ sprint(showerror, err)
    end

    @testset "@\\(; \$(f(x)))" begin
        ex = :(@\(; $(Expr(:$, :(f(x))))))
        err = @test_error @eval $ex
        @test "syntax: invalid named tuple element \"f(x)\"" ⊏ sprint(showerror, err)
    end

    @testset "@\\(; \$(x, y))" begin
        ex = :(@\(; $(Expr(:$, :((x, y))))))
        err = @test_error @eval $ex
        @test "Only single-argument \$ is supported. Got: " ⊏ sprint(showerror, err)
    end

    @testset "@\\_" begin
        ex = :(@\_)
        err = @test_error @eval $ex
        @test "`@\\_` not supported. Use `identity`." ⊏ sprint(showerror, err)
    end

    @testset "Unsupported expression" begin
        ex = :(@\ import Base import Base)
        err = @test_error @eval $ex
        @test "Unsupported expression:" ⊏ sprint(showerror, err)
    end
end

@testset "doctest" begin
    doctest(AtBackslash; manual=false)
end

end  # module


text:
## REGISTERING LABELS OF OBJECTS DURING ASSIGNMENT

const HANDLE_GIVEN_ID = Dict{UInt64,Symbol}()

macro colon(p)
    Expr(:quote, p)
end

"""
    @constant x = value

Equivalent to `const x = value` but registers the binding thus:

    MLJBase.HANDLE_GIVEN_ID[objectid(value)] = :x

Registered objects get displayed using the variable name to which it was bound in calls to `show(x)`, etc.

WARNING: As with any `const` declaration, binding `x` to new value of
the same type is not prevented and the registration will not be updated.

"""
macro constant(ex)
    ex.head == :(=) || throw(error("Expression must be an assignment."))
    handle = ex.args[1]
    value = ex.args[2]
    quote
        const $(esc(handle)) = $(esc(value))
        id = objectid($(esc(handle)))
        HANDLE_GIVEN_ID[id] = @colon $handle
        $(esc(handle))
    end
end

"""to display abbreviated versions of integers"""
function abbreviated(n)
    as_string = string(n)
    return as_string[1]*"…"*as_string[end-1:end]
end

"""return abbreviated object id (as string)  or it's registered handle (as string) if this exists"""
function handle(X)
    id = objectid(X)
    if id in keys(HANDLE_GIVEN_ID)
        return string(HANDLE_GIVEN_ID[id])
    else
        return abbreviated(id)
    end
end

## EXPOSED SHOW METHODS

# short version of showing a `MLJType` object:
function Base.show(stream::IO, object::MLJType)
    id = objectid(object)
    description = string(typeof(object).name.name)
    parameters = typeof(object).parameters
    if length(parameters) == 1
        p_string = string(parameters[1])
        if length(p_string) > 25
            p_string = p_string[1:24]*"…"
        end
        description *= "{"*p_string*"}"
    end
    #   description = string(typeof(object))
    str = "$description @ $(handle(object))"
    if !isempty(fieldnames(typeof(object)))
        printstyled(IOContext(stream, :color=> true), str, bold=true)#color=:blue
    else
        print(stream, str)
    end
end

# longer version of showing a `MLJType` object (used at REPL):
function Base.show(stream::IO, ::MIME"text/plain", object::MLJType)
    _recursive_show(stream, object, 1, DEFAULT_SHOW_DEPTH)
end

# version showing a `MLJType` object to arbitrary depth:
function Base.show(stream::IO, object::MLJType, depth::Int)
    _recursive_show(stream, object, 1, depth)
end

# for convenience:
Base.show(object::MLJType, depth::Int) = show(stdout, object, depth)

"""
    @more

Entered at the REPL, equivalent to `show(ans, 100)`. Use to get a
recursive description of all fields of the last REPL value.

"""
macro more()
    esc(quote
        show(Main.ans, 100)
    end)
end

## METHODS TO SUPRESS THE DISPLAY OF LARGE NON-BASETYPE OBJECTS

istoobig(::Any) = true
istoobig(::DataType) = false
istoobig(::UnionAll) = false
istoobig(::Union) = false
istoobig(::Number) = false
istoobig(::Char) = false
istoobig(::Function) = false
istoobig(s::Symbol) = false
istoobig(str::AbstractString) = length(str) > 50

## THE `_show` METHOD

# Note: The `_show` method controls how field values are displayed in
# the table generated by `_recursive_show`. See top of file.

# _show fallback:
function _show(stream::IO, object)
    if !istoobig(object)
        show(stream, MIME("text/plain"), object)
        println(stream)
    else
        println(stream, "(omitted ", typeof(object), ")")
    end
end

_show(stream::IO, object::MLJType) = println(stream, object)

# _show for other types:

# function _show(stream::IO, df::AbstractDataFrame)
#     println(stream, "(omitted $(typeof(df)) of size $(size(df)))")
# end

istoobig(t::Tuple{Vararg{T}}) where T<:Union{Number,Symbol,Char,MLJType} =
    length(t) > 5
function _show(stream::IO, t::Tuple)
    if !istoobig(t)
        show(stream, MIME("text/plain"), t)
        println(stream)
    else
        println(stream, "(omitted $(typeof(t)) of length $(length(t)))")
    end
end

istoobig(A::AbstractArray{T}) where T<:Union{Number,Symbol,Char,MLJType} =
    maximum(size(A)) > 5
function _show(stream::IO, A::AbstractArray)
    if !istoobig(A)
        show(stream, MIME("text/plain"), A)
        println(stream)
    else
        println(stream, "(omitted $(typeof(A)) of size $(size(A)))")
    end
end

istoobig(d::Dict{T,Any}) where T <: Union{Number,Symbol,Char,MLJType} =
    length(keys(d)) > 5
function _show(stream::IO, d::Dict{T, Any}) where T <: Union{Number,Symbol}
    if isempty(d)
        println(stream, "empty $(typeof(d))")
    elseif !istoobig(d)
        println(stream, "omitted $(typeof(d)) with keys: ")
        show(stream, MIME("text/plain"), collect(keys(d)))
        println(stream)
    else
        println(stream, "(omitted $(typeof(d)))")
    end
end

function _show(stream::IO, v::Array{T, 1}) where T
    if !istoobig(v)
        show(stream, MIME("text/plain"), v)
        println(stream)
    else
        println(stream, "(omitted Vector{$T} of length $(length(v)))")
    end
end

_show(stream::IO, T::DataType) = println(stream, T)

_show(stream::IO, ::Nothing) = println(stream, "nothing")

## THE RECURSIVE SHOW METHOD

"""
    _recursive_show(stream, object, current_depth, depth)

Generate a table of the field values of the `MLJType` object,
dislaying each value by calling the method `_show` on it. The
behaviour of `_show(stream, f)` is as follows:

1. If `f` is itself a `MLJType` object, then its short form is shown
and `_recursive_show` generates as separate table for each of its
field values (and so on, up to a depth of argument `depth`).

2. Otherwise `f` is displayed as "(omitted T)" where `T = typeof(f)`,
unless `istoobig(f)` is false (the `istoobig` fall-back for arbitrary
types being `true`). In the latter case, the long (ie,
MIME"plain/text") form of `f` is shown. To override this behaviour,
overload the `_show` method for the type in question.

"""
function _recursive_show(stream::IO, object::MLJType, current_depth, depth)
    if depth == 0 || isempty(fieldnames(typeof(object)))
        println(stream, object)
    elseif current_depth <= depth
        fields = fieldnames(typeof(object))
        print(stream, "#"^current_depth, " ")
        show(stream, object)
        println(stream, ": ")
#        println(stream)
        if isempty(fields)
            println(stream)
            return
        end
        for fld in fields
            fld_string = string(fld)*" "^(max(0,COLUMN_WIDTH - length(string(fld))))*"=>   "
            print(stream, fld_string)
            if isdefined(object, fld)
                _show(stream, getfield(object, fld))
                #           println(stream)
            else
                println(stream, "(undefined)")
                #           println(stream)
            end
        end
        println(stream)
        for fld in fields
            if isdefined(object, fld)
                subobject = getfield(object, fld)
                if isa(subobject, MLJType) && !isempty(fieldnames(typeof(subobject)))
                    _recursive_show(stream, getfield(object, fld),
                                    current_depth + 1, depth)
                end
            end
        end
    end
end


text:
######
##FEM Heat Nonlinear Test
######
using FiniteElementDiffEq, DiffEqProblemLibrary, DiffEqPDEBase

#Define a parabolic problem
prob = prob_femheat_birthdeath

#Solve it with a bunch of different algorithms, plot solution
println("Euler")
sol = solve(prob,FEMDiffEqHeatEuler())

println("Semi-implicit Euler")
sol = solve(prob,FEMDiffEqHeatSemiImplicitEuler())

println("Semi-implicit Crank Nicholson")
sol = solve(prob,FEMDiffEqHeatSemiImplicitCrankNicholson())

dx = 1//2^(2)
dt = 1//2^(3)
T = 1
fem_mesh = parabolic_squaremesh([0 1 0 1],dx,dt,T,:neumann)
println("Implicit Euler")
sol = solve(prob,FEMDiffEqHeatImplicitEuler(),autodiff=true)
#println("Implicit Euler Autodiff False")
#sol = solve(prob,FEMDiffEqHeatImplicitEuler(),autodiff=false)
#plot(sol)

#Returns true if nonlinear solver is correct
bool1 = maximum(abs.(sol[end] .- .777))<.01

### Stochastic Tests

#Define a parabolic problem
T = 1
dx = 1//2^(3)
dt = 1//2^(7)
fem_mesh = parabolic_squaremesh([0 1 0 1],dx,dt,T,:neumann)
prob = prob_femheat_stochasticbirthdeath

#Solve it with a bunch of different algorithms, plot solution
println("Euler")
sol = solve(prob,FEMDiffEqHeatEuler())

println("Semi-implicit Euler")
sol = solve(prob,FEMDiffEqHeatSemiImplicitEuler())

#=
# CG and GMRES require size 1 vector, breaks with numvars change
println("Semi-implicit Crank Nicholson")
sol = solve(prob,alg=:SemiImplicitCrankNicholson,solver=:CG)

println("Semi-implicit Crank Nicholson GMRES")
sol = solve(prob,alg=:SemiImplicitCrankNicholson,solver=:GMRES)
=#

#Define a quicker problem
prob =  prob_femheat_stochasticbirthdeath_fast
println("Implicit Euler")
sol = solve(prob,FEMDiffEqHeatImplicitEuler(),autodiff=true)
sol = solve(prob,FEMDiffEqHeatImplicitEuler(),autodiff=false)
#plot(sol)

bool1


text:
@testset "Deprecations" begin

    @testset "NCC (#63)" begin
        type_list = generate_test_types([Bool, Float32, N0f8], [Gray])
        A = [0.0 0.0 0.0; 1.0 1.0 1.0; 0.0 0.0 0.0]
        B = [1.0 1.0 1.0; 0.0 0.0 0.0; 1.0 1.0 1.0]
        for T in type_list
            a, b = T.(A), T.(B)
            @test nearlysame(ncc(a, b), zncc(a, b))
            @test nearlysame(NCC()(a, b), ZNCC()(a, b))
        end

        type_list = generate_test_types([Float32, N0f8], [RGB])
        A = [RGB(0.0, 0.0, 0.0) RGB(0.0, 1.0, 0.0) RGB(0.0, 1.0, 1.0)
            RGB(0.0, 0.0, 1.0) RGB(1.0, 0.0, 0.0) RGB(1.0, 1.0, 0.0)
            RGB(1.0, 1.0, 1.0) RGB(1.0, 0.0, 1.0) RGB(0.0, 0.0, 0.0)]
        B = [RGB(0.0, 0.0, 0.0) RGB(0.0, 0.0, 1.0) RGB(1.0, 1.0, 1.0)
            RGB(0.0, 1.0, 0.0) RGB(1.0, 0.0, 0.0) RGB(1.0, 0.0, 1.0)
            RGB(0.0, 1.0, 1.0) RGB(1.0, 1.0, 0.0) RGB(0.0, 0.0, 0.0)]
        for T in type_list
            a, b = T.(A), T.(B)
            @test nearlysame(ncc(a, b), zncc(a, b))
            @test nearlysame(NCC()(a, b), ZNCC()(a, b))
        end
    end
end


text:
using DiffPrecTest
using Statistics, StatsBase, LinearAlgebra
using SparseArrays
using ProximalBase, CoordinateDescent, CovSel
using Random, Distributions
using JLD

rep   = parse(Int,ARGS[1])
dir = ARGS[2]

p = 100
n = 300

# generate model
Ωx = Matrix{Float64}(I, p, p)
mf = [1., 0.5, 0.4]
# block 1
α = 1.
bp = 1
ep = 50
for k=0:2
    v = α * mf[k+1]
    for l=bp:ep-k
        Ωx[l  , l+k] = v
        Ωx[l+k, l  ] = v
    end
end
# block 2
α = 2.
bp = 51
ep = 75
for k=0:2
    v = α * mf[k+1]
    for l=bp:ep-k
        Ωx[l  , l+k] = v
        Ωx[l+k, l  ] = v
    end
end
# block 3
α = 4.
bp = 76
ep = 100
for k=0:2
    v = α * mf[k+1]
    for l=bp:ep-k
        Ωx[l  , l+k] = v
        Ωx[l+k, l  ] = v
    end
end
Ωy = Matrix(Diagonal(Ωx))
Δ = Ωx - Ωy

Σx = inv(Symmetric(Ωx))
Σy = inv(Symmetric(Ωy))

dist_X = MvNormal(convert(Matrix, Σx))
dist_Y = MvNormal(convert(Matrix, Σy))

# generate data
Random.seed!(2345 + rep)
X = rand(dist_X, n)'
Y = rand(dist_Y, n)'
Xtest = rand(dist_X, n)'    # validation data for
Ytest = rand(dist_Y, n)'

#################
#
#   our method
#
#################

τArr = collect( range(10,0,length=500) )

Sx = Symmetric( cov(X) )
Sy = Symmetric( cov(Y) )
Sxtest = Symmetric( cov(Xtest) )
Sytest = Symmetric( cov(Ytest) )

@time eS = DiffPrecTest.__initSupport(Sx, Sy, X, Y)
@time eΔNormal, eΔNormalP, _ = supportEstimate(ANTSupport(), X, Y, τArr; estimSupport=eS)
@time eΔBoot  , eΔBootP, _ = supportEstimate(BootStdSupport(), X, Y, τArr; estimSupport=eS)

# compete method

maxλ = maximum( abs.(Sx - Sy) )
minλ = maxλ * 0.04
Λarr = exp.(range(log(maxλ), log(minλ), length=50))

eΔDTr, i2, iInf, loss2arr, lossInfarr = supportEstimate(DTraceValidationSupport(), Sx, Sy, Sxtest, Sytest, Λarr)

@save "$(dir)/res_$(rep).jld" eΔNormal eΔNormalP eΔBoot eΔBootP eS eΔDTr i2 iInf loss2arr lossInfarr


text:
@with_kw mutable struct TermCache{T <: AbstractFloat}
    best_resid::T = Inf
    non_improv_iters::Int64 = 0
end

function unpack_sol(iter::Iterate{T}, prob::Problem{T}, stats::Stats{T}, settings::Settings{T}) where T
    sol = Solution(prob)
    prim_obj = prob.c'*iter.x
    dual_obj = prob.b'*iter.y + prob.h'*iter.z
    if iter.τ ≤ EPS
        sol.x = iter.x
        sol.y = iter.y
        sol.z = iter.z
        sol.s = iter.s
        if prim_obj < EPS
            sol.status = unbounded_status
            sol.optval = -Inf
        elseif dual_obj < EPS
            sol.status = infeasible_status
            sol.optval = Inf
        else
            sol.status = numerical_error_status
            sol.optval = NaN
        end
    else
        sol.status = optimal_status
        sol.x = iter.x / iter.τ
        sol.y = iter.y / iter.τ
        sol.z = iter.z / iter.τ
        sol.s = iter.s / iter.τ
        sol.optval = prim_obj / iter.τ
    end
    if stats.iters == settings.max_iters
        sol.status = max_iters_status
    end
    return sol
end

function check_term_conds!(iter::Iterate{T}, best_iter::Iterate{T}, term_cache::TermCache{T},
                           resid::Residual{T}, prob::Problem{T}, settings::Settings{T},
                           stats::Stats{T}) where T
    norm_resid = norm(stack(resid), Inf)
    terminate = norm_resid ≤ settings.tol

    if norm_resid ≥ settings.non_improv_ratio_thres * term_cache.best_resid
        term_cache.non_improv_iters += 1
    else
        term_cache.non_improv_iters = 0
    end
    if term_cache.non_improv_iters > settings.max_non_improv_iters
        terminate = true
    end

    if norm_resid < term_cache.best_resid
        term_cache.best_resid = norm_resid
        copy_to!(iter, best_iter)
        stats.best[stats.iters] = true
    else
        stats.best[stats.iters] = false
    end
    stats.residual[stats.iters] = norm_resid

    return terminate
end


text:
using Documenter, GtkReactive

makedocs(sitename = "GtkReactive",
         format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
         pages    = ["index.md", "controls.md", "drawing.md", "zoom_pan.md", "reference.md"]
         )

deploydocs(repo   = "github.com/JuliaGizmos/GtkReactive.jl.git",
           )


text:
"""
    bars(values::Array;
            yheight = 200,
            xwidth = 25,
            labels = true,
            barfunction = f,
            labelfunction = f,
        )

Draw some bars where each bar is the height of a value in the array. The bars
will fit in a box `yheight` high (even if there are negative values).

To control the drawing of the text and bars, define functions that process the
end points:

`mybarfunction(bottom::Point, top::Point, value; extremes=[a, b], barnumber=0,
bartotal=0)`

`mylabelfunction(bottom::Point, top::Point, value; extremes=[a, b], barnumber=0,
bartotal=0)`

and pass them like this:

```julia
bars(v, yheight=10, xwidth=10, barfunction=mybarfunction)
bars(v, xwidth=15, yheight=10, labelfunction=mylabelfunction)
```

or:

```
bars(v, labelfunction = (args...; extremes=[], barnumber=0, bartotal=0) ->  setgray(rand()))
```

To suppress the text labels, use optional keyword `labels=false`.
"""
function bars(values::Array;
    yheight = 200,
    xwidth = 25,
    barfunction   = (bottom::Point, top::Point, value;
        extremes=extrema(values), barnumber=0, bartotal=0) -> begin
                setline(xwidth)
                line(bottom, top, :stroke)
            end,
    labels::Bool=true,
    labelfunction = (bottom::Point, top::Point, value;
        extremes=extrema(values), barnumber=0, bartotal=0) -> begin
            t = string(round(value, digits=2))
            textoffset = textextents(t)[4]
            fontsize(10)
            if top.y < 0
                tp = Point(top.x, min(top.y, bottom.y) - textoffset)
            else
                tp = Point(top.x, max(top.y, bottom.y) + textoffset)
            end
            text(t, tp, halign=:center, valign=:middle)
        end)

    x = O.x
    mn, mx = extrema(values)
    isapprox(mn, mx, atol=0.00001) && (mx = mn + 100) # better show something than nothing
    for (n, v) in enumerate(values)
        # remember y increases downwards by default
        bottom = Point(x, -rescale(min(v, 0), min(0, mn), mx, 0, yheight))
        top    = Point(x, -rescale(max(v, 0), min(0, mn), mx, 0, yheight))
        barfunction(bottom, top, v, extremes=extrema(values), barnumber=n, bartotal=length(values))
        labels && labelfunction(bottom, top, v, extremes=extrema(values), barnumber=n, bartotal=length(values))
        x += xwidth
    end
end


text:
using Plots
using StatsBase
import StatsBase.std

function pysafelabel(s)
    if s isa String
        replace(s, r"_" => s"\\_")
    else
        return s
    end
end
export pysafelabel

include("Plotting/ParameterEstimate.jl")
include("Plotting/DimensionalityEstimate.jl")
include("Plotting/TimeSeriesPlots.jl")
include("Plotting/FeatureDistributions.jl")
include("Plotting/LowDimensionalPlots.jl")
include("Plotting/StatisticalFilter.jl")
include("Plotting/ToyVarianceMaps.jl")
include("Plotting/FeatureImage.jl")
include("Plotting/transparentdark.jl")
include("Plotting/FeatureResponses.jl")
#include("Plotting/WordCloud.jl")


text:
# Chemfiles.jl, a modern library for chemistry file reading and writing
# Copyright (C) Guillaume Fraux and contributors -- BSD license

__ptr(property::Property) = __ptr(property.__handle)
__const_ptr(property::Property) = __const_ptr(property.__handle)

@enum PropertyKind begin
    PropertyBool = lib.CHFL_PROPERTY_BOOL
    PropertyDouble = lib.CHFL_PROPERTY_DOUBLE
    PropertyString = lib.CHFL_PROPERTY_STRING
    PropertyVector3d = lib.CHFL_PROPERTY_VECTOR3D
end

PropertyValue = Union{Bool,Float64,String,Vector{Float64}}

# Create a `Bool` `Property`.
function Property(value::Bool)
    ptr = @__check_ptr(lib.chfl_property_bool(convert(UInt8, value)))
    return Property(CxxPointer(ptr, is_const=false))
end

# Create a `Float64` `Property`.
function Property(value::Float64)
    ptr = @__check_ptr(lib.chfl_property_double(value))
    return Property(CxxPointer(ptr, is_const=false))
end

# Create a `String` `Property`.
function Property(value::String)
    ptr = @__check_ptr(lib.chfl_property_string(pointer(value)))
    return Property(CxxPointer(ptr, is_const=false))
end

# Create a `Vector3D` `Property`.
function Property(value::Vector{Float64})
    @assert size(value) == (3,)
    ptr = @__check_ptr(lib.chfl_property_vector3d(value))
    return Property(CxxPointer(ptr, is_const=false))
end

# Get the kind of a property.
function kind(property::Property)
    result = Ref{lib.chfl_property_kind}(0)
    __check(lib.chfl_property_get_kind(__const_ptr(property), result))
    return PropertyKind(result[])
end

function __extract_double(property::Property)
    result = Ref{Cdouble}(0)
    __check(lib.chfl_property_get_double(__const_ptr(property), result))
    return result[]
end

function __extract_bool(property::Property)
    result = Ref{UInt8}(0)
    __check(lib.chfl_property_get_bool(__const_ptr(property), result))
    return convert(Bool, result[])
end

function __extract_string(property::Property)
    return __call_with_growing_buffer(
        (buffer, size) -> __check(lib.chfl_property_get_string(__const_ptr(property), buffer, size))
    )
end

function __extract_vector3d(property::Property)
    result = Float64[0, 0, 0]
    __check(lib.chfl_property_get_vector3d(__const_ptr(property), result))
    return result
end

# Obtain the value stored by a property.
function extract(property::Property)::PropertyValue
    property_kind = kind(property)
    if property_kind == PropertyBool
        return __extract_bool(property)
    elseif property_kind == PropertyDouble
        return __extract_double(property)
    elseif property_kind == PropertyString
        return __extract_string(property)
    elseif property_kind == PropertyVector3d
        return __extract_vector3d(property)
    else
        throw(ChemfilesError(
            "Invalid property kind '$property_kind'. This is a bug"
        ))
    end
end


text:
module AlgebraicRelations
  include("DB.jl")
  include("Queries.jl")
  include("Interface.jl")
  include("Presentations.jl")
end


text:
using Hashpipe
using Test

@testset "Hashpipe.jl" begin
    @testset "Status" begin
        instance_id = 98
        @test status_exists(instance_id) == 0
        status = status_t(0,0,0,0)
        @test status_attach(instance_id, Ref(status)) == 0

        @test status_exists(instance_id) == 1
        @test status_unlock(Ref(status)) == 0
        update_status(status, "JLTEST", "RUNNING")
        status_buf_lock_unlock(Ref(status)) do
            update_status(status, "JLTEST", "COMPLETE")
        end
        @test status_lock(Ref(status)) == 0

        # Detach from and delete shared memory segment
        status_detach(Ref(status))
        ipcrm_cmd = `ipcrm -m $(status.shmid)`
        run(ipcrm_cmd)
        @test status_exists(instance_id) == 0
    end

    @testset "Databuf" begin
        # TODO: Test creation, data retrieval, updating, locking, unlocking
        # TODO: Test hput/get functions
        #db = Hashpipe.databuf_create()
    end

    @testset "Thread" begin
        # Test registration and finding
    end
end


text:
# Utils for filter performance
function fill_a1!(a::Matrix{T}) where T
    @inbounds for i in axes(a, 2)
        a[1, i] = zero(T)
    end
end
function fill_P1!(P::AbstractArray{T}; bigkappa::T = 1e6) where T
    @inbounds for i in axes(P, 1), j in axes(P, 2)
        if i == j
            P[i, j, 1] = bigkappa
        else
            P[i, j, 1] = zero(T)
        end
    end
end

function repeat_matrix_t_plus_1!(mat::AbstractArray{T}, t::Int) where T
    @inbounds for j in axes(mat, 2), i in axes(mat, 1)
          mat[i, j, t+1] = mat[i, j, t]
    end
    return
end
function repeat_vector_t_plus_1!(mat::Matrix{T}, t::Int) where T
    @inbounds for i in axes(mat, 1)
        mat[i, t+1] = mat[i, t]
    end
    return
end

function big_update_a!(a::Matrix{Typ}, att::Matrix{Typ}, T::Matrix{Typ}, c::Matrix{Typ}, t::Int) where Typ
    @views @inbounds mul!(a[t+1, :], T, att[t, :])
    for i in axes(a, 2)
        a[t + 1, i] += c[t, i]
    end
    return
end
function small_update_a!(a::Matrix{Typ}, att::Matrix{Typ}, T::Matrix{Typ}, c::Matrix{Typ}, t::Int) where Typ
    @inbounds for i in axes(a, 2)
        a[t+1, i] = c[t, i]
        for j in axes(T, 2)
            a[t+1, i] += T[i, j]*att[t, j]
        end
    end
    return
end
function update_a!(a::Matrix{Typ}, att::Matrix{Typ}, T::Matrix{Typ}, c::Matrix{Typ}, t::Int) where Typ
    # Here there is a trade-off between memory and speed, usually if the dimension of a is smaller than 15
    # it is more performant to do the hard-coded version (small_update_a)
    if size(a, 2) < 16
        small_update_a!(a, att, T, c, t)
    else
        big_update_a!(a, att, T, c, t)
    end
    return
end

function update_ZP!(ZP::Matrix{T}, Z::Array{T, 3}, P::Array{T, 3}, t::Int) where T
    @inbounds @views mul!(ZP, Z[:, :, t], P[:, :, t])
    return
end
function update_ZP!(ZP::Vector{T}, Z::Array{T, 3}, P::Array{T, 3}, t::Int) where T
    @inbounds for i in axes(ZP, 1)
        ZP[i] = zero(T)
        for j in axes(ZP, 1)
            ZP[i] += Z[1, j, t]*P[j, i, t]
        end
    end
    return
end

function update_K!(K::AbstractArray{Typ}, P_Ztransp_invF::AbstractMatrix{Typ}, T::AbstractMatrix{Typ}, t::Int) where Typ
    @inbounds @views mul!(K[:, :, t], T, P_Ztransp_invF)
    return
end

function update_K!(K::AbstractMatrix{Typ}, P_Ztransp_invF::AbstractVector{Typ}, T::AbstractMatrix{Typ}, t::Int) where Typ
    # Here there is a trade-off memory-speed, usually if the dimension of a is smaller than 15
    # it is more performant to do the hard coded version (small_update_K!)
    if size(K, 1) < 16
        small_update_K!(K, P_Ztransp_invF, T, t)
    else
        big_update_K!(K, P_Ztransp_invF, T, t)
    end
    return
end

function small_update_K!(K::AbstractMatrix{Typ}, P_Ztransp_invF::AbstractVector{Typ}, T::AbstractMatrix{Typ}, t::Int) where Typ
    @inbounds for i in axes(K, 1)
        K[i, t] = zero(Typ)
        for j in axes(K, 1)
            K[i, t] += T[i, j]*P_Ztransp_invF[j]
        end
    end
    return
end

function big_update_K!(K::AbstractMatrix{Typ}, P_Ztransp_invF::AbstractVector{Typ}, T::AbstractMatrix{Typ}, t::Int) where Typ
    @inbounds @views K[:, t] = T * P_Ztransp_invF
    return
end

function invF(F::AbstractArray{T}, t::Int) where T
    return @inbounds @views invertF(F, t)
end

function update_P!(P::Array{Typ, 3}, T::Matrix{Typ}, Ptt::Array{Typ, 3}, RQR::Matrix{Typ}, t::Int)  where Typ
    @views @inbounds LinearAlgebra.BLAS.gemm!('N', 'T', Typ(1.0), T * Ptt[:, :, t], T, Typ(0.0), P[:, :, t+1]) # P[:, :, t+1] = T * Ptt[:, :, t] * T'
    sum_matrix!(P, RQR, t, 1) # P[:, :, t+1] .+= RQR
    ensure_pos_sym!(P, t + 1)
    return
end

function update_v!(v::Matrix{T}, y::Matrix{T}, Z::Array{T, 3}, d::Matrix{T}, a::Matrix{T}, t::Int) where T
    # v[t, :] = y[t, :] - Z[:, :, t]*a[t, :]
    @inbounds for i in axes(Z, 1)
        v[t, i] = y[t, i] - d[t, i]
        for j in axes(Z, 2)
            v[t, i] -= Z[i, j, t]*a[t, j]
        end
    end
    return
end

function update_v!(v::Vector{T}, y::Matrix{T}, Z::Array{T, 3}, d::Matrix{T}, a::Matrix{T}, t::Int) where T
    # v[t, :] = y[t, :] - Z[:, :, t]*a[t, :]
    v[t] = y[t, 1] - d[t, 1]
    @inbounds for j in axes(Z, 2)
        v[t] -= Z[1, j, t]*a[t, j]
    end
    return
end

function update_att!(att::Matrix{T}, a::Matrix{T}, P_Ztransp_invF::Matrix{T}, v::Matrix{T}, t::Int) where T
    # att[t, :] = a[t, :] + P_Ztransp_invF*v[t, :]
    @inbounds for i in axes(P_Ztransp_invF, 1)
        att[t, i] = a[t, i]
        for j in axes(P_Ztransp_invF, 2)
            att[t, i] += P_Ztransp_invF[i, j]*v[t, j]
        end
    end
    return
end

function update_att!(att::Matrix{T}, a::Matrix{T}, P_Ztransp_invF::Vector{T}, v::Vector{T}, t::Int) where T
    # att[t, :] = a[t, :] + P_Ztransp_invF*v[t, :]
    @inbounds for i in axes(P_Ztransp_invF, 1)
        att[t, i] = a[t, i] + P_Ztransp_invF[i]*v[t]
    end
    return
end

function update_att!(att::Matrix{T}, a::Matrix{T}, t::Int) where T
    @inbounds for j in axes(att, 2)
        att[t, j] = a[t, j]
    end
    return
end

function update_Ptt!(Ptt::Array{T, 3}, P::Array{T, 3}, t::Int) where T
    @inbounds for i in axes(Ptt, 1), j in axes(Ptt, 2)
        Ptt[i, j, t] = P[i, j, t]
    end
    return
end

function update_Ptt!(Ptt::Array{T, 3}, P::Array{T, 3}, P_Ztransp_invF::Matrix{T},
                    ZP::Matrix{T}, t::Int) where T
    @views @inbounds mul!(Ptt[:, :, t], -P_Ztransp_invF, ZP) # Ptt_t = - P_Ztransp_invF * ZP
    sum_matrix!(Ptt, P, t, 0) # Ptt_t += P_t
    return
end
function update_Ptt!(Ptt::Array{T, 3}, P::Array{T, 3}, P_Ztransp_invF::Vector{T},
                    ZP::Vector{T}, t::Int) where T
    @inbounds for i in axes(Ptt, 1), j in axes(Ptt, 2)
        Ptt[i, j, t] = P[i, j, t] - P_Ztransp_invF[i]*ZP[j] # Ptt_t = - P_Ztransp_invF * ZP
    end
    return
end

function update_F!(F::Array{T, 3}, ZP::Matrix{T}, Z::Array{T, 3}, H::AbstractArray{T}, t::Int) where T
    @views @inbounds LinearAlgebra.BLAS.gemm!('N', 'T', T(1.0), ZP, Z[:, :, t], T(0.0), F[:, :, t]) # F_t = Z_t * P_t * Z_t
    sum_matrix!(F, H, t, 0) # F[:, :, t] .+= H
    return
end

function update_F!(F::Vector{T}, ZP::Vector{T}, Z::Array{T, 3}, H::Matrix{T}, t::Int) where T
    F[t] = H[1, 1]
    @inbounds for i in axes(ZP, 1)
        F[t] += ZP[i]*Z[1, i, t] # F_t = Z_t * P_t * Z_t
    end
    return
end

function update_P_Ztransp_Finv!(P_Ztransp_invF::Matrix{T}, ZP::Matrix{T}, F::Array{T, 3}, t::Int) where T
    # P_Ztransp_invF = (ZP)' * F^-1
    if size(F, 1) == 1
        for p in axes(ZP, 1), m in axes(ZP, 2)
            P_Ztransp_invF[m, p] = ZP[p, m]/F[1, 1, t]
        end
    else
        LinearAlgebra.BLAS.gemm!('T', 'N', T(1.0), ZP, invF(F, t), T(0.0), P_Ztransp_invF)
    end
    return
end

function update_P_Ztransp_Finv!(P_Ztransp_invF::Vector{T}, ZP::Vector{T}, F::Vector{T}, t::Int) where T
    @inbounds for i in axes(P_Ztransp_invF, 1)
        P_Ztransp_invF[i] = ZP[i]/F[t]
    end
    return
end