text:
function solveDiffusion_threaded(
    color_list,
    CMidx,
    KT,
    MT,
    FT,
    DoF_T,
    coordinates,
    VarT::ThermalParameters,
    ScratchDifussion::ScratchThermalGlobal,
    SF_Diffusion::ShapeFunctionsThermal,
    valA,
    ρ,
    Δt,
    T,
    T0,
    TBC,
    to,
)

    ∂Ωt = TBC.Ω
    tfix = TBC.vfix
    tfree = TBC.ifree
    copyto!(T0, T)

    @timeit to "Thermal diffusion threaded" begin
        # Reset Matrices
        fill!(KT.nzval, 0.0)
        fill!(MT.nzval, 0.0)
        fill!(FT, 0.0)

        # Element assembly
        assembly_threaded!(
            color_list,
            CMidx,
            KT,
            MT,
            FT,
            DoF_T,
            coordinates,
            VarT,
            ScratchDifussion,
            SF_Diffusion,
            valA,
            ρ,
            Δt,
        )

        # Apply Boundary conditions
        _prepare_matrices!(MT, FT, KT, T, Δt)
        _apply_bcs!(T, MT, FT, ∂Ωt, tfix)

        # Solve temperature
        # T[tfree] .= MT[tfree, tfree] \ FT[tfree]
        T[tfree], = Krylov.gmres(MT[tfree, tfree], FT[tfree])

        # Temperature increment
        ΔT = @tturbo T .- T0

    end
    return T, T0, ΔT, to
end

"""
assembly_threaded!(color_list, K, M, F, dof, coordinates::ElementCoordinates,
                           VarT::ThermalParameters, A::ScratchThermalGlobal,
                           B::ShapeFunctionsThermal, valA, ρ, Δt)

Threaded assembly of FEM matrices for thermal diffusion. Based on a colored mesh.
"""
function assembly_threaded!(
    color_list,
    CMidx,
    KT,
    MT,
    FT,
    DoF_T,
    coordinates::ElementCoordinates,
    VarT::ThermalParameters,
    ScratchDifussion::ScratchThermalGlobal,
    SF_Diffusion::ShapeFunctionsThermal,
    valA,
    ρ,
    Δt,
)

    ## Shape function and their derivatives
    ni, nn, nn3 = Val(ScratchDifussion.nip), Val(ScratchDifussion.nnodel), Val(3)
    N, ∇N, dN3ds, w_ip = _get_SF(ni, nn)
    N3, = _get_SF(ni, nn3)

    for color in color_list
        Threads.@threads for element in color
            assemble_element!(
                element,
                CMidx,
                KT,
                MT,
                FT,
                DoF_T,
                coordinates.θ,
                coordinates.r,
                VarT,
                ScratchDifussion,
                N,
                ∇N,
                dN3ds,
                w_ip,
                N3,
                valA,
                ρ,
                Δt,
            )
        end
    end
    return KT, MT, FT
end

function assemble_element!(
    iel,
    CMidx,
    KT,
    MT,
    FT,
    DoF_T,
    θ,
    r,
    VarT::ThermalParameters,
    A::ScratchThermalGlobal,
    N,
    ∇N,
    dN3ds,
    w_ip,
    N3,
    valA,
    ρ,
    Δt,
)

    # Unpack physical parameters
    κ, Cp, dQdT = VarT.κ, VarT.Cp, VarT.dQdT
    # Polar coordinates of element nodes
    θ_el = @SVector [θ[i, iel] for i in 1:3]
    r_el = @SVector [r[i, iel] for i in 1:3]
    coords = SMatrix{3,2}([θ_el r_el])

    # Jacobian n. 1 (p:=polar, l:=local): reference element --> current element
    J_pl = dN3ds * coords
    # detJ_pl = det(J_pl) # fast |-> unrolled
    # detJ_pl = J_pl[1] * J_pl[4] - J_pl[2] * J_pl[3] # Determinant of Jacobi matrix
    detJ_pl = muladd(J_pl[1], J_pl[4], -J_pl[2]*J_pl[3]) # Determinant of Jacobi matrix

    R_21 = r_el[2] - r_el[1] # = -detJa_PL*deta_dth
    R_31 = r_el[3] - r_el[1] # =  detJa_PL*dxi_dth
    Th_31 = θ_el[3] - θ_el[1]   # = -detJa_PL*dxi_dr
    Th_21 = θ_el[2] - θ_el[1]   # =  detJa_PL*deta_dr

    # the Jacobian ∂ξ∂θ to transform local (ξ, η) into global (θ,r) derivatives
    #     ∂ξ∂θ = [ R_31    -R_21
    #             -Th_31    Th_21] / detJa_PL
    # ∂ξ∂θ = @SMatrix [(r_el[3]-r_el[1]) (r_el[2]-r_el[1])
    #                  (θ_el[3]-θ_el[1])  (θ_el[2]-θ_el[1])]

    Ke = empty_element_matrices(valA)
    Me = empty_element_matrices(valA)
    fe = empty_element_force_vector(valA)

    el_dofs = DoF_T[iel]
    ρ_el = @SVector [ρ[el_dofs[i]] for i in 1:3]
    dQdT_el = @SVector [dQdT for _ in 1:6]

    # INTEGRATION LOOP
    @inbounds for ip in 1:(A.nip)
        # Unpack shape functions
        N_ip = N[ip]
        N3_ip = N3[ip]
        ∇N_ip = ∇N[ip]

        # ρ at ith integration point
        ρ_ip = mydot(ρ_el, N3_ip)

        # Polar coordinates of the integration points
        θ_ip = mydot(θ_el, N3_ip)
        r_ip = mydot(r_el, N3_ip)
        cos_ip = cos(θ_ip)
        sin_ip = sin(θ_ip)
        cos_ip_r_ip = cos_ip / r_ip
        sin_ip_r_ip = sin_ip / r_ip
        # Build inverse of the 2nd Jacobian
        # invJ_double = @SMatrix [(∂ξ∂θ[1,1]*cos_ip_r_ip-∂ξ∂θ[2,1]*sin_ip) (-∂ξ∂θ[1,2]*cos_ip_r_ip+∂ξ∂θ[2,2]*sin_ip);
        #                        -(∂ξ∂θ[1,1]*sin_ip_r_ip+∂ξ∂θ[2,1]*cos_ip) (∂ξ∂θ[1,2]*sin_ip_r_ip+∂ξ∂θ[2,2]*cos_ip) ]
        invJ_double = @SMatrix [
            R_31 * cos_ip_r_ip-Th_31 * sin_ip -R_21 * cos_ip_r_ip+Th_21 * sin_ip
            -R_31 * sin_ip_r_ip-Th_31 * cos_ip R_21 * sin_ip_r_ip+Th_21 * cos_ip
        ]

        # Partial derivatives
        ∂N∂x = invJ_double * ∇N_ip / detJ_pl
        # Integration weight
        ω = r_ip * detJ_pl * w_ip[ip]

        # Update elemental matrices
        NxN = N_ip' * N_ip
        Ke += (∂N∂x' * κ * ∂N∂x) * ω
        Me += NxN * ω * ρ_ip * Cp

        # Force vector -- right hand side
        if dQdT != 0
            fe += vec(N_ip) .* dQdT_el * (Δt*ω)
        end
    end

    # Update stiffness matrix
    return update_stiffness_matrix!(KT, MT, el_dofs, Ke, Me)

    # Update force vector
    if dQdT != 0
        update_force_vector!(FT, fe, el_dofs)
    end
end

empty_element_matrices(::Val{sz}) where {sz} = zeros(SMatrix{sz,sz})

empty_element_matrices(::Val{sI}, ::Val{sJ}) where {sI, sJ} = zeros(SMatrix{sI, sJ})

empty_element_force_vector(::Val{sz}) where {sz} = zeros(SVector{sz})

function update_stiffness_matrix!(K::SparseMatrixCSC{Float64,Int64}, el_dofs, Ke)
    @inbounds for (_j, j) in enumerate(el_dofs), (_i, i) in enumerate(el_dofs)
        K[i, j] += Ke[_i, _j]
    end
end

function update_stiffness_matrix!(
    K::SparseMatrixCSC{Float64,Int64}, M::SparseMatrixCSC{Float64,Int64}, el_dofs, Ke, Me
)
    for (_j, j) in enumerate(el_dofs), (_i, i) in enumerate(el_dofs)
        @inbounds K[i, j] += Ke[_i, _j]
        @inbounds M[i, j] += Me[_i, _j]
    end
end

function update_force_vector!(F::Vector{Float64}, fe, el_dofs)
    @inbounds for (_i, i) in enumerate(el_dofs)
        F[i] += fe[_i]
    end
end

get_NxN(N) = N' * N

function matrices_diffusion(els, nnodel, nnod)
    _, K = sparsitythermal(els, nnodel) # stiffness matrix
    M = deepcopy(K) # mass matrix
    F = fill(0.0, nnod) # force vector

    return K, M, F
end

function diffusion_immutables(gr, ndim, nvert, nnodel, nip)

    EL2NOD, nnod, nel = gr.e2n, gr.nnod, gr.nel

    ## Constant mesh values
    ScratchDifussion = ScratchThermalGlobal(
        ndim, # := model dimension
        nvert, # := number of vertices per element
        nnodel, # := number of nodes per element
        nel, # := number of elements
        nip,
    ) # := number of vertices per element

    # Shape function and their derivatives
    ni, nn, nn3 = Val(ScratchDifussion.nip), Val(ScratchDifussion.nnodel), Val(3)
    N, ∇N, dN3ds, w_ip = _get_SF(ni, nn)
    N3, = _get_SF(ni, nn3)

    SF_Diffusion = ShapeFunctionsThermal{Float64}(N, ∇N, dN3ds, w_ip, N3)

    # Degrees of freedom per element
    DoF_T = DoFs_Thermal(EL2NOD, ScratchDifussion.nnodel)

    # Diffusion sparse matrices and force vector
    KT, MT, FT = matrices_diffusion(EL2NOD, ScratchDifussion.nnodel, nnod)

    return KT, MT, FT, DoF_T, nn, SF_Diffusion, ScratchDifussion
end

max_stars_repo_name: albert-de-montserrat/Persephone

max_issues_repo_name: albert-de-montserrat/Persephone

max_forks_repo_name: albert-de-montserrat/Persephone


text:
 #returns the specification symbol in multicomponent
#returns :singlecomponent if there are no specifications
#throws error if more than one specification is found

const AMOUNT_CONST =Dict{Int,Any}(
     0000 => OneMol()
    ,0001 => MaterialAmount{MOLAR}()
    ,0002 => MaterialAmount{MASS}()
    ,0099 => SingleComponent()
    ,1100 => MaterialCompounds{MOLAR,FRACTION}()
    ,1200 => MaterialCompounds{MASS,FRACTION}()
    ,2100 => MaterialCompounds{MOLAR,TOTAL_AMOUNT}()
    ,2200 => MaterialCompounds{MASS,TOTAL_AMOUNT}()
    ,3100 => HumiditySpec{WetBulbTemperature}()
    ,3200 => HumiditySpec{HumidityRatio}()
    ,3300 => HumiditySpec{MolarHumidity}()
    ,3400 => HumiditySpec{MassHumidity}()
    ,3500 => HumiditySpec{RelativeHumidity}()
    ,3600 => HumiditySpec{HumidityDewPoint}()
)

#the following functions are essencial in the ordering
#do not change the relative order!!

#equilibria properties are always first
_stateorder(x::PhaseFractions) = 1
_stateorder(x::TwoPhaseEquilibrium) = 2
_stateorder(x::VaporFraction) = 3

#pressure is always first in the scalar properties
_stateorder(x::Pressure) = 4
_stateorder(x::Helmholtz) = 5
_stateorder(x::Gibbs) = 6
_stateorder(x::InternalEnergy) = 7
_stateorder(x::VolumeAmount) = 8
_stateorder(x::Enthalpy) = 9
_stateorder(x::Entropy) = 10
#temperature always for last of the scalar properties
_stateorder(x::Temperature) = 11

#From here, it doesnt matter the order, as those are treated specially

_stateorder(x::HumiditySpec) = 1000

_stateorder(x::MaterialCompounds) = 10002

_stateorder(x::MaterialAmount ) = 10001
_stateorder(x::PhaseTag) = 10000
_stateorder(x::Options) = 10004
_stateorder(x::MolecularWeight) = 10005
_stateorder(x::VariableSpec) = 10006

#using specs here allows to delete all reduce_check_mass(::Type)
function _static_specs(x::Type{ThermodynamicState{S,C}}) where {S,C}
    @nospecialize S,C
    tvec = S.parameters
    res1 = [tvec[i].parameters[1]() for i in 1:length(tvec)]
    if C === Tuple{}
        return res1
    else
        tvec2 = C.parameters
        res2 = [tvec2.parameters[i]() for i in 1:length(tvec)]
        append!(res1,res2)
        return res1
    end
end

@generated function static_amount_type(x::T) where T<: ThermodynamicState
    sps = _static_specs(x)
    val = mapreduce(_reduce_check_mass,+,sps)
    modval = mod(val,MATERIAL_SINGLE_MAX)
    if val in (2100,2200)
        _1 = val
        _2 = val
    else
        _1 = max(val-modval,99)
        _2 = modval
    end
    res =  (AMOUNT_CONST[_1],AMOUNT_CONST[_2])

    return :($res)
end

function _state_from_type(x::Type{T}) where T<: ThermodynamicState
    @nospecialize T
    state_specs =  _static_specs(T)
    has_compound = false
    for spec in state_specs
        if spec isa MaterialCompounds
            has_compound = true
            comp_spec = spec
            break
        end
    end

    #double iteration, used in moist air to create an invalid
    #state where compound and humidity spec are present
    #in this case, the resulting type is MultiPT instead of HumPT
    for spec in state_specs
        if (spec isa HumiditySpec) & !has_compound
            has_compound = true
            comp_spec = spec
            break
        end
    end

    if !has_compound
        comp_spec = SingleComponent()
    end
    sorted_specs = sort(state_specs,by= _stateorder)
    res = (sorted_specs[1],sorted_specs[2],comp_spec)
end

function _static_get_spec_idx(@nospecialize(x::Type{ThermodynamicState{S,C}}),@nospecialize(sp::Type{S2}))::Int where {S,C,S2}
    @nospecialize S,C,S2
    state_specs =  _static_specs(x)
    for (i,spec) in pairs(state_specs)
        if spec isa sp
            return i
        end
    end
    return 0
end

"""
    state_type(x::ThermodynamicState)

Returns a sorted tuple of thermodynamic specifications of the input state,
 the tuple has 3 elements, corresponding to two general specifications and a last amount specification,
  indicating if the state is single or multicomponent.

"""
@generated function state_type(x::ThermodynamicState)
    res = _state_from_type(x)
    return :($res)
end

function get_spec(val::T,st::Tuple) where T<:AbstractSpec
    for spec in st
        if specification(spec) isa T
            return spec
        end
    end
    return nothing
end

function get_spec(val::Type{T},st::Tuple) where T<: AbstractSpec
    for spec in st
        if specification(spec) isa T
            return spec
        end
    end
    return nothing
end

get_spec(val,st::ThermodynamicState) = get_spec(val,st.specs)

function throw_get_spec(val,st::ThermodynamicState)
    return static_get_spec(st,val)
end

function has_spec(val::Type{T},st::ThermodynamicState) where T<: AbstractSpec
    static_has_spec(st,val)
end

function amount_type(st::T) where T<:ThermodynamicState
    return static_amount_type(st)
end

@generated function static_get_spec(st::ThermodynamicState,val::Type{T}) where T<:AbstractSpec
    idx =_static_get_spec_idx(st,T)
    if iszero(idx)
        return :(nothing)
    else
        return :(st.specs[$idx])
    end
end

@generated function static_has_spec(st::ThermodynamicState,val::Type{T}) where T<:AbstractSpec
    idx =_static_get_spec_idx(st,T)
    res = !iszero(idx)
    return :($res)
end

ext: jl


text:
"""
GBD iterations for Ground Structure with Discrete Diameters.

Solver object to store parameter values.
"""
struct IterGBD <: GroundStructureSolver
    addnogoods::Bool
    addcritpath::Bool
    maxiter::Int
    timelimit::Float64 # seconds
    writemodels::Bool
    mastersolver
    subsolver

    function IterGBD(mastersolver, subsolver;
                     addnogoods=false, addcritpath=true, maxiter::Int=100,
                     timelimit=Inf, writemodels=false)
        new(addnogoods, addcritpath, maxiter, timelimit, writemodels,
            mastersolver, subsolver)
    end
end

"Data type for dual solution of subproblem"
struct SubDualSol
    μ::Array{Float64}
    λl::Array{Float64}
    λu::Array{Float64}
end

"""
Data type for master problem.

Variable meanings:
  y: select arc
  z: select diameter for arc
  q: flow through arc
  ϕ: squared flow through arc (ϕ = q²)
"""
struct Master
    model::JuMP.Model
    y::Vector{JuMP.VariableRef}
    z::Array{JuMP.VariableRef, 2}
    q::Vector{JuMP.VariableRef}
    ϕ::Vector{JuMP.VariableRef}
end

"Build model for master problem (ground structure with discrete diameters)."
function make_master(inst::Instance, topo::Topology, optimizer)
    nodes, nnodes = topo.nodes, length(topo.nodes)
    arcs, narcs = topo.arcs, length(topo.arcs)
    terms, nterms = inst.nodes, length(inst.nodes)
    termidx = termindex(nodes, terms)
    ndiams = length(inst.diameters)

    # demand for all nodes, including junctions
    dem = fill(0.0, nnodes)
    for t in 1:nterms
        dem[termidx[t]] = inst.demand[t]
    end

    # adjacency lists
    inarcs, outarcs = [Int[] for v in 1:nnodes], [Int[] for v in 1:nnodes]
    for a in 1:narcs
        tail, head = arcs[a]
        push!(inarcs[head], a)
        push!(outarcs[tail], a)
    end

    # "big-M" bound for flow on arcs
    maxflow = 0.5 * sum(abs.(inst.demand))

    # always use direct mode for SCIP
    model = JuMP.direct_model(MOI.instantiate(optimizer))

    # select arcs from topology with y
    @variable(model, y[1:narcs], Bin)

    # select single diameter for arc with z
    @variable(model, z[1:narcs, 1:ndiams], Bin)

    # flow through arcs
    @variable(model, 0 <= q[1:narcs] <= maxflow)
    @variable(model, 0 <= ϕ[1:narcs] <= maxflow^2)

    # secant cut for ϕ = q²
    @constraint(model, secant[a=1:narcs], ϕ[a] ≤ maxflow*q[a])

    # mass flow balance at nodes
    @constraint(model, flowbalance[v=1:nnodes],
                sum(q[a] for a=inarcs[v]) - sum(q[a] for a=outarcs[v]) == dem[v])

    # allow flow only for active arcs
    @constraint(model, active[a=1:narcs], q[a] <= maxflow*y[a])

    # choose diameter for active arcs
    @constraint(model, choice[a=1:narcs], sum(z[a,d] for d=1:ndiams) == y[a])

    # exclude antiparallel arcs.
    antiidx = antiparallelindex(topo)
    for a in 1:narcs
        if a < antiidx[a] # only at most once per pair
            @constraint(model, y[a] + y[antiidx[a]] ≤ 1)
        end
    end

    L = pipelengths(topo)
    c = [diam.cost for diam in inst.diameters]
    @objective(model, Min, sum(c[i] * L[a] * z[a,i] for a=1:narcs for i=1:ndiams))

    model, y, z, q, ϕ
end

"""
Build model for subproblem (ground structure with discrete diameters).

Corresponds to the domain relaxation with pressure loss overestimation, which
can be turned off via the flag `relaxed`.
"""
function make_sub(inst::Instance, topo::Topology, cand::CandSol, optimizer;
                  relaxed::Bool=true)
    nodes, nnodes = topo.nodes, length(topo.nodes)
    arcs, narcs = topo.arcs, length(topo.arcs)
    terms, nterms = inst.nodes, length(inst.nodes)
    termidx = termindex(nodes, terms)
    ndiams = length(inst.diameters)

    candarcs = filter(a -> any(cand.zsol[a,:]), 1:narcs)
    ncandarcs = length(candarcs)
    tail = [arcs[a].tail for a in candarcs]
    head = [arcs[a].head for a in candarcs]

    model = JuMP.Model(optimizer)

    # unconstrained variable for squared pressure, the bounds are added with
    # inequalities having slack vars.
    # for now, adding variables for all nodes, even if disconnected.
    @variable(model, π[1:nnodes])

    # overestimated pressure loss inequalities
    Dm5 = [diam.value^(-5) for diam in inst.diameters]
    C = inst.ploss_coeff * pipelengths(topo)[candarcs]
    α = C .* cand.qsol[candarcs].^2 .* (cand.zsol[candarcs,:] * Dm5)
    if relaxed
        @constraint(model, ploss[ca=1:ncandarcs],
                    π[tail[ca]] - π[head[ca]] ≥ α[ca])
    else
        @constraint(model, ploss[ca=1:ncandarcs],
                    π[tail[ca]] - π[head[ca]] == α[ca])
    end

    # slack variables to relax the lower and upper bounds for π.
    termlb = [b.lb^2 for b in inst.pressure]
    lb = fill(minimum(termlb), nnodes)
    lb[termidx] = termlb
    termub = [b.ub^2 for b in inst.pressure]
    ub = fill(maximum(termub), nnodes)
    ub[termidx] = termub

    @variable(model, Δl[1:nnodes] ≥ 0)
    @variable(model, Δu[1:nnodes] ≥ 0)

    @constraint(model, pres_lb[v=1:nnodes], π[v] + Δl[v] ≥ lb[v])
    @constraint(model, pres_ub[v=1:nnodes], π[v] - Δu[v] ≤ ub[v])

    @objective(model, Min, sum(Δl[v] + Δu[v] for v=1:nnodes))

    model, π, Δl, Δu, ploss, pres_lb, pres_ub
end

"Add tangent cut to quadratic inequality (ϕ ≥ q^2) if violated."
function quadratic_tangent(model, q, ϕ, qsol, ϕsol; cb=cb)
    violated = ϕsol < qsol^2 - ɛ
    if !violated
        return 0
    end

    # add 1st order Taylor approx:
    # q^2 ≈ 2*qsol*(q - qsol) + qsol^2 = 2*qsol*q - qsol^2
    if cb === nothing
        @constraint(model, ϕ ≥ 2*qsol*q - qsol^2)
    else
        @cb_constraint(cb, ϕ ≥ 2*qsol*q - qsol^2)
    end

    return 1
end

"""
Find tightest linear overestimator for given values v.

It's of the form: sum_i a_i x_i + sum_j b_j y_j + c
where x_i and y_j take values 0 or 1, with exactly one term active per sum.

The coefficients a, b, c are returned.
"""
function linear_overest(values::Matrix{Float64}, cand_i::Int, cand_j::Int,
                        optimizer)
    m, n = size(values)

    @assert 1 <= cand_i <= m
    @assert 1 <= cand_j <= n

    model = JuMP.Model(optimizer)

    # coefficients to be found
    @variable(model, a[1:m])
    @variable(model, b[1:n])
    @variable(model, c)

    # slack to minimize
    @variable(model, t[1:m,1:n] ≥ 0)

    # overestimate at all points
    @constraint(model, overest[i=1:m, j=1:n], a[i] + b[j] + c == values[i,j] + t[i,j])

    # be exact at candidate solution
    @constraint(model, t[cand_i, cand_j] == 0)

    # be as tight as possible everywhere else
    @objective(model, Min, sum(t[i,j] for i=1:m for j=1:n))

    # solve it
    JuMP.optimize!(model)
    status = JuMP.termination_status(model)
    @assert status == MOI.OPTIMAL

    JuMP.value.(a), JuMP.value.(b), JuMP.value(c)
end

"Linearized & reformulated cut based on single path."
function pathcut(inst::Instance, topo::Topology, master::Master, cand::CandSol,
                 path::Vector{Arc}, solver; cb=nothing)
    aidx = arcindex(topo)
    pathidx = [aidx[arc] for arc in path]
    npath = length(path)
    ndiam = length(inst.diameters)
    zsol = cand.zsol[pathidx,:] # sparse solution
    qsol = cand.qsol[pathidx,:] # sparse solution
    D = [diam.value for diam in inst.diameters]

    # set loose pressure bounds for non-terminal nodes
    terminals = indexin(inst.nodes, topo.nodes)
    nnodes = length(topo.nodes)
    πlb_min = minimum([b.lb for b in inst.pressure])^2
    πub_max = maximum([b.ub for b in inst.pressure])^2
    πlb = fill(πlb_min, nnodes)
    πub = fill(πub_max, nnodes)
    πlb[terminals] = [b.lb^2 for b in inst.pressure]
    πub[terminals] = [b.ub^2 for b in inst.pressure]

    # coefficients of z in supremum expression
    ν = zsol * D.^(-5)
    @assert length(ν) == npath
    β = ν * (D.^5)'
    @assert size(β) == (npath, ndiam)
    @assert all(β .> 0)

    # linearize the supremum, build up dense coefficient matrix
    coeffs = fill(0.0, (npath, ndiam))
    offset = 0.0
    @assert size(coeffs) == size(β)

    # - tail of path
    tail = path[1].tail
    coeffs[1,:] += πub[tail] * β[1,:]

    # - intermediate nodes, with in- and outgoing arcs
    for v in 2:npath
        uv, vw = path[v-1], path[v]
        @assert uv.head == vw.tail
        node = uv.head

        # prepare all possible values that need to be overestimated
        supvalues = zeros(ndiam + 1, ndiam + 1)
        # the case where no diameter is selected on the first arc
        supvalues[1, 2:end] = β[v, :] * πub[node]
        # similarly for no diameter on the second arc
        supvalues[2:end, 1] = - β[v-1, :] * πlb[node]
        # finally, when both arcs are active
        β1st, β2nd = β[v-1, :], β[v, :]
        @assert size(β1st) == (ndiam,)
        @assert size(β2nd) == (ndiam,)
        βdiff = - repeat(β1st, 1, ndiam) + repeat(β2nd', ndiam, 1)
        @assert size(βdiff) == (ndiam, ndiam)
        supvalues[2:end, 2:end] = max.(βdiff * πub[node], βdiff * πlb[node])

        # the current values should be met exactly
        cand_i = findfirst(zsol[v-1,:])
        cand_j = findfirst(zsol[v,:])
        @assert cand_i ≠ 0 && cand_j ≠ 0

        # get coeffs of overestimation, assuming aux vars z_uv,0 and z_vw,0
        fix_i, fix_j = cand_i + 1, cand_j + 1
        cuv, cvw, c = linear_overest(supvalues, fix_i, fix_j, solver)

        # need to transform the coefficients to remove aux vars
        coeffs[v-1,:] += cuv[2:end] .- cuv[1]
        coeffs[v,:]   += cvw[2:end] .- cvw[1]
        offset += c + cuv[1] + cvw[1]
    end

    # - head of path
    head = path[end].head
    coeffs[end,:] += -1 * πlb[head] * β[end,:]

    # coefficients of ϕ
    C = inst.ploss_coeff * pipelengths(topo)[pathidx]
    α = ν .* C

    # add cut:  coeffs * z + offset ≥ α * ϕ
    z = master.z[pathidx,:]
    ϕ = master.ϕ[pathidx]

    if cb === nothing
        @constraint(
            master.model,
            sum(coeffs[a,i]*z[a,i] for a=1:npath for i=1:ndiam) + offset
            ≥ sum(α[a]*ϕ[a] for a=1:npath))
    else
        @cb_constraint(
            cb,
            sum(coeffs[a,i]*z[a,i] for a=1:npath for i=1:ndiam) + offset
            ≥ sum(α[a]*ϕ[a] for a=1:npath))
    end

    return 1
end

"Linearized & reformulated cuts based on critical paths."
function critpathcuts(inst::Instance, topo::Topology, master::Master,
                      cand::CandSol, sub::SubDualSol, solver; cb=nothing)
    ncuts = 0

    # compute dense dual flow
    narcs = length(topo.arcs)
    dualflow = fill(0.0, narcs)
    actarcs = vec(sum(cand.zsol, dims=2) .> 0)
    dualflow[actarcs] = sub.μ

    paths, pathflow = flow_path_decomp(topo, dualflow)

    # TODO: test that cuts are valid and separating
    for path in paths
        ncuts += pathcut(inst, topo, master, cand, path, solver, cb=cb)
    end

    for aidx in 1:narcs
        ncuts += quadratic_tangent(master.model, master.q[aidx], master.ϕ[aidx],
                                   cand.qsol[aidx], cand.ϕsol[aidx], cb=cb)
    end

    return ncuts
end

"Construct all Benders cuts from the solution of a subproblem."
function cuts(inst::Instance, topo::Topology, master::Master, cand::CandSol,
              sub::SubDualSol, solver; addnogoods=true, addcritpath=true)
    @assert any([addnogoods, addcritpath]) # must cut off!
    ncuts = 0
    if addnogoods
        ncuts += nogood(master.model, master.z, cand.zsol)
    end
    if addcritpath
        ncuts += critpathcuts(inst, topo, master, cand, sub, solver)
    end
    ncuts
end

"Iteration based implementation of GBD."
function PipeLayout.optimize(inst::Instance, topo::Topology, solver::IterGBD)
    run_gbd(inst, topo, solver.mastersolver, solver.subsolver,
            maxiter=solver.maxiter, timelimit=solver.timelimit,
            addnogoods=solver.addnogoods, addcritpath=solver.addcritpath,
            writemodels=solver.writemodels)
end

function run_gbd(inst::Instance, topo::Topology, mastersolver, subsolver;
                 maxiter::Int=100, timelimit=Inf, addnogoods=false,
                 addcritpath=true, writemodels=false)
    finaltime = time() + timelimit

    # initialize
    master = Master(make_master(inst, topo, mastersolver)...)
    dual, status = 0.0, MOI.OPTIMIZE_NOT_CALLED

    for iter=1:maxiter
        if !stilltime(finaltime)
            status = MOI.TIME_LIMIT
            @debug "Timelimit reached."
            return Result(status, nothing, Inf, dual, iter)
        end

        @debug "Iter $(iter)"

        # resolve (relaxed) master problem, build candidate solution
        writemodels && writeLP(master.model, "master_iter$(iter).lp", genericnames=false)
        settimelimit!(master.model, mastersolver, finaltime - time())
        JuMP.optimize!(master.model)
        status = JuMP.termination_status(master.model)
        if status == MOI.INFEASIBLE
            @debug "  relaxed master is infeasible :-("
            return Result(status, nothing, Inf, Inf, iter)
        elseif status in [MOI.TIME_LIMIT, MOI.NODE_LIMIT]
            return Result(status, nothing, Inf, dual, iter)
        elseif status == MOI.OPTIMAL
            # good, we continue below
        else
            error("Unexpected status: $(status)")
        end

        zsol = JuMP.value.(master.z)
        cand = CandSol(zsol .>= 0.5, JuMP.value.(master.q), JuMP.value.(master.ϕ))

        dual = JuMP.objective_value(master.model)
        @debug begin
            "  dual bound: $(dual)"
            "  cand. sol:$(Tuple.(findall(!iszero, cand.zsol)))"
        end

        # check whether candidate has tree topology
        ysol = JuMP.value.(master.y)
        candtopo = topology_from_candsol(topo, ysol)
        if !is_tree(candtopo)
            fullcandtopo = topology_from_candsol(topo, ysol, true)
            cycle = find_cycle(fullcandtopo)
            if length(cycle) == 0
                # TODO: Actually, this might be optimal, but it could also occur
                # when adding some irrelevant pipe is cheaper than increasing
                # the diameter. How to distinguish these cases?
                nogood(master.model, master.y, ysol)
                @debug "  skip disconnected topology with nogood."
            else
                avoid_topo_cut(master.model, master.y, topo, cycle)
                @debug "  skip non-tree topology, cycle: $(cycle)"
            end
            continue
        end

        # solve subproblem (from scratch, no warmstart)
        submodel, π, Δl, Δu, ploss, plb, pub = make_sub(inst, topo, cand, subsolver)
        writemodels && writeLP(submodel, "sub_relax_iter$(iter).lp", genericnames=false)
        settimelimit!(submodel, subsolver, finaltime - time())
        JuMP.optimize!(submodel)
        substatus = JuMP.termination_status(submodel)
        @assert substatus == MOI.OPTIMAL "Slack model is always feasible"
        totalslack = JuMP.objective_value(submodel)
        if totalslack ≈ 0.0
            # maybe only the relaxation is feasible, we have to check also the
            # "exact" subproblem with equations constraints.
            submodel2, _ = make_sub(inst, topo, cand, subsolver, relaxed=false)
            writemodels && writeLP(submodel2, "sub_exact_iter$(iter).lp", genericnames=false)
            settimelimit!(submodel2, subsolver, finaltime - time())
            JuMP.optimize!(submodel2)
            substatus2 = JuMP.termination_status(submodel2)
            @assert substatus2 == MOI.OPTIMAL "Slack model is always feasible"
            totalslack2 = JuMP.objective_value(submodel2)

            if totalslack2 ≈ 0.0
                @debug "  found feasible solution :-)"
                primal = JuMP.objective_value(master.model)
                return Result(substatus2, cand, primal, dual, iter)
            else
                # cut off candidate with no-good on z
                @debug "  subproblem/relaxation gap!"
                nogood(master.model, master.z, cand.zsol)
                continue
            end
        end

        dualsol = SubDualSol(JuMP.dual.(ploss), JuMP.dual.(plb), JuMP.dual.(pub))

        # generate cuts and add to master
        ncuts = cuts(inst, topo, master, cand, dualsol, subsolver,
                     addnogoods=addnogoods, addcritpath=addcritpath)
        @debug "  added $(ncuts) cuts."
    end

    Result(MOI.ITERATION_LIMIT, nothing, Inf, dual, maxiter)
end


text:
@testset "Midway Histogram Equalization" begin
    for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64})
        # Define two images consisting of luminance gradients with different
        # (non-overlapping) intensities. The "midway" equalization should produce
        # two new images for which the cummulative distribution function of the
        # image intensities is identical.
        img1 = T.(zeros(256,10))
        for i = 32:96
            for j = 1:10
                img1[i,j] = i/255.0
            end
        end
        img2 = T.(zeros(256,10))
        for i = 128:192
            for j = 1:10
                img2[i,j] = i/255.0
            end
        end
        img1o, img2o = adjust_histogram([img1, img2], MidwayEqualization(nbins = 256))
        edges1, counts1 = build_histogram(img1o, 256, minval = 0, maxval = 1)
        edges2, counts2 = build_histogram(img2o, 256, minval = 0, maxval = 1)
        @test sum(cumsum(counts2) - cumsum(counts1)) == 0

        edges, _ = build_histogram(img1, 256, minval = 0, maxval = 1)
        img1o, img2o = adjust_histogram([img1, img2], MidwayEqualization(edges = edges))
        edges1, counts1 = build_histogram(img1o, 256, minval = 0, maxval = 1)
        edges2, counts2 = build_histogram(img2o, 256, minval = 0, maxval = 1)
        @test sum(cumsum(counts2) - cumsum(counts1)) == 0

    end

    for T in (RGB{N0f8}, RGB{N0f16}, RGB{Float32}, RGB{Float64})
        imgg1 = zeros(Gray{Float32},256,10)
        for i = 32:96
            for j = 1:10
                imgg1[i,j] = i/255.0
            end
        end
        img1 = colorview(RGB,imgg1,imgg1,imgg1)
        img1 = T.(img1)

        imgg2 = zeros(Gray{Float32},256,10)
        for i = 128:192
            for j = 1:10
                imgg2[i,j] = i/255.0
            end
        end
        img2 = colorview(RGB,imgg2,imgg2,imgg2)
        img2 = T.(img2)

        img1o, img2o = ImageContrastAdjustment.adjust_histogram([img1, img2], ImageContrastAdjustment.MidwayEqualization(nbins = 256))

        edges1, counts1 = ImageContrastAdjustment.build_histogram(img1o, 256, minval = 0, maxval = 1)
        edges2, counts2 = ImageContrastAdjustment.build_histogram(img2o, 256, minval = 0, maxval = 1)
        # The algorithm equalizes the Y channels from the YIQ decomposition of the RGB images and then
        # constructs new RGB images by combining the equalised Y channels with the IQ components.
        # The build_histogram function then implicitly converts the "midway" RGB images to Gray
        # in order to construct the histogram. After this process the cummulative distribution functions
        # of these luminance gradients are no longer identical but still close.
        @test abs(sum(cumsum(counts2) - cumsum(counts1))) <= 20

        edges, _ = build_histogram(img1, 256, minval = 0, maxval = 1)
        img1o, img2o = adjust_histogram([img1, img2], MidwayEqualization(edges = edges))
        edges1, counts1 = build_histogram(img1o, 256, minval = 0, maxval = 1)
        edges2, counts2 = build_histogram(img2o, 256, minval = 0, maxval = 1)
        @test sum(cumsum(counts2) - cumsum(counts1)) == 0
    end
end


text:
type History{I,V} <: UnivalueHistory{I}
    lastiter::I
    iterations::Vector{I}
    values::Vector{V}

    function History(::Type{V}, ::Type{I})
        new(typemin(I), Array(I, 0), Array(V, 0))
    end
end

function History{I,V}(v::Type{V}, i::Type{I} = Int)
    History{I,V}(v, i)
end

Base.length(history::History) = length(history.iterations)
Base.enumerate(history::History) = zip(history.iterations, history.values)
Base.first(history::History) = history.iterations[1], history.values[1]
Base.last(history::History) = history.iterations[end], history.values[end]
Base.get(history::History) = history.iterations, history.values

function Base.push!{I,V}(
        history::History{I,V},
        iteration::I,
        value::V)
    lastiter = history.lastiter
    iteration > lastiter || throw(ArgumentError("Iterations must increase over time"))
    history.lastiter = iteration
    push!(history.iterations, iteration)
    push!(history.values, value)
    value
end

function Base.push!{I,V}(
        history::History{I,V},
        value::V)
    lastiter = history.lastiter == typemin(I) ? zero(I) : history.lastiter
    iteration = lastiter + one(history.lastiter)
    history.lastiter = iteration
    push!(history.iterations, iteration)
    push!(history.values, value)
    value
end

Base.print{I,V}(io::IO, history::History{I,V}) = print(io, "$(length(history)) elements {$I,$V}")

function Base.show{I,V}(io::IO, history::History{I,V})
    println(io, "History")
    println(io, "  * types: $I, $V")
    print(io,   "  * length: $(length(history))")
end


text:
using AccelerationBenchmark
using JLD, DataFrames,  StatPlots
pyplot()

searchdir(path,key) = filter(x->contains(x,key), readdir(path))

"Performance profiles from the Moré tests"
function more_tests()
    mdf = reduce(append!, [load(fname, "mdf") for fname in ["data/"*fname for fname in searchdir("data", "jld")]])
    rdf = AccelerationBenchmark.createratiodataframe(mdf)

    # TODO: GroupedErrors is using a uniform spread in the x-values (\tau-values),
    # so there is quite low resolution for small x-values.
    # I need to deal with this somehow, and combine that with a downsampling approach.
    # Alternatively, we can create fevals_AG.csv and run it through "mergealldata.jl" with BenchmarkProfiles.

    syms = [:fcalls, :gcalls, :Iterations, :CPUtime]
    for sym in syms
        rdf[sym] .=  NaNMath.log2.(rdf[sym])
    end
    maxfun(x) = 1.1maximum(x)

    pltf = AccelerationBenchmark.createperfprofiles(rdf, :fcalls; xlabel="f-calls, log₂-scale",
                                                    xscale=:identity, maxfun = maxfun)
    plti = AccelerationBenchmark.createperfprofiles(rdf, :Iterations; xlabel="Iterations, log₂-scale",
                                                    xscale=:identity, maxfun = maxfun)
    pltg = AccelerationBenchmark.createperfprofiles(rdf, :gcalls; xlabel="g-calls, log₂-scale",
                                                    xscale=:identity, maxfun = maxfun)
    pltc = AccelerationBenchmark.createperfprofiles(rdf, :CPUtime; xlabel="CPU time, log₂-scale",
                                                    xscale=:identity, maxfun = maxfun)

    xlims!(pltf,0,5);xlims!(pltg,0,5);xlims!(pltc,0,5);xlims!(plti,0,5)
    plt = plot(pltf,pltg,pltc,plti)
end

function cutest_tests(fname::AbstractString = "data/cutest_0_1000.jld",
                      xmax::Number = NaN)
    mdf = load(fname, "mdf")

    rdf = AccelerationBenchmark.createratiodataframe(mdf)

    # TODO: GroupedErrors is using a uniform spread in the x-values (\tau-values),
    # so there is quite low resolution for small x-values.
    # I need to deal with this somehow, and combine that with a downsampling approach.
    syms = [:fcalls, :gcalls, :Iterations, :CPUtime]
    for sym in syms
        rdf[sym] .=  NaNMath.log2.(rdf[sym])
    end
    maxfun(x) = 1.1maximum(x)

    pltf = AccelerationBenchmark.createperfprofiles(rdf, :fcalls; xlabel="f-calls, log₂-scale",
                                                    xscale=:identity, maxfun = maxfun)
    plti = AccelerationBenchmark.createperfprofiles(rdf, :Iterations; xlabel="Iterations, log₂-scale",
                                                    xscale=:identity, maxfun = maxfun)
    pltg = AccelerationBenchmark.createperfprofiles(rdf, :gcalls; xlabel="g-calls, log₂-scale",
                                                    xscale=:identity, maxfun = maxfun)
    pltc = AccelerationBenchmark.createperfprofiles(rdf, :CPUtime; xlabel="CPU time, log₂-scale",
                                                    xscale=:identity, maxfun = maxfun)

    if isfinite(xmax)
        xlims!(pltf,0,xmax);xlims!(pltg,0,xmax);xlims!(pltc,0,xmax);xlims!(plti,0,xmax)
    end
    plt = plot(pltf,pltg,pltc,plti)
end

plt = more_tests()
savefig(plt, "perfprof_more_tests.svg")
savefig(plt, "perfprof_more_tests.png")

minvar = 0; maxvar = 1000
pltc = cutest_tests("data/cutest_$(minvar)_$(maxvar).jld")
savefig(pltc, "perfprof_cutest_$(minvar)_$(maxvar).svg")
savefig(pltc, "perfprof_cutest_$(minvar)_$(maxvar).png")


text:
using Documenter
using LightLattices

makedocs(
    sitename = "LightLattices.jl",
    format = Documenter.HTML(; prettyurls = get(ENV, "CI", nothing) == "true"),
    modules = [LightLattices],
    authors = "Grigory Starkov",
    pages = ["Contents" => "index.md",
             "Manual" => "manual.md",
             "Examples" => "examples.md",
             "Index" => "list.md"
            ]
)

deploydocs(
    repo = "github.com/Gregstrq/LightLattices.jl.git",
    devbranch = "main"
)


text:
using Test
using Forecast

@testset "hma" begin
    r_data = rand(24)

    # asserting return of correct weights
    # as listed here: https://www.mathworks.com/help/econ/seasonal-adjustment-using-snxd7m-seasonal-filters.html
    x = hmaSymmetricWeights(13)
    @test length(x) == 13
    @test x         == Float16.([-0.019, -0.028, 0.0, .066, .147, .214, .24, .214, .147, .066, 0.0, -0.028, -0.019])

    # asserting correct application of the 13-term moving average
    x = hma(r_data, 13)
    @test length(x)             == 24
    @test first(x)              == 0.5508635070863032
    @test last(x)               == 0.5843536698022159
    @test x[length(x)//2|>Int]  == 0.32899917067680573
end


text:
using LazySets: _isapprox

@testset "ORBIT algorithm" begin

    # --------------------------------------------------------------------------
    # harmonic oscillator x'' = -x, s.t. x(0) = 1, x'(0) = 0
    # --------------------------------------------------------------------------
    prob, dt = harmonic_oscillator();
    δ = 1e-2
    sol = solve(prob, tspan=(0.0, 20.0), alg=ORBIT(δ=δ));

    # check result
    @test sum((set(sol[k]).element[1] - cos(δ * (k-1)))^2 for k in 1:length(sol)) < eps(Float64)
    @test sum((set(sol[k]).element[2] + sin(δ * (k-1)))^2 for k in 1:length(sol)) < eps(Float64)

    # test number of steps
    sol = solve(prob, tspan=(0.0, 3δ), alg=ORBIT(δ=δ));
    @test _isapprox(tspan(sol), interval(0, 3δ))
    x = sol[end]

    @test _isapprox(tspan(x), interval(3δ))
    @test element(set(x)) ≈ [cos(3δ), -sin(3δ)]

    # test time span sequence
    @test all(_isapprox(tspan(sol[i]), interval(δ*(i-1))) for i in eachindex(sol))

    # --------------------------------------------------------------------------
    # forced oscillator Mx'' + Kx = R, s.t. x(0) = [0, 0], x'(0) = [0, 0]
    # --------------------------------------------------------------------------
    prob, dt = forced_oscillator();
    δ = 1e-2
    sol = solve(prob, tspan=(0.0, 20.0), alg=ORBIT(δ=δ));

    # check result
    U = forced_oscillator_solution()
    U1 = k -> U((k-1) * δ)[1]
    U2 = k -> U((k-1) * δ)[2]
    @test sum((set(sol[k]).element[1] - U1(k))^2 for k in 1:length(sol)) < eps(Float64)
    @test sum((set(sol[k]).element[2] - U2(k))^2 for k in 1:length(sol)) < eps(Float64)
end


text:
using JLD2, FileIO, NetCDF
using Elefridge, ZfpCompression
using TranscodingStreams,CodecZstd

ZstdCompressorL3 = ZstdCompressor(level=3)
TranscodingStreams.initialize(ZstdCompressorL3)

ZstdCompressorL10 = ZstdCompressor(level=10)
TranscodingStreams.initialize(ZstdCompressorL10)

ZstdCompressorL22 = ZstdCompressor(level=22)
TranscodingStreams.initialize(ZstdCompressorL22)

path = "/network/aopp/chaos/pred/kloewer/esowc/"
member1files = filter(x->endswith(x,".nc"),readdir(joinpath(path,"member1")))
steps = [parse(Int,split(split(file,"step")[2],".")[1]) for file in member1files]
sort!(steps)

# parameter
Nens = 25
Nbits = 32
Nlon = 1800
Nlat = 901
Nvert = 91
Ntsteps = length(steps)
chunk_size = 8
nchunks = Nvert ÷ chunk_size
keepbits = [8,7,6,5]

# preallocate
Xzfp,Xzstd = try load(joinpath(path,"analysis/ensemble_compression_$(Nens)members_s5-8.jld2"),"Xzfp","Xzstd") catch; zeros(Float64,Ntsteps,length(keepbits)),zeros(Float64,Ntsteps,length(keepbits)) end

for t in 13:Ntsteps
    println("Time step $(steps[t])h")
    temp = fill(0f0,Nens,Nlon,Nlat,Nvert)

    print("Loading ensemble member ")
    # read the data
    for ie in 1:Nens
        print("$ie,")
        ncfile = NetCDF.open(joinpath(path,"member$ie","ensemble.t.member$ie.step$(steps[t]).ll.nc"))
        temp[ie,:,:,:] = ncfile.vars["t"][1:Nlon,1:Nlat,1:Nvert]
    end
    println("\nCompressing...")

    # ZFP COMPRESSION
    for (ik,keepbit) in enumerate(keepbits)

        precision = keepbit+7   # For 4D compression

        # ZFP: compress in chunks
        compressed_size = 0

        for ichunk in 1:nchunks-1    # compress the last one after loop
            print("$ichunk,")
            start = (ichunk-1)*chunk_size+1
            stop = start+chunk_size-1
            temp_chunk = temp[:,:,:,start:stop]
            compressed_size += sizeof(zfp_compress(temp_chunk,precision=precision))
        end

        # compress last chunk
        start = (nchunks-1)*chunk_size+1
        temp_chunk = temp[:,:,:,start:end]
        compressed_size += sizeof(zfp_compress(temp_chunk,precision=precision))

        Xzfp[t,ik] = 2*sizeof(temp)/compressed_size
        println("$(keepbit)bits, CF ZFP: $(Xzfp[t,ik])")
        save(joinpath(path,"analysis/ensemble_compression_$(Nens)members_s5-8.jld2"),"Xzfp",Xzfp,"Xzstd",Xzstd)
    end

    # ZSTD COMPRESSION
    for (ik,keepbit) in enumerate(keepbits)
        round!(temp,keepbit)
        Xr8 = copy(reinterpret(UInt8,vec(temp)))
        Xzstd[t,ik] = 2*sizeof(Xr8)/sizeof(transcode(ZstdCompressorL10,Xr8))
        println("$(keepbit)bits, CF ZSTD: $(Xzstd[t,ik])")
        save(joinpath(path,"analysis/ensemble_compression_$(Nens)members_s5-8.jld2"),"Xzfp",Xzfp,"Xzstd",Xzstd)
    end

    println("Compression done.")
    println("---")
end


text:
import MitosisStochasticDiffEq as MSDE
using Test, Random
using StaticArrays
using Statistics
using LinearAlgebra

@testset "static array tests" begin
seed = 1234
Random.seed!(seed)

# define SDE function
f(u,p,t) = p[1]*u .+ p[2]
g(u,p,t) = p[3]

function fstat(u,p,t)
  dx = p[1]*u[1] + p[2]
  @SVector [dx]
end

function gstat(u,p,t)
  @SVector [p[3]]
end

# time span
tstart = 0.0
tend = 1.0
dt = 0.01
trange = tstart:dt:tend

# initial condition
u0 = [1.1]
u0stat = @SVector [1.1]

# set true model parameters
par = [-0.2, 0.1, 0.9]

# set of linear parameters Eq.~(2.2)
plin = copy(par)

sdekernel1 = MSDE.SDEKernel(f,g,trange,par)
sdekernel2 = MSDE.SDEKernel(fstat,gstat,trange,par)

Random.seed!(seed)
samples1 = MSDE.sample(sdekernel1, u0, save_noise=false)
Random.seed!(seed)
samples2 = MSDE.sample(sdekernel2, u0stat, save_noise=false)

@test isapprox(samples1[1], samples2[1], rtol=1e-10)
@test isapprox(samples1[2][3], samples2[2][3], rtol=1e-10)
@test typeof(samples2[2][2][end]) <: SArray

# initial values for ODE
mynames = (:logscale, :μ, :Σ);
myvalues1 =  [0.0, 1.5, 0.1];
NT1 = NamedTuple{mynames}(myvalues1)

#myvalues2 = @SVector [0.0, 1.5, 0.1] # will be just converted to a standard array
myvalues2 = [0.0, (@SVector [1.5]), (@SMatrix[0.1])]
NT2 = NamedTuple{mynames}(myvalues2)

sdetildekernel1 = MSDE.SDEKernel(f,g,trange,plin)
sdetildekernel2 = MSDE.SDEKernel(fstat,gstat,trange,plin)
message1, backward1 = MSDE.backwardfilter(sdetildekernel1, NT1)
message2, backward2 = MSDE.backwardfilter(sdetildekernel2, NT2)

@test isapprox(backward1, backward2, rtol=1e-10)
@test typeof(backward2.x[1]) <: SArray
@test typeof(backward2.x[2]) <: SArray
@test typeof(backward2.x[3]) <: SArray

x0 = [1.34]
x0stat = @SVector [1.34]
ll0 = randn()

Random.seed!(seed)
samples1 = MSDE.forwardguiding(sdekernel1, message1,
  (x0, ll0), inplace=true)
Random.seed!(seed)
samples2 = MSDE.forwardguiding(sdekernel2, message2,
  (x0stat, ll0), inplace=false)

@test isapprox(samples1[2][1], samples2[2][1], rtol=1e-10)
@test isapprox(samples1[2][2], samples2[2][2], rtol=1e-10)
@test_broken typeof(samples2[2][end]) <: SArray
@test_broken typeof(samples2[3]) <: SArray
end

@testset "static tilde parameter tests" begin
  d = 2

  𝕏 = SVector{d,Float64}
  𝕏_ = SVector{d+1,Float64}
  bθ = @SVector [.5, 0.9]
  σ0 = @SVector [1.25, 1.35]
  θ0 = (bθ, σ0)

  M = @SMatrix [-1.0 1.0; 1.0 -1.0]
  f(u,θ,t) = 𝕏(tanh.(Diagonal(θ[1]) * M * u))  # f(u,θ,t) = Diagonal(θ[1]) * M * u
  g(u,θ,t) = θ[2]

  dt0 = 0.001
  trange = 0.0:dt0:1.0

  u0 = zero(𝕏)  # value at root node

  B(θ) = (Diagonal(θ[1]) * M)
  Σ(θ) = SMatrix{2,2}(Diagonal(θ[2]))
  b̃ = @SVector [.1, 0.1]  #θ̃ = θ0
  θ̃ = (b̃, σ0)

  θlinstat = ((B(θ̃), @SVector(zeros(d)), Σ(θ̃)))
  θlin =  Array.(θlinstat)
  @show typeof(θlin), typeof(θlinstat)

  # initial values for ODE
  Random.seed!(123)
  logscale = randn()
  ν = @SVector randn(d)
  P = @SMatrix randn(d,d)
  myvalues = [logscale, ν, P]
  mynames = (:logscale, :μ, :Σ)
  NTstat = NamedTuple{mynames}(myvalues)
  NT = NamedTuple{mynames}([logscale, Array(ν), Array(P)])

  tildekernelstat = MSDE.SDEKernel(f,g,trange,θlinstat)
  messagestat, backwardstat = MSDE.backwardfilter(tildekernelstat, NTstat)

  @test typeof(backwardstat.x[1]) <: SArray
  @test typeof(backwardstat.x[2]) <: SArray
  @test typeof(backwardstat.x[3]) <: SArray

  tildekernel = MSDE.SDEKernel(f,g,trange,θlin)
  message, backward = MSDE.backwardfilter(tildekernel, NT)

  @test !(typeof(backward.x[1]) <: SArray)
  @test !(typeof(backward.x[2]) <: SArray)
  @test !(typeof(backward.x[3]) <: SArray)

  @test isapprox(backwardstat.x[1], backward.x[1], rtol=1e-10)
  @test isapprox(backwardstat.x[2], backward.x[2], rtol=1e-10)
  @test isapprox(backwardstat.x[3], backward.x[3], rtol=1e-10)

  # mixed cases
  message1, backward1 = MSDE.backwardfilter(tildekernel, NTstat)
  @test typeof(backward1.x[1]) <: SArray
  @test typeof(backward1.x[2]) <: SArray
  @test typeof(backward1.x[3]) <: SArray

  @test isapprox(backward1.x[1], backward.x[1], rtol=1e-10)
  @test isapprox(backward1.x[2], backward.x[2], rtol=1e-10)
  @test isapprox(backward1.x[3], backward.x[3], rtol=1e-10)

  message2, backward2 = MSDE.backwardfilter(tildekernelstat, NT)
  @test !(typeof(backward2.x[1]) <: SArray)
  @test !(typeof(backward2.x[2]) <: SArray)
  @test !(typeof(backward2.x[3]) <: SArray)

  @test isapprox(backward2.x[1], backward.x[1], rtol=1e-10)
  @test isapprox(backward2.x[2], backward.x[2], rtol=1e-10)
  @test isapprox(backward2.x[3], backward.x[3], rtol=1e-10)

  # backward
  (c, μ, Pmat) = NT
  p = tildekernel.p
  u0 = MSDE.mypack(μ, Pmat, c)
  MSDE.CovarianceFilterODE(u0, p, 0.0)
  b, β, σtil = p
  MSDE.compute_dν(b,μ,β)
  MSDE.compute_dP(b,Pmat,σtil)
  @test typeof(MSDE.compute_dν(b,μ,β)) == typeof(μ)
  @test typeof(MSDE.compute_dP(b,Pmat,σtil)) == typeof(Pmat)

  # backwardstat
  (c, μ, Pmat) = NTstat
  p = tildekernelstat.p
  u0 = MSDE.mypack(μ, Pmat, c)
  MSDE.CovarianceFilterODE(u0, p, 0.0)
  b, β, σtil = p
  MSDE.compute_dν(b,μ,β)
  MSDE.compute_dP(b,Pmat,σtil)
  @test typeof(MSDE.compute_dν(b,μ,β)) == typeof(μ)
  @test typeof(MSDE.compute_dP(b,Pmat,σtil)) == typeof(Pmat)

  # backward1
  (c, μ, Pmat) = NTstat
  p = tildekernel.p
  u0 = MSDE.mypack(μ, Pmat, c)
  MSDE.CovarianceFilterODE(u0, p, 0.0)
  b, β, σtil = p
  MSDE.compute_dν(b,μ,β)
  MSDE.compute_dP(b,Pmat,σtil)
  @test typeof(MSDE.compute_dν(b,μ,β)) == typeof(μ)
  @test typeof(MSDE.compute_dP(b,Pmat,σtil)) == typeof(Pmat)

  # backward2
  (c, μ, Pmat) = NT
  p = tildekernelstat.p
  u0 = MSDE.mypack(μ, Pmat, c)
  MSDE.CovarianceFilterODE(u0, p, 0.0)
  b, β, σtil = p
  @test typeof(MSDE.compute_dν(b,μ,β)) == typeof(μ)
  @test typeof(MSDE.compute_dP(b,Pmat,σtil)) == typeof(Pmat)
end


text:
"""
    KNNKernel
Abstract super type for all weighting kernels
"""
abstract type KNNKernel end

sort_idxs(::KNNKernel) = false

"""
    list_kernels()
Lists all implemented KNN weighting kernels
"""
list_kernels() = subtypes(KNNKernel)

####
# User Defined Kernel (UDK)
####

_nothing(x) = nothing

"""
    UserDefinedKernel(;func::Function = x->nothing, sort::Bool=false)

Wrap a user defined nearest neighbors weighting function `func` as a `KNNKernel`.

# Keywords

- `func` : user-defined nearest neighbors weighting function. The function
   should have the signature `func(dists_matrix)::Union{Nothing, <:AbstractMatrix}`.
   The `dists_matrix` is a `n` by `K` nearest neighbors distances matrix where
   `n` is the number of samples in the test dataset and `K` is number of neighbors.
   `func` should either output `nothing` or an `AbstractMatrix` of the same shape
   as `dists_matrix`.  If `func(dists_matrix)` returns nothing then all k-nearest
   neighbors in each row are assign equal weights.
- `sort` : if true requests that the `dists_matrix` be sorted before being passed
   to `func`. The sort is done in a manner that puts the k-nearest neighbors in
   each row  of `dists_matrix` in acesending order .

"""
struct UserDefinedKernel{T<:Function} <: KNNKernel
   func::T
   sort::Bool
   UserDefinedKernel(func::F, sort::Bool) where {F<:Function} = new{F}(func, sort)
end

"""
    UDK
Alias for `UserDefinedKernel`
"""
const UDK = UserDefinedKernel

UserDefinedKernel(;func= _nothing, sort=false) = UserDefinedKernel(func, sort)
sort_idxs(kernel::UserDefinedKernel) = kernel.sort

function get_weights(kernel::UserDefinedKernel, dists_matrix)
    func = kernel.func
    nsamples, K = size(dists_matrix)
    weights = func(dists_matrix)
    weights === nothing && return nothing
    check_onebased_indexing("weights", weights)
    (nsamples, K) == size(weights) || throw(
        DimensionMismatch(
            "Expected a $((nsamples, K)) knn weights, got a $(size(weights)) knn "*
            "weights matrix"
        )
    )
    eltype(weights) <: AbstractFloat || throw(
        error(
            "The element type of knn weights matrix must be a subtype of `AbstractFloat`"
        )
    )
    return weights
end

####
# Uniform Kernel (Uniform)
####

"""
    Uniform()

assigns uniform weights to all k-nearest neighbors.

see also:
[`Inverse`](@ref), [`ISquared`](@ref), [`Zavreal`](@ref)
"""
struct Uniform <: KNNKernel end

get_weights(::Uniform, dists_matrix) = nothing

####
# Dudani
####

"""
    Dudani()

Assigns the closest neighbor a weight of `1`, the furthest neighbor weight `0` and the
others are scaled between by a linear mapping.

For more information see the paper by Geler et.al [Comparison of different weighting
schemes for the kNN classifier on time-series data](https://perun.pmf.uns.ac.rs/
radovanovic/publications/2016-kais-knn-weighting.pdf).

see also:
[`Macleod`](@ref), [`DualU`](@ref), [`DualD`](@ref)
"""
struct Dudani <: KNNKernel end

sort_idxs(::Dudani) = true

function _dudani_weights(dists_matrix, inds)
    i, j = inds.I
    dij = dists_matrix[i, j]
    dik = dists_matrix[i, end]
    di1 = dists_matrix[i, 1]
    w = (dik - dij)/ (dik - di1)
    return ifelse(isfinite(w), w, oneunit(w))
end

function get_weights(::Dudani, dists_matrix)
    return _dudani_weights.(Ref(dists_matrix), eachindex(dists_matrix))
end

####
# Inverse
####

"""
    Inverse()

Assigns each neighbor a weight equal to the inverse of the corresponsting distance of the
neighbor.

For more information see the paper by Geler et.al [Comparison of different weighting
schemes for the kNN classifier on time-series data](https://perun.pmf.uns.ac.rs/
radovanovic/publications/2016-kais-knn-weighting.pdf).

see also:
[`ISquared`](@ref)
"""
struct Inverse <: KNNKernel end

function _inverse_weights(dist)
    w = 1/dist
    return w
end

function replace_infinite_weights!(weights)
    inf_mask = (!isfinite).(weights)
    inf_row = @view(any(inf_mask, dims=2)[:])
    @inbounds(weights[inf_row, :] .= @view(inf_mask[inf_row, :]))
    return weights
end

function get_weights(::Inverse, dists_matrix)
   weights =  _inverse_weights.(dists_matrix)
   # Check and correct for non-finite rows
   replace_infinite_weights!(weights)
   return weights
end

####
# Inverse Squared (ISquared)
####

## Takes SKLearn's approah of removing non finite values,
## instead of adding a constant for numerical stability.
"""
    ISquared()

Assigns each neighbor a weight equal to the inverse of the corresponsting squared-distance
of the neighbor.

For more information see the paper by Geler et.al [Comparison of different weighting
schemes for the kNN classifier on time-series data](https://perun.pmf.uns.ac.rs/
radovanovic/publications/2016-kais-knn-weighting.pdf).

"""
struct ISquared <: KNNKernel end

function _inversesquared_weights(dist)
    w = 1/dist^2
    return w
end

function get_weights(::ISquared, dists_matrix)
    weights = _inversesquared_weights.(dists_matrix)
    # Check and correct for non-finite rows
    replace_infinite_weights!(weights)
    return weights
end

####
# Rank
####

"""
    Rank()

Assigns each neighbor a weight as a rank such that the closest neighbor get's a weight of
`1` and the Kth closest neighbor gets a weight of `K`.

For more information see the paper by Geler et.al [Comparison of different weighting
schemes for the kNN classifier on time-series data](https://perun.pmf.uns.ac.rs/
radovanovic/publications/2016-kais-knn-weighting.pdf).

see also:
[`ReciprocalRank`](@ref)
"""
struct Rank <: KNNKernel end

sort_idxs(::Rank) = true

function get_weights(kernel::Rank, dists_matrix)
    k = size(dists_matrix, 2)
    weights = similar(dists_matrix, Int)
    for i in axes(weights, 2)
        weights[:, i] .= k - i + 1
    end
    return weights
end

####
# Macleod
####

#implemented with `s` = `k`
"""
    Macleod(;a::Real= 0.0)

Assigns the closest neighbor a weight of `1`, the furthest neighbor weight `0` and the
others are scaled between by a linear mapping.

For more information see the paper by Geler et.al [Comparison of different weighting
schemes for the kNN classifier on time-series data](https://perun.pmf.uns.ac.rs/
radovanovic/publications/2016-kais-knn-weighting.pdf).

see also:
[`Dudani`](@ref), [`DualU`](@ref), [`DualD`](@ref)

"""
struct Macleod{T <: AbstractFloat} <: KNNKernel
    a::T
    Macleod(a::AbstractFloat) = new{typeof(a)}(a)
end

Macleod(a::Real) = Macleod(float(a))
Macleod(;a=1.0) = Macleod(a)
sort_idxs(::Macleod) = true

function _macleod_weights(dists_matrix, inds, s, a)
    i, j = inds.I
    dij = dists_matrix[i, j]
    dis = dists_matrix[i, s]
    di1 = dists_matrix[i, 1]
    w = ((dis - dij) + a*(dis - di1)) / (1 + a)*(dis - di1)
    return ifelse(isfinite(w), w, oneunit(w))
end

function get_weights(kernel::Macleod, dists_matrix::Mat{T}) where {T}
    k = size(dists_matrix, 2)
    s = k
    #s = ifelse(iszero(kernel.s), k, kernel.s)
    #s >= 0 || throw(ArgumentError("`s` can't be less than number of neighbors `K`"))
    return _macleod_weights.(Ref(dists_matrix), eachindex(dists_matrix), s, T(kernel.a))
end

####
# Zavreal
####

"""
    Zavreal(;s::Real = 0.0, a::Real=1.0)

Assigns each neighbor an exponential weight given by
    ``e^{ - α ⋅ d_i^{\beta}}``
where `α` and `β` are constants and `dᵢ` is the distance of the given neighbor.

For more information see the paper by Geler et.al [Comparison of different weighting
schemes for the kNN classifier on time-series data](https://perun.pmf.uns.ac.rs/
radovanovic/publications/2016-kais-knn-weighting.pdf).

"""
struct Zavreal{T <: AbstractFloat} <: KNNKernel
    a::T
    b::T
    Zavreal(a::K, b::K) where {K <: AbstractFloat} = new{K}(a, b)
end

Zavreal(;a=0.0, b= 1.0) = Zavreal(a, b)
Zavreal(a::Real, b::Real) = Zavreal(float(a), float(b))

function get_weights(kernel::Zavreal, dists_matrix::Mat{T}) where {T}
    return exp.(T(kernel.a) .* dists_matrix .^ T(kernel.b))
end

####
# ReciprocalRank
####
"""
    ReciprocalRank(;a::Real= 0.0)

Assigns each closest neighbor a weight which is equal to the reciprocal of it's rank.
i.e the closest neighbor get's a weight of `1` and the Kth closest weight get's a weight
of `1/K`

For more information see the paper by Geler et.al [Comparison of different weighting
schemes for the kNN classifier on time-series data](https://perun.pmf.uns.ac.rs/
radovanovic/publications/2016-kais-knn-weighting.pdf).

see also:
[`Rank`](@ref)

"""
struct ReciprocalRank <: KNNKernel end

sort_idxs(::ReciprocalRank) = true

function get_weights(kernel::ReciprocalRank, dists_matrix)
    rowvector = 1 ./ axes(dists_matrix, 2)
    nrows = size(dists_matrix, 1)
    return RankOneMatrix(rowvector, nrows)
end

####
# DualU
####
"""
    DualU()

Assigns the closest neighbor a weight of `1`, the furthest neighbor weight `0` and the
others are scaled between by a mapping.

For more information see the paper by Geler et.al [Comparison of different weighting
schemes for the kNN classifier on time-series data](https://perun.pmf.uns.ac.rs/
radovanovic/publications/2016-kais-knn-weighting.pdf).

see also:
[`DualD`](@ref)

"""
struct DualU <: KNNKernel end

sort_idxs(::DualU) = true

function _dualu_weights(dists_matrix, inds)
    i, j = inds.I
    dij = dists_matrix[i, j]
    dik = dists_matrix[i, end]
    di1 = dists_matrix[i, 1]
    w = (dik - dij)/ (dik - di1)
    return ifelse(isfinite(w), w, oneunit(w))
end

function get_weights(::DualU, dists_matrix)
    return _dualu_weights.(Ref(dists_matrix), eachindex(dists_matrix))
end

####
# DualD
####

"""
    DualD()

Assigns the closest neighbor a weight of `1`, the furthest neighbor weight `0` and the
others are scaled between by a mapping.

For more information see the paper by Geler et.al [Comparison of different weighting
schemes for the kNN classifier on time-series data](https://perun.pmf.uns.ac.rs/
radovanovic/publications/2016-kais-knn-weighting.pdf).

see also:
[`DualU`](@ref)

"""
struct DualD <: KNNKernel end

sort_idxs(::DualD) = true

function _duald_weights(dists_matrix, inds)
    i, j = inds.I
    dij = dists_matrix[i, j]
    dik = dists_matrix[i, end]
    di1 = dists_matrix[i, 1]
    w = ((dik - dij)/ (dik - di1)) * ((dik + di1) / (dik + dij))
    return ifelse(isfinite(w), w, oneunit(w))
end

# eachindex(dists_matrix) gives `CartesianIndices` because
# `Base.IndexStyle(dists_matrix) == IndexCartesian`
function get_weights(::DualD, dists_matrix)
    return _dualu_weights.(Ref(dists_matrix), eachindex(dists_matrix))
end

####
# Fibonacci
####
"""
    Fibonacci()

Assigns neighbors weights corresponding to fibonacci numbers starting from the furthest
neighbor. i.e the furthest neighbor a weight of `1`, the second furthest neighbor a
weight of `1` and the third furthest neighbor a weight of `2` and so on.

For more information see the paper by Geler et.al [Comparison of different weighting
schemes for the kNN classifier on time-series data](https://perun.pmf.uns.ac.rs/
radovanovic/publications/2016-kais-knn-weighting.pdf).

"""
struct Fibonacci <: KNNKernel end

sort_idxs(::Fibonacci) = true

function get_weights(::Fibonacci, dists_matrix)
    nrows, ncols = size(dists_matrix)
    rowvector = similar(Array{Int}, ncols)
    rowvector[end] = rowvector[end - 1] = 1
    for i in  StepRange(ncols - 2, -1, 1)
       rowvector[i] = rowvector[i+1] + rowvector[i+2]
    end
    return RankOneMatrix(rowvector, nrows)
end


text:
#comment adding these packages in and out as needed
#import Pkg;
#Pkg.add("LightGraphs")
#Pkg.add("SimpleWeightedGraphs")

using LightGraphs
using SimpleWeightedGraphs

# Finds the node with the lowest edge weight currently in key
function minKey(key, visited, numberOfNodes)
    # Creates global variables
    globalMin = Inf;
    globalIndex = 1;

    # Create a mutex
    m = Threads.Mutex();

    # Goes through each unvisited node and checks to see if its edge weight is small than the current local minimum weight
    Threads.@threads for i in 1:numberOfNodes
        localMin = globalMin[];
        localIndex = globalIndex[];
        if ((visited[i] == 0) && (key[i] < localMin))
            localMin = key[i];
            localIndex = i;
        end

        # Updates the global variables
        lock(m);
        if (localMin < globalMin[])
            globalMin = localMin;
            globalIndex = localIndex;
        end
        unlock(m);
    end

    return globalIndex[];
end

function parallelPrims(G, numberOfNodes)
    # Create adjacency matrix
    weights_matrix = weights(G);

    # Intialise arrays
    from = zeros(numberOfNodes); # tracks how to travel to each node
    key = zeros(numberOfNodes); # distance vector containing the current smallest weight to reach each node
    visited = zeros(numberOfNodes); # tracks if a node has been visited

    # Populate key array with Infinity
    for i in 1:numberOfNodes
        key[i] = Inf;
    end

    key[1] = 0;
    from[1] = -1; # -1 indicates that this is the starting point of the MST

    # Chooses which node to travel to based off what has the lowest value in key
    while 0 in visited
        nextNode = minKey(key, visited, numberOfNodes);
        visited[nextNode] = 1;

        # Updates the key array so that it contains the lowest weights of all possible paths
        Threads.@threads for possibleNode in 1:numberOfNodes
            # checks that there is an edge between the nodes, that the node hasn't been visited already and that the weight is smaller than the current edge weight in key
            if ((weights_matrix[nextNode, possibleNode] != 0) && (visited[possibleNode] == 0) && (weights_matrix[nextNode, possibleNode] < key[possibleNode]))
                from[possibleNode] = nextNode;
                key[possibleNode] = weights_matrix[nextNode, possibleNode];
            end
        end
    end

    # Creates an array containing the edges that make up the MST
    edges = [];
    for i in 2:numberOfNodes
        currentEdge = Edge(Int(from[i]), i);
        push!(edges, currentEdge);
    end
    return edges;
end


text:
using Libz
println(string(Libz.crc32(UInt8.(b"The quick brown fox jumps over the lazy dog")), base=16))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# %% Recurrences
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, DynamicalBilliards, PyPlot
w = 0.2
pc = Particle(-0.01, 0.2, sqrt(3)/2)
pr = Particle(0.0, 1.9, 0.0)

ps = [pc, pr]
cs = ["C0", "C3"]
bd = billiard_mushroom(1, w, 1)

zbox = ((420, 130), (460, 170))

fig, axs = subplots(1, 3)
εs = [0.1, 0.1]
for i in 1:2
    ax = axs[i==1 ? 1 : 3]
    p = ps[i]
    N = 500
    bmap, = boundarymap(p, bd, N)
    tr = standardize(Dataset(bmap))
    R = RecurrenceMatrix(tr, εs[i])
    x, y = coordinates(R)
    ax.scatter(x, y, s = 1, color = cs[i])
    ax.set_aspect("equal")
    ℓ = dl_entropy(R)
    r = rt_entropy(R)
    ax.set_xlim(0, N)
    ax.set_ylim(0, N)
    ax.set_xticks(0:250:500)
    ax.set_yticks(0:250:500)
    ax.grid()
    ax.set_title("\$H_\\ell = $(round(ℓ;digits= 1)),\\, H_r =$(round(r;digits=1))\$")
end
axs[3].set_yticklabels([])

# add zoomin
axis_zoomin!(axs[2], axs[1], zbox, zbox, "C1")
axs[2].axis("off")
axs[2].set_aspect("equal")
bmap, = boundarymap(pc, bd, 500)
tr = standardize(Dataset(bmap))
R = RecurrenceMatrix(tr, εs[1])
x, y = coordinates(R)
axs[2].scatter(x, y, s = 10, color = cs[1])
axs[2].set_xlim(zbox[1][1], zbox[2][1])
axs[2].set_ylim(zbox[1][2], zbox[2][2])

# Create indicative arrows
y0 = 137
yspan = 12
yc = (2y0 + yspan)/2
nice_arrow!(axs[2], 427, yc, 0, yspan; tex = "\$r=$(yspan)\$", color = "C2")

x0 = 445
y0 = 140
xspan = 10
yspan = 10
yc = (2y0 + yspan)/2
xc = (2x0 + xspan)/2
nice_arrow!(axs[2], xc, yc, xspan, yspan; tex = "\$\\ell=10\$", color = "C4")

add_identifiers!(fig, (axs[1], axs[3]))
fig.subplots_adjust(left = 0.08, bottom = 0.01, right = 0.95, top = 0.97, hspace = 0.1)
wsave(plotsdir("8","recurrence"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
using LuxurySparse

using Yao
using Yao.Intrinsics
using Random
using BenchmarkTools
using Statistics
using LinearAlgebra
using SparseArrays

vreg = randn(ComplexF64, 1<<16)
mreg = randn(ComplexF64, 1<<16, 2)

randmat(::Type{Val{:Dense}}, n) = randn(ComplexF64,n,n)
randmat(::Type{Val{:Diag}}, n) = Diagonal(randn(ComplexF64,n))
randmat(::Type{Val{:Perm}}, n) = pmrand(ComplexF64, n)
randmat(::Type{Val{:CSC}}, n) = SparseMatrixCSC(pmrand(ComplexF64, n))

vs = [vreg, mreg]

for m in [:Diag, :Dense, :Perm, :Diag]
    for v in [1,2]
        token = "$m * $(v ==1 ? "Vec" : "Mat")"
        for ng in [1,2]
            n=1<<ng
            un = randmat(Val{m}, n)
            sun = un |> staticize
            inds =  randperm(n) .+18
            sinds = SVector{n}(inds)

            if ng==1 && m == :Dense
                println("$token, u1rows($ng) => ", median(@benchmark u1rows!($(vs[v]), $(inds[1]), $(inds[2]), $(un[1]), $(un[3]), $(un[2]), $(un[4]))))
            end
            #println("$token, dynamic($ng) => ", median(@benchmark unrows!($(vs[v]), $inds, $un) seconds=1))
            res = @benchmark unrows!($(vs[v]), $sinds, $sun)
            println("$token, static($ng) => ", median(res))
            display(res)
        end
    end
end

using StaticArrays: SVector, SMatrix, SDiagonal

@inline function unrows!(state::Vector, inds::AbstractVector, A::Union{SSparseMatrixCSC, SparseMatrixCSC}, work::Vector)
    work .= 0
    @inbounds for col = 1:length(inds)
        xj = state[inds[col]]
        @inbounds @simd for j = A.colptr[col]:(A.colptr[col + 1] - 1)
            work[A.rowval[j]] += A.nzval[j]*xj
        end
    end
    state[inds] = work
    state
end

n=2
inds =  randperm(n) .+18
sinds = SVector{n}(inds)
un = Diagonal(randn(ComplexF64, n))
#un = randmat(Val{:CSC}, n)
sun = un |> staticize;

work = Vector{ComplexF64}(undef, n)
#@benchmark unrows!($vreg, $sinds, $sun, $work)
@benchmark unrows!($mreg, $sinds, $sun)


text:
# P is the type of polynomials in the rhs of the ODE system
"""
The main structure that represents input ODE system.

Stores information about states (`x_vars`), outputs (`y_vars`), inputs (`u_vars`), parameters (`parameters`) and the equations.

This structure is constructed via `@ODEmodel` macro.
"""
struct ODE{P}
    poly_ring::MPolyRing
    x_vars::Array{P, 1}
    y_vars::Array{P, 1}
    u_vars::Array{P, 1}
    parameters::Array{P, 1}
    x_equations::Dict{P, <: Union{P, Generic.Frac{P}}}
    y_equations::Dict{P, <: Union{P, Generic.Frac{P}}}
    function ODE{P}(
            x_eqs::Dict{P, <: Union{P, Generic.Frac{P}}},
            y_eqs::Dict{P, <: Union{P, Generic.Frac{P}}},
            inputs::Array{P, 1}
        ) where {P <: MPolyElem{<: FieldElem}}
        # Initialize ODE
        # x_eqs is a dictionary x_i => f_i(x, u, params)
        # y_eqs is a dictionary y_i => g_i(x, u, params)

        num, den = unpack_fraction(collect(values(x_eqs))[1])
        poly_ring = parent(num)
        if !all(isascii.(string.(gens(poly_ring))))
            nonascii_chars = filter(g->!isascii(g), string.(gens(poly_ring)))
            st = join(nonascii_chars, ", ")
            @warn "Non-ascii characters are not supported by Singular: " * st
        end
        x_vars = collect(keys(x_eqs))
        y_vars = collect(keys(y_eqs))
        u_vars = inputs
        parameters = filter(v -> (!(v in x_vars) && !(v in u_vars) && !(v in y_vars)), gens(poly_ring))
        new{P}(poly_ring, x_vars, y_vars, u_vars, parameters, x_eqs, y_eqs)
    end
end

#------------------------------------------------------------------------------

function add_outputs(ode::ODE{P}, extra_y::Dict{String, <: RingElem}) where P <: MPolyElem
    new_var_names = vcat(collect(map(var_to_str, gens(ode.poly_ring))), collect(keys(extra_y)))
    new_ring, new_vars = Nemo.PolynomialRing(base_ring(ode.poly_ring), new_var_names)

    new_x_eqs = Dict{P, Union{P, Generic.Frac{P}}}(parent_ring_change(x, new_ring) => parent_ring_change(f, new_ring) for (x, f) in ode.x_equations)
    new_y_eqs = Dict{P, Union{P, Generic.Frac{P}}}(parent_ring_change(y, new_ring) => parent_ring_change(g, new_ring) for (y, g) in ode.y_equations)
    extra_y_eqs = Dict{P, Union{P, Generic.Frac{P}}}(str_to_var(y, new_ring) => parent_ring_change(g, new_ring) for (y, g) in extra_y)
    merge!(new_y_eqs, extra_y_eqs)
    new_us = map(v -> switch_ring(v, new_ring), ode.u_vars)
    return ODE{P}(new_x_eqs, new_y_eqs, new_us)
end

#------------------------------------------------------------------------------

"""
    set_parameter_values(ode, param_values)

Input:
- `ode` - an ODE as above
- `param_values` - values for (possibly, some of) the parameters as dictionary `parameter` => `value`

Output:
- new ode with the parameters in param_values plugged with the given numbers
"""
function set_parameter_values(ode::ODE{P}, param_values::Dict{P, T}) where {T <: FieldElem, P <: MPolyElem{T}}
    new_vars = map(var_to_str, [v for v in gens(ode.poly_ring) if !(v in keys(param_values))])
    small_ring, small_vars = Nemo.PolynomialRing(base_ring(ode.poly_ring), new_vars)
    eval_dict = Dict(str_to_var(v, ode.poly_ring) => str_to_var(v, small_ring) for v in new_vars)
    merge!(eval_dict, Dict(p => small_ring(val) for (p, val) in param_values))

    return ODE{P}(
        Dict{P, Union{P, Generic.Frac{P}}}(eval_at_dict(v, eval_dict) => eval_at_dict(f, eval_dict) for (v, f) in ode.x_equations),
        Dict{P, Union{P, Generic.Frac{P}}}(eval_at_dict(v, eval_dict) => eval_at_dict(f, eval_dict) for (v, f) in ode.y_equations),
        [eval_at_dict(u, eval_dict) for u in ode.u_vars]
    )
end

#------------------------------------------------------------------------------

"""
    power_series_solution(ode, param_values, initial_conditions, input_values, prec)

Input:
- `ode` - an ode to solve
- `param_values` - parameter values, must be a dictionary mapping parameter to a value
- `initial_conditions` - initial conditions of `ode`, must be a dictionary mapping state variable to a value
- `input_values` - power series for the inpiuts presented as a dictionary variable => list of coefficients
- `prec` - the precision of solutions

Output:
- computes a power series solution with precision prec presented as a dictionary variable => corresponding coordiante of the solution
"""
function power_series_solution(
        ode::ODE{P},
        param_values::Dict{P, T},
        initial_conditions::Dict{P, T},
        input_values::Dict{P, Array{T, 1}},
        prec::Int
    ) where {T <: FieldElem, P <: MPolyElem{T}}
    new_varnames = map(var_to_str, vcat(ode.x_vars, ode.u_vars))
    append!(new_varnames, map(v -> var_to_str(v) * "_dot", ode.x_vars))

    new_ring, new_vars = Nemo.PolynomialRing(base_ring(ode.poly_ring), new_varnames)
    equations = Array{P, 1}()
    evaluation = Dict(k => new_ring(v) for (k, v) in param_values)
    for v in vcat(ode.x_vars, ode.u_vars)
        evaluation[v] = switch_ring(v, new_ring)
    end
    for (v, eq) in ode.x_equations
        num, den = map(p -> eval_at_dict(p, evaluation), unpack_fraction(eq))
        push!(equations, den * str_to_var(var_to_str(v) * "_dot", new_ring) - num)
    end
    new_inputs = Dict(switch_ring(k, new_ring) => v for (k, v) in input_values)
    new_ic = Dict(switch_ring(k, new_ring) => v for (k, v) in initial_conditions)
    result = ps_ode_solution(equations, new_ic, new_inputs, prec)

    # Evaluation outputs
    result = Dict(v => result[switch_ring(v, new_ring)] for v in vcat(ode.x_vars, ode.u_vars))
    ps_ring = parent(first(values(result)))
    for p in ode.parameters
        result[p] = ps_ring(param_values[p])
    end
    eval_outputs = []
    for (y, g) in ode.y_equations
        result[y] = eval_at_dict(g, result)
    end
    return result
end

#------------------------------------------------------------------------------

function power_series_solution(
        ode::ODE{P},
        param_values::Dict{P, Int},
        initial_conditions::Dict{P, Int},
        input_values::Dict{P, Array{Int, 1}},
        prec::Int
    ) where P <: MPolyElem{<: FieldElem}
    bring = base_ring(ode.poly_ring)
    return power_series_solution(
        ode,
        Dict(p => bring(v) for (p, v) in param_values),
        Dict(x => bring(v) for (x, v) in initial_conditions),
        Dict(u => map(v -> bring(v), vv) for (u, vv) in input_values),
        prec
    )
end
#------------------------------------------------------------------------------

"""
    _reduce_mod_p(f, p)

Reduces a polynomial/rational function over Q modulo p
"""
function _reduce_mod_p(poly::fmpq_mpoly, p::Int)
    den = denominator(poly)
    num = change_base_ring(Nemo.ZZ, den * poly)
    if Nemo.GF(p)(den) == 0
        throw(Base.ArgumentError("Prime $p divides the denominator of $poly"))
    end
    return change_base_ring(Nemo.GF(p), num) * (1 // Nemo.GF(p)(den))
end

function _reduce_mod_p(rat::Generic.Frac{fmpq_mpoly}, p::Int)
    num, den = map(poly -> _reduce_mod_p(poly, p), [numerator(rat), denominator(rat)])
    if den == 0
        throw(Base.ArgumentError("Prime $p divides the denominator of $rat"))
    end
    return num // den
end

#--------------------------------------

"""
    reduce_ode_mod_p(ode, p)

Input: ode is an ODE over QQ, p is a prime number
Output: the reduction mod p, throws an exception if p divides one of the denominators
"""
function reduce_ode_mod_p(ode::ODE{<: MPolyElem{Nemo.fmpq}}, p::Int)
    new_ring, new_vars = Nemo.PolynomialRing(Nemo.GF(p), map(var_to_str, gens(ode.poly_ring)))
    new_type = typeof(new_vars[1])
    new_inputs = map(u -> switch_ring(u, new_ring), ode.u_vars)
    new_x_eqs = Dict{new_type, Union{new_type, Generic.Frac{new_type}}}()
    new_y_eqs = Dict{new_type, Union{new_type, Generic.Frac{new_type}}}()
    for (old, new) in Dict(ode.x_equations => new_x_eqs, ode.y_equations => new_y_eqs)
        for (v, f) in old
            new_v = switch_ring(v, new_ring)
            new[new_v] = _reduce_mod_p(f, p)
        end
    end
    return ODE{new_type}(new_x_eqs, new_y_eqs, new_inputs)
end

#------------------------------------------------------------------------------

function macrohelper_extract_vars(equations::Array{Expr, 1})
    funcs, x_vars, all_symb = Set(), Set(), Set()
    aux_symb = Set([:(+), :(-), :(=), :(*), :(^), :t, :(/), :(//)])
    for eq in equations
        MacroTools.postwalk(
            x -> begin
                if @capture(x, f_'(t))
                    push!(x_vars, f)
                    push!(all_symb, f)
                elseif @capture(x, f_(t))
                    push!(funcs, f)
                elseif (x isa Symbol) && !(x in aux_symb)
                    push!(all_symb, x)
                end
                return x
            end,
            eq
        )
    end
    io_vars = setdiff(funcs, x_vars)
    all_symb = collect(all_symb)
    return collect(x_vars), collect(io_vars), collect(all_symb)
end

function macrohelper_extract_vars(equations::Array{Symbol, 1})
    return macrohelper_extract_vars(map(Expr, equations))
end

#------------------------------------------------------------------------------

function macrohelper_clean(ex::Expr)
    ex = MacroTools.postwalk(x -> @capture(x, f_'(t)) ? f : x, ex)
    ex = MacroTools.postwalk(x -> @capture(x, f_(t)) ? f : x, ex)
    ex = MacroTools.postwalk(x -> x == :(/) ? :(//) : x, ex)
    return ex
end

#------------------------------------------------------------------------------

"""
Macro for creating an ODE from a list of equations.
Also injects all variables into the global scope.

This macro accepts a sybolically written ODE system and generates an `ODE` structure instance:
```julia
ode = @ODEmodel(
    x1'(t) = -k1 * x1(t),
    y1(t) = x1(t)
)
```
"""
macro ODEmodel(ex::Expr...)
    equations = [ex...]
    x_vars, io_vars, all_symb = macrohelper_extract_vars(equations)

    # creating the polynomial ring
    vars_list = :([$(all_symb...)])
    R = gensym()
    vars_aux = gensym()
    exp_ring = :(($R, $vars_aux) = StructuralIdentifiability.Nemo.PolynomialRing(
        StructuralIdentifiability.Nemo.QQ,
        map(string, $all_symb)
    ))
    assignments = [:($(all_symb[i]) = $vars_aux[$i]) for i in 1:length(all_symb)]

    # preparing equations
    equations = map(macrohelper_clean, equations)
    x_dict = gensym()
    y_dict = gensym()
    y_vars = Set()
    x_dict_create_expr = :($x_dict = Dict{
        StructuralIdentifiability.Nemo.fmpq_mpoly,
        Union{StructuralIdentifiability.Nemo.fmpq_mpoly, StructuralIdentifiability.AbstractAlgebra.Generic.Frac{StructuralIdentifiability.Nemo.fmpq_mpoly}}
    }())
    y_dict_create_expr = :($y_dict = Dict{
        StructuralIdentifiability.Nemo.fmpq_mpoly,
        Union{StructuralIdentifiability.Nemo.fmpq_mpoly, StructuralIdentifiability.AbstractAlgebra.Generic.Frac{StructuralIdentifiability.Nemo.fmpq_mpoly}}
    }())
    eqs_expr = []
    for eq in equations
        if eq.head != :(=)
            throw("Problem with parsing at $eq")
        end
        lhs, rhs = eq.args[1:2]
        loc_all_symb = macrohelper_extract_vars([rhs])[3]
        to_insert = undef
        if lhs in x_vars
            to_insert = x_dict
        elseif lhs in io_vars
            to_insert = y_dict
            push!(y_vars, lhs)
        else
            throw("Unknown left-hand side $lhs")
        end
        if isempty(loc_all_symb)
            push!(eqs_expr, :($to_insert[$lhs] = $R($rhs)))
        else
            push!(eqs_expr, :($to_insert[$lhs] = ($rhs)))
        end
    end

    u_vars = setdiff(io_vars, y_vars)
    params = setdiff(all_symb, union(x_vars, y_vars, u_vars))
    @info "Summary of the model:"
    @info "State variables: " * join(map(string, collect(x_vars)), ", ")
    @info "Parameters: " * join(map(string, collect(params)), ", ")
    @info "Inputs: " * join(map(string, collect(u_vars)), ", ")
    @info "Outputs: " * join(map(string, collect(y_vars)), ", ")

    # creating the ode object
    ode_expr = :(ODE{StructuralIdentifiability.Nemo.fmpq_mpoly}($x_dict, $y_dict, Array{StructuralIdentifiability.Nemo.fmpq_mpoly}([$(u_vars...)])))

    result = Expr(
        :block,
        exp_ring, assignments...,
        x_dict_create_expr, y_dict_create_expr, eqs_expr...,
        ode_expr
    )
    return esc(result)
end

#------------------------------------------------------------------------------

function Base.show(io::IO, ode::ODE)
    varstr = Dict(x => var_to_str(x) * "(t)" for x in vcat(ode.x_vars, ode.u_vars, ode.y_vars))
    merge!(varstr, Dict(p => var_to_str(p) for p in ode.parameters))
    R_print, vars_print = Nemo.PolynomialRing(base_ring(ode.poly_ring), [varstr[v] for v in gens(ode.poly_ring)])
    for (x, eq) in ode.x_equations
        print(io, var_to_str(x) * "'(t) = ")
        print(io, evaluate(eq, vars_print))
        print(io, "\n")
    end
    for (y, eq) in ode.y_equations
        print(io, var_to_str(y) * "(t) = ")
        print(io, evaluate(eq, vars_print))
        print(io, "\n")
    end
end

#------------------------------------------------------------------------------
"""
    function PreprocessODE(de::ModelingToolkit.ODESystem, inputs)

Input:
- `diff_eqs` - array of ModelingToolkit differential equations
- `out_eqs` - array of output equations
- `states` - array of state variables
- `outputs` - array of output function names
- `inputs` - array of input function names
- `parameters` - array of parameter names

Output:
- `ODE` object containing required data for identifiability assessment
"""
function PreprocessODE(de::ModelingToolkit.ODESystem, measured_quantities::Array{ModelingToolkit.Equation})
    @info "Preproccessing `ModelingToolkit.ODESystem` object"
    diff_eqs = filter(eq->!(ModelingToolkit.isoutput(eq.lhs)), ModelingToolkit.equations(de))
    y_functions = [each.lhs for each in measured_quantities]
    inputs = filter(v->ModelingToolkit.isinput(v), ModelingToolkit.states(de))
    state_vars = filter(s->!(ModelingToolkit.isinput(s) || ModelingToolkit.isoutput(s)), ModelingToolkit.states(de))
    params = ModelingToolkit.parameters(de)
    t = ModelingToolkit.arguments(measured_quantities[1].lhs)[1]
    params_from_measured_quantities = ModelingToolkit.parameters(ModelingToolkit.ODESystem(measured_quantities, t, name=:DataSeries))
    params = union(params, params_from_measured_quantities)

    input_symbols = vcat(state_vars, y_functions, inputs, params)
    generators = string.(input_symbols)
    generators = map(g->replace(g, "(t)"=>""), generators)
    R, gens_ = Nemo.PolynomialRing(Nemo.QQ, generators)
    state_eqn_dict = Dict{StructuralIdentifiability.Nemo.fmpq_mpoly,Union{StructuralIdentifiability.Nemo.fmpq_mpoly,StructuralIdentifiability.Nemo.Generic.Frac{StructuralIdentifiability.Nemo.fmpq_mpoly}}}()
    out_eqn_dict = Dict{StructuralIdentifiability.Nemo.fmpq_mpoly,Union{StructuralIdentifiability.Nemo.fmpq_mpoly,StructuralIdentifiability.Nemo.Generic.Frac{StructuralIdentifiability.Nemo.fmpq_mpoly}}}()

    for i in 1:length(diff_eqs)
        state_eqn_dict[substitute(state_vars[i], input_symbols.=>gens_)] = eval_at_nemo(diff_eqs[i].rhs, Dict(input_symbols.=>gens_))
    end
    for i in 1:length(measured_quantities)
        out_eqn_dict[substitute(y_functions[i], input_symbols.=> gens_)] = eval_at_nemo(measured_quantities[i].rhs, Dict(input_symbols.=>gens_))
    end

    inputs_ = [substitute(each,  input_symbols .=> gens_) for each in inputs]
    if isequal(length(inputs_), 0)
        inputs_ = Vector{StructuralIdentifiability.Nemo.fmpq_mpoly}()
    end
    return (ODE{StructuralIdentifiability.Nemo.fmpq_mpoly}(state_eqn_dict, out_eqn_dict, inputs_), input_symbols, gens_)
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -42.99754755332367 <= q <= 40.38858433100043)

                     add_NL_constraint(m, :(swish(-0.8925450194703983 + 0.6110890427800215*swish(-0.16626426792027038 + 0.44845618413236066*swish(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + 0.6799320945085388*swish(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + -0.7514544713734312*swish(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + 0.795720190818876*swish(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + 0.1740395626916258*swish(0.8063063231762757 + 0.4905142959916642*swish(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + 0.7092497630992938*swish(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + 0.8897979792487347*swish(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + -0.11224890387090136*swish(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + -0.8017072131372296*swish(0.9210510999614945 + -0.4935968211971775*swish(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + -0.33547745483400737*swish(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + 0.1703969445706579*swish(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + -0.2742027660464488*swish(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + -0.7966067754659663*swish(-0.7328670881809516 + 0.6102429954796835*swish(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + -0.6353392761755483*swish(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + -0.7830650903096812*swish(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + 0.1627621376749886*swish(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5])))) + swish(-0.7008872510227278 + -0.258173144464545*swish(-0.16626426792027038 + 0.44845618413236066*swish(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + 0.6799320945085388*swish(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + -0.7514544713734312*swish(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + 0.795720190818876*swish(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + 0.49550787111793193*swish(0.8063063231762757 + 0.4905142959916642*swish(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + 0.7092497630992938*swish(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + 0.8897979792487347*swish(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + -0.11224890387090136*swish(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + -0.8462943069272191*swish(0.9210510999614945 + -0.4935968211971775*swish(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + -0.33547745483400737*swish(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + 0.1703969445706579*swish(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + -0.2742027660464488*swish(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + -0.05869131996649157*swish(-0.7328670881809516 + 0.6102429954796835*swish(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + -0.6353392761755483*swish(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + -0.7830650903096812*swish(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + 0.1627621376749886*swish(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5])))) + swish(-0.2643911691309153 + -0.31660511055293883*swish(-0.16626426792027038 + 0.44845618413236066*swish(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + 0.6799320945085388*swish(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + -0.7514544713734312*swish(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + 0.795720190818876*swish(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + -0.09136517532701127*swish(0.8063063231762757 + 0.4905142959916642*swish(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + 0.7092497630992938*swish(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + 0.8897979792487347*swish(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + -0.11224890387090136*swish(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + -0.36400896662330684*swish(0.9210510999614945 + -0.4935968211971775*swish(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + -0.33547745483400737*swish(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + 0.1703969445706579*swish(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + -0.2742027660464488*swish(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + 0.4876638787119827*swish(-0.7328670881809516 + 0.6102429954796835*swish(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + -0.6353392761755483*swish(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + -0.7830650903096812*swish(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + 0.1627621376749886*swish(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5])))) + swish(0.6679332747384503 + -0.6061375406574494*swish(-0.16626426792027038 + 0.44845618413236066*swish(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + 0.6799320945085388*swish(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + -0.7514544713734312*swish(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + 0.795720190818876*swish(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + -0.7532640510772808*swish(0.8063063231762757 + 0.4905142959916642*swish(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + 0.7092497630992938*swish(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + 0.8897979792487347*swish(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + -0.11224890387090136*swish(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + 0.9959397069462161*swish(0.9210510999614945 + -0.4935968211971775*swish(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + -0.33547745483400737*swish(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + 0.1703969445706579*swish(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + -0.2742027660464488*swish(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))) + -0.8544693505755525*swish(-0.7328670881809516 + 0.6102429954796835*swish(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])) + -0.6353392761755483*swish(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])) + -0.7830650903096812*swish(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])) + 0.1627621376749886*swish(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions