text:
# Test script to test solver with JuMP on a closest correlation matrix problem
using COSMO, JuMP, LinearAlgebra, SparseArrays, Test, Random
rng = Random.MersenneTwister(12345);

# Original problem has the following format:
# min_X   1/2 ||X-C||^2
# s.t.    Xii = 1
#         X ⪴ 0

# create a random test matrix C
n = 8
C = -1 .+ rand(rng, n, n) .* 2;
c = vec(C);

# define problem in JuMP
q = -vec(C);
r = 0.5 * vec(C)' * vec(C);
m = JuMP.Model(with_optimizer(COSMO.Optimizer, verbose=true, eps_abs = 1e-4));
@variable(m, X[1:n, 1:n], PSD);
x = vec(X);
@objective(m, Min, 0.5 * x' * x  + q' * x + r)
for i = 1:n
  @constraint(m, X[i, i] == 1.)
end

# solve and get results
status = JuMP.optimize!(m)
obj_val = JuMP.objective_value(m)
X_sol = JuMP.value.(X)

known_opt_val = 12.5406
known_solution =  [
  1.0         0.732562   -0.319491   -0.359985   -0.287543   -0.15578     0.0264044  -0.271438;
  0.732562    1.0         0.0913246  -0.0386357   0.299199   -0.122733    0.126612   -0.187489;
 -0.319491    0.0913246   1.0        -0.0863377   0.432948    0.461783   -0.248641   -0.395299;
 -0.359985   -0.0386357  -0.0863377   1.0         0.503379    0.250601    0.141151    0.286088;
 -0.287543    0.299199    0.432948    0.503379    1.0        -0.0875199   0.137518    0.0262425;
 -0.15578    -0.122733    0.461783    0.250601   -0.0875199   1.0        -0.731556    0.0841783;
  0.0264044   0.126612   -0.248641    0.141151    0.137518   -0.731556    1.0        -0.436274;
 -0.271438   -0.187489   -0.395299    0.286088    0.0262425   0.0841783  -0.436274    1.0  ];

@testset "Closest correlation matrix example" begin
  @test isapprox(obj_val, known_opt_val , atol=1e-3)
  @test norm(X_sol - known_solution, Inf) < 1e-3
end
nothing


text:
using Pkg
Pkg.activate(".")

import DotEnv
DotEnv.config()

using MOSInterface

# Interface
interface = Interface()

# Delete model
delete_model_with_name(interface, "Portfolio Model")

# New model
model = new_model(interface, "./examples/cvxpy/portfolio/portfolio_model.py")

# Get model
model = get_model_with_name(interface, "Portfolio Model")

# Set inputs
set_interface_object(model, "L", 2)
set_interface_object(model, "gamma", 0.1)
set_interface_file(model, "stockdata", "examples/cvxpy/portfolio/stockdata.csv")

@assert(get_name(model) == "Portfolio Model")
@assert(get_system(model) == "cvxpy")
@assert(get_status(model) == "created")

# Solve
MOSInterface.run(model)

@assert(get_status(model) == "success")

ext: jl

max_stars_repo_name: Fuinn/mos-examples

max_issues_repo_name: Fuinn/mos-examples

max_forks_repo_name: Fuinn/mos-examples


text:
using Avro.Schemas

using Test

const PRIMITIVE_EXAMPLES = [
    ["\"$type_name\"" for type_name in PRIMITIVE_TYPES];
    ["{\"type\": \"$type_name\"}" for type_name in PRIMITIVE_TYPES]
    ]

const RECORD_EXAMPLES = [
    """
    {
        "type": "record",
        "name": "LongList",
        "aliases": ["LinkedLongs"],
        "fields" : [
            {
                "name": "value",
                "type": "long"
            },
            {
                "name": "next",
                "type": ["null", "LongList"]
            }
        ]
    }
    """
    ]

ENUM_EXAMPLES = [
    """
    {
        "type": "enum",
        "name": "Suit",
        "symbols": ["SPADES", "HEARTS", "DIAMONDS", "CLUBS"]
    }
    """
    ]

ARRAY_EXAMPLES = [
    """
    {
        "type": "array",
        "items": "string"
    }
    """
    ]

MAP_EXAMPLES = [
    """
    {
        "type": "map",
        "values": "long"
    }
    """
    ]

UNION_EXAMPLES = [
    "[\"null\", \"string\"]"
    ]

FIXED_EXAMPLES = [
    """
    {
        "type": "fixed",
        "name": "md5",
        "size": 256
    }
    """
    ]

@testset "Parsing schemas" begin
    @testset "Primitive" for example in PRIMITIVE_EXAMPLES
        Schemas.parse(example)
    end

    @testset "Record" for example in RECORD_EXAMPLES
        Schemas.parse(example)
    end

    @testset "Enum" for example in ENUM_EXAMPLES
        Schemas.parse(example)
    end

    @testset "Array" for example in ARRAY_EXAMPLES
        Schemas.parse(example)
    end

    @testset "Map" for example in MAP_EXAMPLES
        Schemas.parse(example)
    end

    @testset "Union" for example in UNION_EXAMPLES
        Schemas.parse(example)
    end

    @testset "Fixed" for example in FIXED_EXAMPLES
        Schemas.parse(example)
    end
end

@testset "Schema equality" begin
    @testset "Primitive" begin
        @test Schemas.NULL == NullSchema()
        @test Schemas.NULL == Schemas.PrimitiveSchema("null")

        @test Schemas.BOOLEAN == BooleanSchema()
        @test Schemas.BOOLEAN == Schemas.PrimitiveSchema("boolean")

        @test Schemas.INT == IntSchema()
        @test Schemas.INT == Schemas.PrimitiveSchema("int")

        @test Schemas.LONG == LongSchema()
        @test Schemas.LONG == Schemas.PrimitiveSchema("long")

        @test Schemas.FLOAT == FloatSchema()
        @test Schemas.FLOAT == Schemas.PrimitiveSchema("float")

        @test Schemas.DOUBLE == DoubleSchema()
        @test Schemas.DOUBLE == Schemas.PrimitiveSchema("double")

        @test Schemas.BYTES == BytesSchema()
        @test Schemas.BYTES == Schemas.PrimitiveSchema("bytes")

        @test Schemas.STRING == StringSchema()
        @test Schemas.STRING == Schemas.PrimitiveSchema("string")
    end
end


text:
function SeisSmooth2(A;L1=7,L2=7,Nrepeat=1)
    a1 = isodd(L1) ? 0 : 1
    a2 = isodd(L2) ? 0 : 1
    kernel = ones(eltype(A),L1,L2)/(L1*L2);
    for irepeat = 1 : Nrepeat
	A = conv2(A,kernel)
	A = A[floor(Int,L1/2) + 1:end-floor(Int,L1/2) + a1,floor(Int,L2/2) + 1:end-floor(Int,L2/2) + a2]
	end
    return A
end

function SeisSmooth2(A,h::Array{Header,1};L1=7,L2=7)
    A = SeisSmooth2(A,L1=L1,L2=L2,Nrepeat=Nrepeat)
    return A,h
end

function SeisSmooth2(in::String,out::String;key=["iaz";"iang"],L1=7,L2=7,Nrepeat=1)
    d, h, ext = SeisRead(in)
    nang = size(d,3)
    for iang = 1 : nang
	A = d[:,1,iang,1,:]
	A = SeisSmooth2(A,L1=L1,L2=L2,Nrepeat=Nrepeat)
	d[:,1,iang,1,:] = A
    end
    SeisWrite(out,d,h,ext)

    #@compat parameters = Dict(:L1=>L1,:L2=>L2)
    #SeisProcess(in,out,[SeisSmooth2],[parameters];group="gather",key=key)
end


text:
using Base.UTF8proc: category_abbrev

const HALF = ' '
const CAP = '\ue302'
const SPACE_PRIORITY = Dict('L' => 7, 'M' => 7, 'N' => 5, 'S' => 3,
                            'P' => 1, 'Z' => -1, 'C' => -3)

function space_priority(char)
    return SPACE_PRIORITY[category_abbrev(char)[1]]
end

"Simple reversible tokenizer"
function tokenize(s::AbstractString, decap::Bool=false, split_punctuation=True)
    toks = [Vector{Char}()]
    current_cat = 0 # or -2?
    for c in s
        cat = space_priority(c)
        if c == ' '
            push!(toks[end], HALF)
            push!(toks, [HALF])
            current_cat = 0
            continue
        elseif current_cat == 0
            push!(toks[end], c)
        elseif cat == current_cat && (cat > 2 || !split_punctuation)
            push!(toks[end], c) # HALF + c
        elseif cat <= 0 && current_cat <= 0
            push!(toks, [c])
        elseif cat < current_cat
            push!(toks[end], HALF)
            push!(toks, [c])
        else
            push!(toks, [HALF, c])
        end
        current_cat = cat
    end
    if toks[1] == Vector{Char}()
        toks = toks[2:end]
    end
    if current_cat > 0
        push!(toks[end], HALF)
    end
    toks = map(arr -> convert(String, arr), toks)
    if decap
        return decapitalize.(toks)
    end
    return toks
end

function decapitalize(tok)
    pre, tok = tok[1] == HALF ? (HALF, tok[2:end]) : ("", tok)
    if length(tok) < 1 || tok[1] == lowercase(tok[1])
        return string(pre, tok)
    end
    if (tok[1] == uppercase(tok[1]) &&
        (length(tok) == 1 || tok[2] != uppercase(tok[2])))
        return string(CAP, pre, lcfirst(tok))
    end
    return string(pre, tok)
end

function detokenize(l)
    text = replace(join(l, ""), "$CAP$HALF", "$HALF$CAP")
    text = replace(text, Regex("$HALF+"), s->" " ^ (length(s) ÷ 2))
    return replace(text, Regex("\\Q$CAP\\E.", "s"), s->uppercase(s[end]))
end


text:
"""
# module laminart

- Julia version: 1.4
- Author: niallcullinane
- Date: 2020-06-07

# Examples

```jldoctest
julia>
```
"""
module Laminartv1
    include("./LamKernels.jl")

using NNlib, ImageFiltering, Images, OffsetArrays
# , MEngProject.LamKernels

function f!(du, u, p, t)
    @inbounds begin
        x = @view u[:, :, 1:p.K]
        y = @view u[:, :, p.K+1:2*p.K]
        m = @view u[:, :, 2*p.K+1:3*p.K]
        z = @view u[:, :, 3*p.K+1:4*p.K]
        s = @view u[:, :, 4*p.K+1:5*p.K]

#    C = @view uu[:, :, 5*p.K+1:6*p.K]
#    H_z = @view uu[:, :, 6*p.K+1:7*p.K]

        v_p = @view u[:, :, 5*p.K+1]
        v_m = @view u[:, :, 5*p.K+2]
#    x_lgn = @view uu[:, :, 7*p.K+3]

        dx = @view du[:, :, 1:p.K]
        dy = @view du[:, :, p.K+1:2*p.K]
        dm = @view du[:, :, 2*p.K+1:3*p.K]
        dz = @view du[:, :, 3*p.K+1:4*p.K]
        ds = @view du[:, :, 4*p.K+1:5*p.K]

        dv_p = @view du[:, :, 5*p.K+1]
        dv_m = @view du[:, :, 5*p.K+2]

        x_lgn = fun_x_lgn(x, p)
        C = fun_v_C(v_p, v_m, p)
        H_z = fun_H_z(z, p)

        temp_dv_p = fun_dv(v_p, p.r, x_lgn, p)
        temp_dv_m = fun_dv(v_m, .-p.r, x_lgn, p)
        temp_dx = fun_dx_V1(x, C, z, p.x_V2, p)
        temp_dy = fun_dy(y, C, x, m, p)
        temp_dm = fun_dm(m, x, p)
        temp_dz = fun_dz(z, y, H_z, s, p)
        temp_ds = fun_ds(s, z, H_z, p)

        @. dv_p = temp_dv_p
        @. dv_m = temp_dv_m
        @. dx = temp_dx
        @. dy = temp_dy
        @. dm = temp_dm
        @. dz = temp_dz
        @. ds = temp_ds
   end
   nothing
end

function kernels(img::AbstractArray, p::NamedTuple)
    C_Q_temp = reshape(Array{eltype(img)}(undef, p.C_AB_l, p.C_AB_l * p.K), p.C_AB_l, p.C_AB_l, p.K)
    C_P_temp = copy(C_Q_temp)
    H_temp = reshape(Array{eltype(img)}(undef, p.H_l, p.H_l * p.K), p.H_l, p.H_l, p.K)
    T_temp = reshape(Array{eltype(img)}(undef, 1, 1 * p.K), 1, 1, p.K)     #ijk,  1x1xk,   ijk
    W_temp = reshape(Array{eltype(img)}(undef, 19, 19 * p.K * p.K), 19, 19, p.K, p.K)     #ijk,  1x1xk,   ijk
    for k ∈ 1:p.K
        θ = π*(k-1)/p.K
        C_Q_temp[:,:,k] = reflect(centered(LamKernels.kern_d(p.σ_2, θ)))           #ij ijk ijk
        C_P_temp[:,:,k] = reflect(centered(LamKernels.kern_b(p.σ_2, θ)))               #ij ijk ijk
        H_temp[:,:,k] = reflect(p.H_fact .* LamKernels.gaussian_rot(p.H_σ_x, p.H_σ_y, θ, p.H_l))  #ijk, ij for each k; ijk
        T_temp[:,:,k] = reshape([p.T_fact[k]],1,1)
#todo: generalise T and W for higher K
#         T_temp[:,:,k] = KernelFactors.gaussian(p.T_σ, p.K)
#         for l ∈ 1:p.K
#             W_temp[:,:,l,k] =
#         end
    end
    W_temp[:,:,1,1] = reflect(5 .* LamKernels.gaussian_rot(3,0.8,0,19) .+ LamKernels.gaussian_rot(0.4,1,0,19))
    W_temp[:,:,2,2] = reflect(5 .* LamKernels.gaussian_rot(3,0.8,0,19) .+ LamKernels.gaussian_rot(0.4,1,π/2,19))
    W_temp[:,:,1,2] = reflect(relu.(0.2 .- LamKernels.gaussian_rot(2,0.6,0,19) .- LamKernels.gaussian_rot(0.3,1.2,0,19)))
    W_temp[:,:,2,1] = reflect(relu.(0.2 .- LamKernels.gaussian_rot(2,0.6,0,19) .- LamKernels.gaussian_rot(0.3,1.2,π/2,19)))

# todo fix W kernel
#  W_temp[:,:,1,1] = reflect(LamKernels.gaussian_rot(3,0.8,0,19))
#     W_temp[:,:,2,2] = reflect(LamKernels.gaussian_rot(3,0.8,0,19))
#     W_temp[:,:,1,2] = reflect(LamKernels.gaussian_rot(3,0.8,0,19))
#     W_temp[:,:,2,1] = reflect(LamKernels.gaussian_rot(3,0.8,0,19))

# todo: fix range of W H
#     W_range = -(p.W_size-1)/2:(p.W_size-1)/2
#     H_range = -(p.H_size-1)/2:(p.H_size-1)/2
    W_range = -9:9
    H_range = -9:9

    temp_out = (
    k_gauss_1 = reflect(Kernel.gaussian(p.σ_1)),
    k_gauss_2 = reflect(Kernel.gaussian(p.σ_2)),
    k_C_d = C_Q_temp,
    k_C_b = C_P_temp,
    k_W_p = OffsetArray(W_temp, W_range, W_range, 1:p.K, 1:p.K),
    k_W_m = OffsetArray(W_temp, W_range, W_range, 1:p.K, 1:p.K),
    k_H = OffsetArray(H_temp, H_range, H_range, 1:p.K),
    k_T_p = (T_temp),
    k_T_m = (p.T_p_m .* T_temp),
    k_T_p_v2 = (p.T_v2_fact .* T_temp),
    k_T_m_v2 = (p.T_v2_fact .* p.T_p_m .* T_temp),
    dim_i = size(img)[1],
    dim_j = size(img)[2],
    x_V2 = reshape(zeros(typeof(img[1,1]), size(img)[1], size(img)[2] * p.K), size(img)[1], size(img)[2], p.K))
    return merge(p, temp_out)
end

function varables(I::AbstractArray, p::NamedTuple)
   v_p = zeros(typeof(I[1,1]), size(I)[1], size(I)[2])
   v_m = copy(v_p)
   x_lgn = copy(v_p)
   x = reshape(zeros(typeof(I[1,1]), size(I)[1], size(I)[2] * p.K), size(I)[1], size(I)[2], p.K)
   y = copy(x)
   m = copy(x)
   z = copy(x)
   s = copy(x)
   C = copy(x)
   H_z = copy(x)
   x_V2 = copy(x)
   y_V2 = copy(x)
   m_V2 = copy(x)
   z_V2 = copy(x)
   s_V2 = copy(x)
   H_z_V2 = copy(x)
   return [v_p, v_m, x_lgn, x, y, m, z, s, C, H_z, x_V2, y_V2, m_V2, z_V2, s_V2, H_z_V2]
end

function variables_sep(I::AbstractArray, p::NamedTuple)
   v_p = zeros(typeof(I[1,1]), size(I)[1], size(I)[2])
   v_m = copy(v_p)
   x_lgn = copy(v_p)
   x = reshape(zeros(typeof(I[1,1]), size(I)[1], size(I)[2] * p.K), size(I)[1], size(I)[2], p.K)
   y = copy(x)
   m = copy(x)
   z = copy(x)
   s = copy(x)
   C = copy(x)
   H_z = copy(x)
   x_V2 = copy(x)
   y_V2 = copy(x)
   m_V2 = copy(x)
   z_V2 = copy(x)
   s_V2 = copy(x)
   H_z_V2 = copy(x)
   return v_p, v_m, x_lgn, x, y, m, z, s, C, H_z, x_V2, y_V2, m_V2, z_V2, s_V2, H_z_V2
end

function variables_dict(I::AbstractArray, p::NamedTuple)
   v_p = zeros(typeof(I[1,1]), size(I)[1], size(I)[2])
   v_m = copy(v_p)
   x_lgn = copy(v_p)
   x = reshape(zeros(typeof(I[1,1]), size(I)[1], size(I)[2] * p.K), size(I)[1], size(I)[2], p.K)
   y = copy(x)
   m = copy(x)
   z = copy(x)
   s = copy(x)
   C = copy(x)
   H_z = copy(x)
   x_V2 = copy(x)
   y_V2 = copy(x)
   m_V2 = copy(x)
   z_V2 = copy(x)
   s_V2 = copy(x)
   H_z_V2 = copy(x)
   return v_p, v_m, x_lgn, x, y, m, z, s, C, H_z, x_V2, y_V2, m_V2, z_V2, s_V2, H_z_V2
end

# function varables(I::AbstractArray, p::NamedTuple)
#    ij = zeros(typeof(I[1,1]), size(I)[1], size(I)[2])
#    ijk = reshape(zeros(typeof(I[1,1]), size(I)[1], size(I)[2] * p.K), size(I)[1], size(I)[2], p.K)
#    temp_out = [v_p = copy(ij),
#    v_m = copy(ij),
#    x_lgn = copy(ij),
#    x = copy(ijk),
#    y = copy(ijk),
#    m = copy(ijk),
#    z = copy(ijk),
#    s = copy(ijk),
#    C = copy(ijk),
#    H_z = copy(ijk),
#    x_V2 = copy(ijk),
#    y_V2 = copy(ijk),
#    m_V2 = copy(ijk),
#    z_V2 = copy(ijk),
#    s_V2 = copy(ijk),
#    H_z_V2 = copy(ijk)]
#    return temp_out
# end

function add_I_u_p(I::AbstractArray, p::NamedTuple)
    temp_out = (I = I, r = I_u(I, p))
    return merge(p, temp_out)
end

# retina

function I_u(I::AbstractArray, p::NamedTuple)
    return I - imfilter(I, p.k_gauss_1, p.filling)
end

# lgn feedback

# todo: test
# todo: should lgn_a be normalized, ie divide by k??

function fun_x_lgn(x::AbstractArray, p::NamedTuple)
   # todo: change to abstract array? or is eltype doing that??
#     x_lgn = Array{eltype(x)}(0, p.dim_i, p.dim_j)
    x_lgn = zeros(p.dim_i, p.dim_j)
#     todo: change to map function?
    for k ∈  1:p.K
        x_lgn += x[:,:,k]
    end
    return x_lgn
end

# # lgn_A
# # todo: test
# function x_lgn_A(x_lgn, C1)
#     return C_1 .* x_A
# end
#
# # lgn_B
# # todo: test
# function fun_lgn_B(x_lgn::AbstractArray, σ_1, C_2)
# #     todo: alocate kernel
#     return C_2 .* imfilter(x_A, Kernel.gaussian(σ_1))
# end

function fun_F(value::Real, p::NamedTuple)
    max.(value - p.Γ, 0)
end

# williomson uses differnt F, relu with threshold
function fun_F_willimson(value::Real, p::NamedTuple)
    value < p.Γ ? zero(value) : value
end

# todo: check
function fun_f(x::AbstractArray, p::NamedTuple)
    return (p.μ .* x .^p.n) ./ (p.ν^p.n .+ x.^p.n)
end

function func_filter_W(img::AbstractArray, W::AbstractArray, p::NamedTuple)
#     out = reshape(fill(0/img[1,1,1], size(img)[1], size(img)[2] * p.K * p.K), size(img)[1], size(img)[2], p.K, p.K)
#     out = reshape(Array{eltype(img)}(undef, size(img)[1], size(img)[2] * p.K), size(img)[1], size(img)[2], p.K)
    out = copy(img)
    for k ∈ 1:p.K
#     todo fix W
        out[:,:,k] = imfilter(img[:,:,k], (centered(W[:,:,k,k]),), p.filling)
#         out[:,:,k] = img[:,:,k]
#         out[:,:,k] = conv(img[:,:,k], W[:,:,k,k])
        for l ∈ 1:p.K
            if l ≠ k
#                 out[:,:,k] += img[:,:,l]
                out[:,:,k] += imfilter(img[:,:,l], (centered(W[:,:,k,l]),), p.filling)
#                 out[:,:,k] .+= conv(img[:,:,l], W[:,:,k,l])
            end
        end
    end
    return out
end

# LGN
function fun_dv(v::AbstractArray, u::AbstractArray, x_lgn::AbstractArray, p::NamedTuple)
#     x_lgn = fun_x_lgn(x)
    return p.δ_v .* ( .- v .+
            ((1 .- v) .* max.(u,0) .* (1 .+ p.C_1 .* x_lgn)) .-
            ((1 .+ v) .* p.C_2 .* imfilter(x_lgn, (centered(p.k_gauss_1),), p.filling)))
#             ((1 .+ v) .* p.C_2 .* conv(x_lgn, p.k_gauss_1)))
end

function fun_dv!(du::AbstractArray, v::AbstractArray, u::AbstractArray, x_lgn::AbstractArray, p::NamedTuple)
#     x_lgn = fun_x_lgn(x)
    du .= p.δ_v .* ( .- v .+
            ((1 .- v) .* max.(u,0) .* (1 .+ p.C_1 .* x_lgn)) .-
            ((1 .+ v) .* p.C_2 .* imfilter(x_lgn, (centered(p.k_gauss_1),), p.filling)))
#             ((1 .+ v) .* p.C_2 .* conv(x_lgn, p.k_gauss_1)))
end

# lgn to l6 and l4
function fun_v_C(v_p::AbstractArray, v_m::AbstractArray, p::NamedTuple)
    # isodd(l) || throw(ArgumentError("length must be odd"))

    V = exp(-1/8) .* (imfilter((max.(v_p,0)-max.(v_m,0)), (centered(p.k_gauss_2),), p.filling))
#     V = exp(-1/8) .* (conv((max.(v_p,0) .- max.(v_m,0)), p.k_gauss_2))

# todo: change to abstract array? or is eltype doing that??
    A = reshape(Array{eltype(V)}(undef, size(V)[1], size(V)[2]*p.K),size(V)[1],size(V)[2],p.K)
    B = copy(A)

    for k ∈ 1:p.K
#         θ = π*(k-1)/p.K
        A[:,:,k] = imfilter(V, (centered(p.k_C_d[:,:,k]),), p.filling)
        B[:,:,k] = abs.(imfilter(V, (centered(p.k_C_b[:,:,k]),), p.filling))
#         A[:,:,k] = conv(V, p.k_C_d[:,:,k])
#         B[:,:,k] = abs.(conv(V, p.k_C_b[:,:,k]))
    end

    return p.γ .* (max.(A .- B, 0) .+ max.(.- A .- B, 0))
end

# L6
function fun_dx_V1(x::AbstractArray, C::AbstractArray, z::AbstractArray, x_v2::AbstractArray, p::NamedTuple)
    return p.δ_c .* (-x .+
            ((1 .- x) .*
                ((p.α*C) .+ (p.ϕ .* max.(z .- p.Γ,0)) .+ (p.v_21 .* x_v2) .+ p.att)))
end

#     L4 excit
function fun_dy(y::AbstractArray, C::AbstractArray, x::AbstractArray, m::AbstractArray, p::NamedTuple)
    return  p.δ_c .* (   .-y .+
            ((1 .- y) .* (C .+ (p.η_p .* x))) .-
            ((1 .+ y) .* fun_f(m .* func_filter_W(m, p.k_W_p, p), p)))
end

# function fun_dy(y::AbstractArray, C::AbstractArray, x::AbstractArray, m::AbstractArray, p::NamedTuple)
#     return  p.δ_c .* (   .-y .+
#             ((1 .- y) .* (C .+ (p.η_p .* x))) .-
#             ((1 .+ y) .* fun_f(m , p)))
# end

#     l4 inhib
function fun_dm(m::AbstractArray, x::AbstractArray, p::NamedTuple)
    return p.δ_m .* (  -m .+
                     (p.η_m .* x) .-
                     (m .* fun_f(func_filter_W(m, p.k_W_m, p), p)))
end

# function fun_dm(m::AbstractArray, x::AbstractArray, p::NamedTuple)
#     return p.δ_m .* (  -m .+
#                      (p.η_m .* x) .-
#                      (m .* fun_f(m, p)))
# end

#     L2/3 excit
function fun_dz(z::AbstractArray, y::AbstractArray,  H_z::AbstractArray, s::AbstractArray, p::NamedTuple)
    return p.δ_z .*   (-z .+
                    ((1 .- z) .*
#                         ((p.λ .* max.(y,0)) .+ H_z )) .-
                        ((p.λ .* max.(y,0)) .+ H_z .+ (p.a_23_ex .* p.att))) .-
                    ((z .+ p.ψ) .* (imfilter(s, (centered(p.k_T_p),), p.filling))))
#                     ((z .+ p.ψ) .* (conv(s, p.k_T_p))))
end

#     L2/3 inhib
function fun_ds(s::AbstractArray, z::AbstractArray, H_z::AbstractArray, p::NamedTuple)
    return p.δ_s .*   ( -s .+
                    H_z .+ (p.a_23_in .* p.att) .-
                    (s .* imfilter(s, (centered(p.k_T_m),), p.filling)))  #?????
#                     (s .* conv(s, p.k_T_m)))  #?????
end

function fun_H_z(z::AbstractArray, p::NamedTuple)
    H_z_out = copy(z)
    for k ∈ 1:p.K
#     fix H_z temp!!
#         H_z_out[:,:,k] = imfilter((max.(z[:,:,k] .- p.Γ,0)), (centered(p.k_H[:,:,k]),), p.filling)
        H_z_out[:,:,k] = max.(z[:,:,k] .- p.Γ,0)
#         H_z_out[:,:,k] = conv((max.(z[:,:,k] .- p.Γ,0)), p.k_H[:,:,k])
        end
        return H_z_out
end

#     V2 L6
function fun_dx_v2(x_v2::AbstractArray, z_v2::AbstractArray, z::AbstractArray, p::NamedTuple)
    return p.δ_c .*   (  -x_v2 .+
                    ((1 .- x_v2) .*
                        ((p.v_12_6 .* max.(z .- p.Γ, 0)) + (p.ϕ .* max.(z_v2 .- p.Γ, 0)) .+ p.att)))
end

# V2 L4 excit
function fun_dy_v2(y_v2::AbstractArray, z::AbstractArray, x_v2::AbstractArray, m_v2::AbstractArray, p::NamedTuple)
    return δ_c .*   (-y_v2 .+
                    ((1 .- y_v2) .* ((v_12_4 .* max.(z .- p.Γ, 0)) .+ (p.η_p .* x_v2))) .-
                    ((1 .+ y_v2) .* fun_f.(imfilter(m_v2, (centered(p.k_W_p),), p.filling)), p.μ, p.ν, p.n))
#                     ((1 .+ y_v2) .* fun_f.(conv(m_v2, p.k_W_p)), p.μ, p.ν, p.n))
end

# for equilabrium
fun_equ(x) = x/(1+x)

# lgn
function fun_v_equ(u::AbstractArray, x::AbstractArray, p::NamedTuple)
    return fun_equ.(relu.(u) .* (p.lgn_equ_u .+ (p.lgn_equ_A .* fun_A.(x))) .- (p.lgn_equ_B .* fun_B.(x)))
end

# lgn, no L6 feedback, light
function fun_v_equ_noFb(u::AbstractArray, x::AbstractArray, p::NamedTuple)
    return fun_equ.(relu.(u) .* (p.lgn_para_u))
end

# l6
function fun_x_equ(C::AbstractArray, z::AbstractArray,  p::NamedTuple)
    return fun_equ.((p.α .* C) .+ (p.ϕ .* fun_F.(z,p.Γ)) + (p.v_21 .* x_v2) .+ p.att)
end

# l6, no V2 feedback, light
function fun_x_equ_noV2(C::AbstractArray, z::AbstractArray, p::NamedTuple)
    return fun_equ.((p.α .* C) .+ (p.ϕ .* fun_F.(z,p.Γ)))
end

# l4 excit
function fun_y_equ(C::AbstractArray, x::AbstractArray,  m::AbstractArray, p::NamedTuple)
    return fun_equ.(C .+ (η_p .* x - fun_f.(imfilter(m, p.k_W_p), p.μ, p.ν, p.n)))
end

# l4 inhib, needs initial condition of itself
function fun_m_equ(C::AbstractArray, x::AbstractArray, m_init::AbstractArray, p::NamedTuple)
    return (p.η_m .* x ./ (1 .+  fun_f.(imfilter(m_init, p.k_W_m), p.μ, p.ν, p.n)))
end

# l4 inhib - no feedback kernel
function fun_m_equ_noFb(C::AbstractArray, x::AbstractArray, p::NamedTuple)
    return (p.η_m .* x)
end

#  l2/3 excit, needs initial condition of itself
function fun_z_equ(y::AbstractArray, s::AbstractArray,  z_init::AbstractArray,  p::NamedTuple)
    return (p.λ .* relu.(y)) .+ imfilter(fun_F.(z_init, p.Γ), k_H) .+ (p.a_23_ex .* p.att) .- (p.ϕ .* imfilter(s, k_T_p)) ./
    (1 .+ (p.λ .* relu.(y)) .+ imfilter(fun_F.(z_init, p.Γ), k_H) .+ (p.a_23_ex .* p.att) .+ imfilter(s, k_T_p))
    return
end

#  l2/3 excit - no feedback kernel
function fun_z_equ_noFb(y::AbstractArray, s::AbstractArray, p::NamedTuple)
    return (p.λ .* relu.(y)) .- (p.ϕ .* imfilter(s, k_T_p)) ./
    (1 .+ (p.λ .* relu.(y)) .+ imfilter(s, k_T_p))
    return
end

#  l2/3 inhib, needs initial condition of itself
function fun_s_equ(z::AbstractArray, s::AbstractArray,  s_init::AbstractArray, p::NamedTuple)
    return ((imfilter(fun_F.(z, p.Γ), k_H) + p.a_23_in .* p.att ) ./ (1 .+ imfilter(s_init, k_T_m)))           #????????? is T right?
end

#  l2/3 inhib - no feedback kernel
function fun_s_equ_noFb(z::AbstractArray, s::AbstractArray, p::NamedTuple)
    return (imfilter(fun_F.(z, p.Γ), k_H))
end

# V2 L6
function fun_xV2_equ(x_v2::AbstractArray, z::AbstractArray, z_v2::AbstractArray, p::NamedTuple)
    return fun_equ.((p.v_12_6 .* fun_F.(z, p.Γ) .+ (p.ϕ .* fun_F.(z_v2, p.Γ))))
end

#  V2 L4 excit
function fun_yV2_equ(z::AbstractArray, x::AbstractArray,  m::AbstractArray, p::NamedTuple)
    return (p.v_12_4 .* fun_F.(z,p.Γ) .+ (p.η_p .* x) - fun_f.(imfilter(m, k_W_p), p.μ, p.ν, p.n) ./
   (1 .+ p.v_12_4 .* fun_F.(z,p.Γ) .+ (p.η_p .* x) .+ fun_f.(imfilter(m, k_W_p), p.μ, p.ν, p.n)))
end

end

# dv_p = δ_v( -v_p +
#             ((1 - v_p) * relu(u) * (1 + lgn_A)) -
#             ((1 + v_p) * lgn_B))
#
#
# dv_m = δ_v( -v_m +
#             ((1 - v_m) * relu(-u) * (1 + lgn_A)) -
#             ((1 + v_m) * lgn_B))
#
# # feedback to lgn
# # todo: lgn_A
# lgn_A = C_1 * imfilter(x, kern_sumk(K))
#
# lgn_B = C_2 * imfilter(x, Kernel.gaussian(σ_1), p.filling)

#
# # L6
#
# dx = δ_c(   -x +
#             ((1 -x) *
#                 ((α*C) + (ϕ * relu.(z - Γ)) + (V_21 * x_v2) + att)))
#
# # L4
#
# dy = δ_c(   -y +
#             ((1 - y) * (C + (η_p * x))) -
#             ((1 + y) * fun_f((m .* imfilter(m, kern_W_p)))))
#
#
# # williamson - eq 20, relu? unclear
# # dy = δ_c(   -y +
# #             ((1 - y) * relu.((C + (η_p * x)))) -
# #             (y + 1) * sigmoid((m .* imfilter(m, W_p))))
#
#
# dm = δ_m(   -m +
#             (η_m * x) -
#             (m .* fun_f(imfilter(m, kern_W_m))))
#
#
# # L2/3
#
# dz = δ_z(   -z +
#             ((1 - z) *
#                 ((λ * relu(y)) + imfilter((relu(z - Γ)), H) + (a_ex_23 * att))) -
#             ((z + ψ) .* (imfilter(s, T_p))))
#
#
# # todo: sTs
# ds = δ_s(   -s +
#             imfilter((relu. (z - Γ)), H) + (a_23_in .* att) -
#             (s .* imfilter(s, T_m)))  #?????
#
#
#
#
# # V1 to V2
#
# # V2 L6
#
# dx_v2 = δ_c(    -x_v2 +
#                 (1 - x_v2) *
#                     ((v_12_6 * relu.(z - Γ)) + (ϕ * relu.(z_v2 - Γ)) + att))
#
#
# # V2 L4 excit
#
# dy_v2 = δ_c(    -y_v2 +
#                 ((1 - y_v2) * ((v_12_4 * relu.(z - Γ)) + (η_p * x_v2))) -
#                 ((1 + y_v2) * fun_f(imfilter(m_v2, W_p))))
#
#
# # V2 L4 inhib
#
#
# dm_v2 = δ_m(    -m_v2 +
#                 (η_m * x_v2) -
#                 (m_v2 .* fun_f(imfilter(m_v2, W_m))))
#
#
# # L2/3
#
# dz_v2 = δ_z(    -z_v2 +
#                 ((1 - z_v2) *
#                     ((λ * relu(y_v2)) +
#                         imfilter((relu.(z_v2 - Γ)), H) + a_ex_23 * att)) -
#                 ((z_v2 + ψ) .* (imfilter(s_v2, T_p))))
#
#
# # todo: sTs
# ds_v2 = δ_s(    -s_v2 +
#                 imfilter((fun_F(z_v2,Γ)),H_v2) -
#                 (s_v2 .* imfilter(s_v2, T_m)))  #?????


text:
export plotting, plotting!, plotting_colored_mutations, plotting_colored_mutations!, colors_by_mutations
export plotting_sampletumor_pies!, plotting_sampletumor_pies

function plotting!(scene, tumor;
        path="", markersize = 1.0, color=nothing, colormap=distinguishable_colors(68), alpha = 1.,
        plotargs...
        )

    isempty(tumor) && return
    p = (tumor isa DataFrameRow) ? [Point(tumor.position...)] : [Point(pos...) for pos in tumor.position]

    colors = if isnothing(color)
        [isempty(mut) ? 1 : mod(mut[end], length(colormap)) + 1 for mut in tumor.mutations]
    elseif !(color isa Vector)
        (color isa Number ? colormap[color] : color, alpha)
    elseif length(color) == nrow(tumor)
        [(c isa Number ? colormap[c] : c, alpha) for c in color]
    else
        error("Options for color: nothing, single value or value for every cell.")
    end

    meshscatter!(scene, p; markersize = markersize, color = colors, plotargs...)
    isempty(path) || save(path, scene)
    return scene
end
function plotting(tumor;
        size=(500,500), show_axis=false, inline=false, plotargs...)

    Makie.inline!(inline)
    scene = Scene(resolution=size, show_axis=show_axis)
    plotting!(scene, tumor; plotargs...)
end

function colors_by_mutations(tumor; colorpalette = palette(:tab20), sub = 1, autodepth = false, show_warning = true)
    points = [Point(pos...) for pos in tumor.position]
    colors = fill(colorpalette[1], length(points))
    subclones = clones(tumor; sub=sub, autodepth=autodepth, show_warning = show_warning)
    for (i,cl) in enumerate(subclones)
        colors[findall(in(cl.index), tumor.index)] .= colorpalette[(i-1)%length(colorpalette)+1]
    end
    return points, colors
end

function plotting_colored_mutations!(scene, tumor;
        path="", limits =  Makie.automatic,
        markersize = 1.0, colorpalette = palette(:tab20), shading = false,
        sub=1, autodepth = false, show_warning = true, alpha = 1., plotargs...
        )

    isempty(tumor) && return
    points, colors = colors_by_mutations(tumor; colorpalette = colorpalette, sub = sub, autodepth = autodepth, show_warning = show_warning)

    newscene = meshscatter!(scene, points, markersize = markersize, color = [(c,alpha) for c in colors], scale_plot = false, shading=shading, limits = limits, plotargs...)

    isempty(path) || save(path, scene)
    return scene
end

function plotting_colored_mutations(tumor;
        size=(500,500), show_axis=false, inline=false, plotargs...)
    Makie.inline!(inline)
    scene = Scene(resolution=size, show_axis=show_axis)
    plotting_colored_mutations!(scene, tumor; plotargs...)
end

function plotting_sampletumor_pies!(scene, sampletumor;
        path="", colorpalette = palette(:tab20), overdraw=false,
        strokecolor=nothing, strokewidth = 1, sample_r = 0., wt_color = :lightgrey,
        mutations = nothing, first = isnothing(mutations) ? 10 : length(mutations),
        plotargs...
        )

    top_muts = isnothing(mutations) ? sort( sampletumor_mfreqs(sampletumor), :frequency, rev=true).mutation[1:first,:] : mutations
    m_colors = Dict( top_muts .=> colorpalette[1:first])

    for s in eachrow(sampletumor)
            mask = map( m-> m in top_muts, s.mutations)
            center, freqs, colors = s.position, s.frequencies[mask], [m_colors[m] for m in s.mutations[mask]]
            if isempty(freqs)
                    freqs, colors = [1.],[wt_color]
            end
            freqs ./= sum(freqs)
            startphis = cumsum(freqs).-freqs
            sections = map(zip(startphis,freqs)) do (fs,f)
                    r = iszero(sample_r) ? s.sample_r : sample_r
                    circ = Point2f0[ center .+ r .* (cos(2π*p), sin(2π*p)) for p in (0.:0.002:f).+fs]
                    !isone(length(freqs)) && push!(circ,Point2f0(center))
                    circ
            end
            for (sec, c) in zip(sections, colors)
                poly!(scene, sec, color = c, strokewidth = strokewidth,
                    strokecolor = isnothing(strokecolor) ? c : strokecolor, overdraw = overdraw,
                    plotargs...
                )
            end
    end
    isempty(path) || save(path, scene)
    return scene
end

function plotting_sampletumor_pies(sampletumor;
        size=(500,500), show_axis=false, inline=false, plotargs...)
    Makie.inline!(inline)
    scene = Scene(resolution=size, show_axis=show_axis)
    plotting_sampletumor_pies!(scene, sampletumor; plotargs...)
end

max_stars_repo_name: aangaji/TumorGrowth

max_issues_repo_name: aangaji/TumorGrowth

max_forks_repo_name: aangaji/TumorGrowth


text:
"""
    check_blocksize(f::IOStream, position_before::Integer, blocksize_before::Integer)

Checks for integer overflow in the size of the block.
"""
function check_blocksize(f::IOStream, position_before::Integer, blocksize_before::Integer)

    seek(f,position_before+blocksize_before+12)
    blocksize_after = read(f,UInt32)

    # of the numbers are the same everything works as it should
    if blocksize_before == blocksize_after
        return blocksize_before
    else
        blocksize_before_fix = blocksize_before + 4294967296
        seek(f,position_before+blocksize_before_fix+12)
        blocksize_after_fix = read(f,UInt32) + 4294967296

        if blocksize_before_fix == blocksize_after_fix
            return blocksize_before_fix
        else
            error("There is an issue with the snapshot:\n
                   Blocksize_before = $blocksize_before_fix\n
                   Blocksize_after  = $blocksize_after_fix")
        end
    end
end

"""
    read_bockname(f::IOStream)

Reads the name of the Format 2 block.
"""
function read_bockname(f::IOStream)
    name = Char.(read!(f, Array{Int8,1}(undef,4)))
    return strip(String(name))
end

"""
    print_blocks(filename::String; verbose::Bool=true)

Reads the block names of blocks in a snapshot and returns them in an array.
Outputs them to console if `verbose=true`
"""
function print_blocks(filename::String; verbose::Bool=true)

    f = open(filename)
    blocksize = read(f, Int32)

    if blocksize != 8
        error("Block search not possible - use snap_format 2!")
    end

    p = position(f)
    seek(f,4)

    blocks = []

    while eof(f) != true

        # read block name
        blockname = read_bockname(f)

        # store block name in array
        push!(blocks, blockname)

        p = position(f)

        seek(f,p+8)

        skipsize = read(f, UInt32)

        skipsize = check_blocksize(f, p, skipsize)

        seek(f,p+skipsize+20)

    end

    if verbose
        println("Found blocks: ")
        for block ∈ blocks
            println(block)
        end
    end

    close(f)

    return String.(blocks)
end

"""
    block_present(filename::String, blockname::String, blocks::Vector{String}=[""])

Checks if a given block is present in the snapshot file, or in the supplied `blocks` Vector.
"""
function block_present(filename::String, blockname::String, blocks::Vector{String}=[""])

    # if no blocks supplied, read in block information
    if blocks == [""]
        blocks = print_blocks(filename, verbose=false)
    end

    # loop over all blocks
    for block ∈ blocks
        if block == blockname
            return true
        end
    end

    return false
end

"""
    check_blocks(snap_base::String, blocks::Array{String})

Check if all requested blocks are present.
"""
function check_blocks(filename::String, blocks::Array{String})

    # check if requested blocks are present
    no_mass_block = false
    blocks_in_file = print_blocks(filename, verbose=false)
    for blockname in blocks
        if !block_present(filename, blockname, blocks_in_file)
            if blockname == "MASS"

                no_mass_block = true
                deleteat!(blocks, findfirst(blocks .== "MASS"))
            else
                error("Block $blockname not present!")
            end
        end
    end

    return blocks, no_mass_block

end

"""
    get_block_positions(filename::String)

Returns a dictionary with the starting positions of all blocks in a snapshot in bits.
"""
function get_block_positions(filename::String)

    f = open(filename)
    blocksize = read(f, Int32)

    # only works for snap format 2
    if blocksize != 8
        error("Block search not possible - use snap_format 2!")
    end

    # allocate data dict
    d = Dict{String, Integer}()

    while eof(f) != true

        # read block name
        blockname = read_bockname(f)

        p = position(f)

        seek(f,p+8)

        skipsize = read(f, UInt32)

        skipsize = check_blocksize(f, p, skipsize)

        # store blockname and position in Dict
        d[blockname] = p+12

        # skip to the next name block
        seek(f,p+skipsize+20)

    end

    close(f)

    return d
end

"""
    select_file(filebase::String, filenum::Integer)

Checks if the requested files exists and returns the path.
"""
function select_file(filebase::String, filenum::Integer)

    if !isfile(filebase)
        filename = filebase * ".$filenum"
        if !isfile(filename)
            error("File $filename not present!")
        else
            return filename
        end
    else
        return filebase
    end
end

"""
    shift_across_box_border(x::Real, x_halo::Real, boxsize::Real, boxsize_half::Real)

Shift coordinate `x` across the box border if the zero coordinate `x₀` is on the other side.
"""
@inline function shift_across_box_border(x::Real, x₀::Real, boxsize::Real, boxsize_half::Real)
    if x - x₀ > boxsize_half
        return x - boxsize
    elseif x₀ - x > boxsize_half
        return x + boxsize
    end
    return x
end


text:
#=
This is free and unencumbered software released into the public domain.

Anyone is free to copy, modify, publish, use, compile, sell, or
distribute this software, either in source code form or as a compiled
binary, for any purpose, commercial or non-commercial, and by any
means.

In jurisdictions that recognize copyright laws, the author or authors
of this software dedicate any and all copyright interest in the
software to the public domain. We make this dedication for the benefit
of the public at large and to the detriment of our heirs and
successors. We intend this dedication to be an overt act of
relinquishment in perpetuity of all present and future rights to this
software under copyright law.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR
OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
OTHER DEALINGS IN THE SOFTWARE.

For more information, please refer to <http://unlicense.org>
=#

"""
TrackingHeaps offers a heap with a tracking system for the stored values.

Inserting a value into a TrackingHeap returns a tracker for the value. The
tracker can be used to access, update, and delete the value without searching
for it first. Heap order do not allow for `O(log m)` search (where `m` is the
number of values currently stored inside the heap), just for `O(m)` search, so
this feature allow for some performance gain if you need to manipulate values
anywhere in the heap (not just on top of the heap). Besides access, which is
`O(1)`, update and delete are `O(log m)` as they may need to rebalance the
tree.

If the tracking system is not needed, there is little reason to use this heap.

I wrote this package because the MutableBinaryHeap of DataStructures.jl did not
allow some behavior I wanted; behavior as:
1) a non-top value can be deleted without being made top first;
2) convenience methods allow to `pop!`/`delete!` stored values and
   immediatelly `track!` others, avoiding double-rebalancing;
3) after a value is deleted, its tracker can be re-reused to re-insert
   that value or insert a new value (but this is not done automatically);
4) the arity of the heap (binary, trinary, etc..) can be defined by the
   user (by parametric type) and inccur in minimal overhead (I hope);
5) all the stored values are in a `Vector{V}` in heap order, for easy
   backdoor/hacking access;
6) the integer type that is the tracker type can be defined by the user.
"""
module TrackingHeaps
# TODO: To allow automatic reuse create a different heap, as the struct would
#       need an extra vector for the "already used but expired" trackers.
#       Such vector could be reset if the heap becomes empty at any time,
#       and values removed from its back do not need to be added there.
# TODO: doctests (https://docs.julialang.org/en/latest/manual/documentation/)
# TODO: (to be considered) track the largest tracker not expired, when it is
#       removed search for the second largest, and resize tk2ix accordingly.
#       This can make a single removal O(num calls to track! already done), but
#       if we search using tk2ix, then it is amortized (for each extra step in
#       some deletion, there should have been a previous track!, and the extra
#       steps consume such "savings" made by previous track!). This assumes
#       the user will not use setindex! to insert values with large
#       tracker keys.
# TODO: for now, tracker keys cannot be of a !isbits type,
#       because for efficiency/convenience we convert the N-arity to the type
#       of K and then dispatch on it with Val, but Val (or any parametric type)
#       cannot wrap !isbits types.

export TrackingHeap # the main data structure
export MinHeapOrder, MaxHeapOrder, is_higher_than # Order related exports
export SafeFromYourself, NoTrainingWheels # Safety types
export top, next_tracker, renew_tracker!, track!, update!, extract!
export extract_and_push!, extract_and_track!
export pop_and_track!, pop_once_and_track_many!
export empty!!
# This module also exports many other methods but it does so by extending
# the Base module functions, not only that, but the methods provided for
# Base allow many other functions to work over TrackingHeap as it is an
# abstractdict.

# Move these to a separated file or a minimalistic package:
@inline oneless(x) = x - one(x)
@inline onemore(x) = x + one(x)
@inline two(::Type{T}) where {T} = convert(T, 2)
@inline two(x :: T) where {T} = convert(T, 2)
@inline ispositive(x) :: Bool = x > zero(x)
# Anyway, they should not be exported.

# TODO: WHAT NEED TO BE FINISHED FOR FIRST COMMIT
# 3) write tests for everything
# 4) save to git, upload to a new repository
# 5) check version, call register to register on the julia General

"""
    MinHeapOrder

The default type for the type parameter `O` (i.e., Order) of TrackingHeap.
If used the heap will have the minimum value as top value.

See also: [`is_higher_than`](@ref), [`TrackingHeap`](@ref)
"""
struct MinHeapOrder; end

"""
    MaxHeapOrder

The alternative type for the type parameter `O` (i.e., Order) of TrackingHeap.
If used the heap will have the maximum value as top value.

See also: [`is_higher_than`](@ref), [`TrackingHeap`](@ref)
"""
struct MaxHeapOrder; end

"""
    is_higher_than(HeapOrderType, x, y) :: Bool

To define an ordering that is not the stored value type (i.e., `V`) default,
you can create a new empty type and extend the `is_higher_than` function with a
method that takes such empty type and returns if `x` should be higher in the
heap than `y` based on it (so you do not need to wrap the value type inside a
new type for a different ordering).

Note that `is_higher_than` defaults to `<` and `>` instead of `isless` and
`!isless` (that would be `>=`).

Also, if you do not want to provide a `<` and `>` for the stored value type,
you can instead just define this function for `::Type{MinHeapOrder}` (and/or
`::Type{MaxHeapOrder}`) and `x`-and-`y` of the specific types of the values
compared.

See also: [`MinHeapOrder`](@ref), [`MaxHeapOrder`](@ref),
[`TrackingHeap`](@ref)
"""
is_higher_than(::Type{MinHeapOrder}, x, y) :: Bool = x < y
is_higher_than(::Type{MaxHeapOrder}, x, y) :: Bool = x > y

"""
    SafeFromYourself

The default type for the type parameter `S` (i.e., Safety) of TrackingHeap.
Inform the heap methods to use `@assert` to check for any possible
inconsistencies, and throw KeyError tried to access/delete/update an
non-existent key.

See also: [`NoTrainingWheels`](@ref)
"""
struct SafeFromYourself; end

"""
    NoTrainingWheels

The alternative type for the type parameter `S` (i.e., Safety) of TrackingHeap.
Inform the heap methods avoid any checking theoretically giving the best speed.

See also: [`SafeFromYourself`](@ref)
"""
struct NoTrainingWheels; end

"""
    TrackingHeap{K, V, N, O, S} <: AbstractDict{K, V}

The type of a heap that has tracker/keys of type K (an integer), stored values
of type V, that is N-ary (binary, trinary, etc...), with the order defined by
O, and the safety level defined by S.

Equal values are allowed but, if more than one value could be the top value,
then any of the equal values may be there (they can yet be distinguished by
tracker).

The TrackingHeap implements (almost) all methods described in Dict interface,
and can be seen as a Dict (with the special property of always allowing for
fast access to the minimal/maximal value and its key). There is an AbstractHeap
of DataStructures, but it was designed for heaps without key, and then
an AbstractMutableHeap interface (for heaps with keys) was grafted into it.
Also, inheriting such interface would make necessary that this package
used DataStructures just for inheriting its abstract type.
"""
struct TrackingHeap{K, V, N, O, S} <: AbstractDict{K, V}
  svals :: Vector{V} # svals: stored values, define the ordering
  trcks :: Vector{K} # trcks: trck[i] is the tracker/key returned for svals[i]
  tk2ix :: Vector{K} # tk2ix: tk2ix[trck[i]] == i, finds svals based on tracker
end

"""
    TrackingHeap(tracking_heap) # copy constructor
    TrackingHeap(::Type{V}; kwargs = ...)

Constructs a TrackingHeap with default types for all type parameters
except V (i.e., the stored values, that has no obvious default).

This constructor also accepts a variety of keyword arguments which allow
to change the default type parameters and initialize the heap.

The constructor allow to initialize using: a values vector in heap order to
be owned by the heap (no overhead, nor checking); a values iterable
collection (the respective trackers should be assumed to be
`one(K):convert(K, length(number of initial values))`); a pairs iterable
collection. These options are mutually exclusive.

See also: [`MinHeapOrder`](@ref), [`SafeFromYourself`](@ref)

# Arguments
- `K`: The tracker keys type. Default: `typeof(length(T[]))`.
- `N`: The N-arity of the heap. Default: binary (i.e, 2).
- `O`: The ordering of the heap values. Default: MinHeapOrder.
- `S`: The safety level of the heap methods. Default: SafeFromYourself.
- `init_val_heap`: a Vector already in heap order. Default: empty.
- `init_val_coll`: an initial collections of values to be copied.
  Default: empty.
- `init_pairs`: an initial collection of pairs tracker-value to be copied.
  Default: empty.
"""
function TrackingHeap(
  ::Type{V};
  K = typeof(length(V[])),
  N = convert(K, 2),
  O = MinHeapOrder,
  S = SafeFromYourself,
  init_val_heap = Vector{V}(),
  init_val_coll = Vector{V}(),
  init_pairs = Vector{Pair{K, V}}()
) where {V}
  N = convert(K, N)
  if !isempty(init_val_heap)
    init_trcks = collect(one(K):convert(K, length(init_val_heap)))
    init_tk2ix = deepcopy(init_trcks)
    TrackingHeap{K, V, N, O, S}(init_val_heap, init_trcks, init_tk2ix)
  else
    h = TrackingHeap{K, V, N, O, S}(Vector{V}(), Vector{K}(), Vector{K}())
    if !isempty(init_val_coll)
      foreach(v -> track!(h, v), init_val_coll)
    elseif !isempty(init_pairs)
      foreach(((k, v),) -> h[k] = v, init_pairs)
    end
    h
  end
end
function TrackingHeap(
  heap :: TrackingHeap{K, V, N, O, S}
) where {K, V, N, O, S}
  deepcopy(heap) # fallbacks to copy
end

#################### THE INTERNAL HELPER METHODS FOLLOW ####################

"""
    parent_ix

Internal helper method. Do not use. Subject to change.
"""
@inline function parent_ix(::Val{N}, cix :: K) where {K, N}
  n = convert(K, N)
  div(cix + (n - two(K)), n)
end
# TODO: The inline below only work if the type of the tracker is Int, fix that.
@inline parent_ix(::Val{2}, cix) = cix >> 1

"""
    last_child_ix

Internal helper method. Do not use. Subject to change.
"""
last_child_ix(::Val{N}, pix :: K) where {K, N} = onemore(pix * convert(K, N))

"""
    sift_down!(heap, ix)

Internal use. May be documented in the future.
Avoid use unless strictly necessary.
"""
function sift_down!(
  heap :: TrackingHeap{K, V, N, O, S}, ix :: K
) where {K, V, N, O, S}
  svals, trcks, tk2ix = heap.svals, heap.trcks, heap.tk2ix

  # Save the value and its tracker before the sift.
  sval = svals[ix]
  trck = trcks[ix]

  # This will not change. Need to be rigth type to avoid instability.
  last_ix_heap = convert(K, length(svals))

  while true
    lc_ix = last_child_ix(Val(N), ix)
    fc_ix = lc_ix - oneless(N)
    # If value is already in the last layer/level (i.e., have no child), break.
    last_ix_heap < fc_ix && break
    lc_ix = min(lc_ix, last_ix_heap)
    min_child_value, min_child_ix = svals[fc_ix], fc_ix
    for curr_child_ix = onemore(fc_ix):lc_ix
      if is_higher_than(O, svals[curr_child_ix], min_child_value)
        min_child_ix = curr_child_ix
        min_child_value = svals[curr_child_ix]
      end
    end
    # If the highest child cannot raise above sval, then sval cannot go down.
    !is_higher_than(O, min_child_value, sval) && break

    # Otherwise, the smallest child goes up (i.e., takes value place as parent).
    svals[ix] = svals[min_child_ix]
    trcks[ix] = trcks[min_child_ix]
    tk2ix[trcks[ix]] = ix

    # Now we need to check if sval will stay at its new position, or descend
    # even more.
    ix = min_child_ix
  end

  # If sval is at the last layer/level or all of its childs are smaller than
  # it, then we can finally write it at its final position. Note that no time
  # was wasted copying it to candidate positions.
  svals[ix] = sval
  trcks[ix] = trck
  tk2ix[trck] = ix
end

"""
    sift_up!(heap, ix)

Internal use. May be documented in the future.
Avoid use unless strictly necessary.
"""
function sift_up!(
  heap :: TrackingHeap{K, V, N, O, S}, ix :: K
) where {K, V, N, O, S}
  svals, trcks, tk2ix = heap.svals, heap.trcks, heap.tk2ix

  # Save the sval and its tracker before the sift.
  sval = svals[ix]
  trck = trcks[ix]
  while ix > one(K) # while ix is not the root (i.e., we can go up yet)
    pix = parent_ix(Val(N), ix)
    # If sval cannot be above its parent, then sval must stay at ix.
    !is_higher_than(O, sval, svals[pix]) && break

    # Otherwise, sval should be at parent position or above, bring the
    # parent down.
    svals[ix] = svals[pix]
    trcks[ix] = trcks[pix]
    tk2ix[trcks[ix]] = ix

    # Now we need to check if sval will stay at its new position, or ascend
    # even more.
    ix = pix
  end

  # If sval is at root position or have a smaller parent, then we can
  # finally write it at its final position. Note that no time was wasted
  # copying it to candidate positions.
  svals[ix] = sval
  trcks[ix] = trck
  tk2ix[trck] = ix
end

"""
    sift!(heap, ix, previous_sval)

Internal use. May be documented in the future.
Avoid use unless strictly necessary.
"""
function sift!(
  heap :: TrackingHeap{K, V, N, O, S},
  ix :: K,
  previous_sval :: V
) where {K, V, N, O, S}
  if is_higher_than(O, previous_sval, heap.svals[ix])
    sift_down!(heap, ix)
  else
    sift_up!(heap, ix)
  end
end

#################### START OF THE EXPORTED API ####################

"""
    top(heap) -> (sval, trck)

Returns the "highest" value stored inside the heap and its tracker. `O(1)`.

See also: [`is_higher_than`](@ref)
"""
function top(
  heap :: TrackingHeap{K, V, N, O, S}
) :: Pair{K, V} where {K, V, N, O, S}
  Pair(heap.trcks[1], heap.svals[1])
end

"""
    next_tracker(heap) -> trck

Return the same tracker the next call of `track!` would return, without
modifying the heap in any way. `O(1)`.

See also: [`track!`](@ref)
"""
function next_tracker(
  heap :: TrackingHeap{K, V, N, O, S}
) :: K where {K, V, N, O, S}
  onemore(convert(K, length(heap.tk2ix)))
end

"""
    renew_tracker!(heap, old_trck, new_trck = next_tracker(heap)) -> new_trck

Marks `old_trck` as unused, begins using `new_trck` to refer to the
value previously pointed by `old_trck`, and returns new_trck.

Does not need any rebalancing, `O(1)`, but may inccur in memory allocation if
`new_trck` was never used. The `new_trck` can be a tracker already used in the
past, but must not be in use at the moment.

See also: [`pop_and_track!`](@ref), [`update!`](@ref), [`next_tracker`](@ref),
[`extract_and_track!`](@ref)
"""
function renew_tracker!(
  heap :: TrackingHeap{K, V, N, O, S},
  old_trck :: K
) :: K where {K, V, N, O, S}
  S === SafeFromYourself &&
    !haskey(heap, old_trck) &&
    throw(KeyError(old_trck))

  ix = heap.tk2ix[old_trck]
  push!(heap.tk2ix, ix)
  new_trck = convert(K, length(heap.tk2ix))
  heap.trcks[ix] = new_trck

  return new_trck
end
function renew_tracker!(
  heap :: TrackingHeap{K, V, N, O, S},
  old_trck :: K,
  new_trck :: K
) :: K where {K, V, N, O, S}
  S === SafeFromYourself &&
    !haskey(heap, old_trck) &&
    throw(KeyError(old_trck))
  S === SafeFromYourself && @assert !haskey(heap, new_trck)

  ix = heap.tk2ix[old_trck]
  heap.trcks[ix] = new_trck
  max_trck = length(heap.tk2ix)
  if new_trck == onemore(max_trck)
    push!(heap.tk2ix, ix)
  else
    new_trck > max_trck && resize!(heap.tk2ix, new_trck)
    heap.tk2ix[new_trck] = ix
  end

  return new_trck
end

"""
    track!(heap, sval) -> trck

Insert a new sval into the heap and return a new "never used" tracker for it.

The return is the same the method `next_tracker` would give if called before
`track!`.
The heap may need to be rebalanced, so `O(log m)`.

See also: [`pop_and_track!`](@ref), [`next_tracker`](@ref), [`empty!!`](@ref)
"""
function track!(
  heap :: TrackingHeap{K, V, N, O, S},
  sval :: V
) :: K where {K, V, N, O, S}
  S === SafeFromYourself && @assert length(heap.svals) == length(heap.trcks)
  push!(heap.svals, sval)
  ix = convert(K, length(heap.svals))
  push!(heap.tk2ix, ix)
  trck = convert(K, length(heap.tk2ix))
  push!(heap.trcks, trck)

	sift_up!(heap, ix)

  return trck
end

"""
    update!(heap, trck, new_value) -> heap
    update!(heap, trck => new_value) -> heap

Update the value pointed by `trck`. `O(log m)`.

Similar to `setindex!` but assumes `trck` exists, if it does not, gives a
KeyError when `S === SafeFromYourself`, and undefined behaviour when `S
!== SafeFromYourself`.

See also: [`setindex!`](@ref)
"""
function update!(
  heap :: TrackingHeap{K, V, N, O, S},
  trck :: K,
  new_value :: V
) :: TrackingHeap{K, V, N, O, S} where {K, V, N, O, S}
  S === SafeFromYourself && !haskey(heap, trck) && throw(KeyError(trck))

  ix = heap.tk2ix[trck]
  old_value = heap.svals[ix]
  heap.svals[ix] = new_value

  sift!(heap, ix, old_value)

  heap
end
function update!(
  heap :: TrackingHeap{K, V, N, O, S},
  elem :: Pair{K, V},
) :: TrackingHeap{K, V, N, O, S} where {K, V, N, O, S}
  trck, new_value = elem
  update!(heap, trck, new_value)
end

"""
    extract!(heap, trck) -> sval

Delete the value referred by `trck` in the `heap`, and return it.

The `trck` can then be used to re-insert it, or to insert another value.

The heap may need to be rebalanced, so `O(log m)`.

See also: [`delete!`](@ref), [`empty!`](@ref)
"""
function extract!(
  heap :: TrackingHeap{K, V, N, O, S},
  trck :: K
) :: V where {K, V, N, O, S}
  # Asserts and abbreviations.
  S === SafeFromYourself && !haskey(heap, trck) && throw(KeyError(trck))
  svals, trcks, tk2ix = heap.svals, heap.trcks, heap.tk2ix

  ix = tk2ix[trck]
  # Save deleted value for knowing how to sift the one that takes its place,
  # and also for returning it.
  del_sval = svals[ix]
  # If the heap has values after the one removed, then the gap left needs
  # to be filled, we use the last value to do so.
  if ix < length(heap)
    svals[ix] = last(svals)
    trcks[ix] = last(trcks)
    tk2ix[trcks[ix]] = ix
  end
  # Delete the last value, which was the one removed, or was already moved to
  # the position of the one deleted to fill the gap.
  pop!(svals)
  pop!(trcks)

  # There will only be a value to sift if the position of the extracted value
  # was not in the last position of the heap.
	ix < length(heap) && sift!(heap, ix, del_sval)

  del_sval
end

"""
    extract_and_push!(heap, old_trck, new_trck => new_sval) -> old_sval

Similar to `extract_and_track!` but let you define `new_trck`.

CAUTION: guarantee that the value referred by old_trck will be extracted and
the new pair will be pushed but does not guarantee that the internal heap
structure will end up in the same state as calling `extract!` followed by
`push!`.

See also: [`extract_and_track!`](@ref), [`extract!`](@ref)
"""
function extract_and_push!(
  heap :: TrackingHeap{K, V, N, O, S},
  old_trck :: K,
  new_pair :: Pair{K, V}
) :: V where {K, V, N, O, S}
  S === SafeFromYourself &&
    !haskey(heap, old_trck) &&
    throw(KeyError(old_trck))
  new_trck, new_sval = new_pair

  # TODO: this can be more optimized by checking if the new_trck is in use
  # and if so, check which one of old_trck and new_trck is the best to update!
  old_sval = heap[old_trck]
  update!(heap, old_trck, new_sval)
  if old_trck != new_trck
    haskey(heap, new_trck) && delete!(heap, new_trck)
    renew_tracker!(heap, old_trck, new_trck)
  end

  return old_sval
end

"""
    extract_and_track!(heap, old_trck, new_sval) -> (old_sval, new_trck)

Similar to calling `extract!` followed by `track!` but more optimized.

CAUTION: guarantee that the value referred by old_trck will be extracted and
the new value will be tracked but does not guarantee that the internal heap
structure will end up in the same state as calling `extract!` followed by
`track!`.

See also: [`pop_and_track!`](@ref), [`extract!`](@ref), [`track!`](@ref),
"""
function extract_and_track!(
  heap :: TrackingHeap{K, V, N, O, S},
  old_trck :: K,
  new_sval :: V
) :: Tuple{V, K} where {K, V, N, O, S}
  S === SafeFromYourself && !haskey(heap, old_trck) &&
    throw(KeyError(old_trck))

  old_sval = heap[old_trck]
  update!(heap, old_trck, new_sval)
  new_trck = renew_tracker!(heap, old_trck)

  return old_sval, new_trck
end

"""
    pop_and_track!(heap, new_sval) -> ((top_sval, top_trck), new_trck)

Similar to calling `pop!` followed by `track!` but more optimized.

CAUTION: guarantee that the top will be popped and the new value will be
tracked but does not guarantee that the internal heap structure will end up in
the same state as calling `pop!` followed by `track!`.

See also: [`pop!`](@ref), [`track!`](@ref),
"""
function pop_and_track!(
  heap :: TrackingHeap{K, V, N, O, S},
  new_sval :: V
) :: Tuple{Pair{K, V}, K} where {K, V, N, O, S}
  old_pair = top(heap)
  old_trck = first(old_pair)
  update!(heap, old_trck, new_sval)
  new_trck = renew_tracker!(heap, old_trck)

  return old_pair, new_trck
end

"""
    pop_once_and_track_many!(heap, new_svals)::((top_sval,top_trck), new_trcks)

Similar to calling `pop!` once followed by many calls to `track!` but more
optimized.

CAUTION: guarantee that the top will be popped and the new values will be
tracked but does not guarantee that the internal heap structure will end up in
the same state as it would by calling `pop!` followed by calls to `track!` (in
the same order given in the array).

See also: [`pop!`](@ref), [`track!`](@ref), [`pop_and_track!`](@ref)
"""
function pop_once_and_track_many!(
  heap :: TrackingHeap{K, V, N, O, S},
  new_svals
) :: Tuple{Pair{K, V}, Vector{K}} where {K, V, N, O, S}
  if isempty(new_svals)
    popped_sval = pop!(heap)
    (popped_sval, K[])
  else
    S === SafeFromYourself && @assert eltype(new_svals) == V
    new_trcks = K[]
    sizehint!(new_trcks, length(new_svals))
    min_new_sval, min_new_sval_ix = findmin(new_svals)
    popped_sval, min_new_sval_trck = pop_and_track!(heap, min_new_sval)
    for new_sval_ix in eachindex(new_svals)
      if new_sval_ix == min_new_sval_ix
        new_trck = min_new_sval_trck
      else
        new_trck = track!(heap, new_svals[new_sval_ix])
      end
      push!(new_trcks, new_trck)
    end
    S === SafeFromYourself && @assert length(new_trcks) == length(new_svals)
    S === SafeFromYourself && @assert !isnothing(popped_sval)
    (popped_sval, new_trcks)
  end
end

function Base.empty(
  heap :: TrackingHeap{K, V, N, O, S}
) where {K, V, N, O, S}
  TrackingHeap(V, K = K, N = N, O = O, S = S)
end

# BELOW ARE ITERATION INTERFACE METHODS: iterate
# https://docs.julialang.org/en/v1/base/collections/#lib-collections-iteration-1
# IteratorSize and IteratorEltype are default.

function Base.iterate(
  heap :: TrackingHeap{K, V, N, O, S}
) :: Union{Nothing, Tuple{Pair{K, V}, K}} where {K, V, N, O, S}
  isempty(heap) && return nothing
  first_ix = one(K)
  (heap.trcks[first_ix] => heap.svals[first_ix], first_ix)
end

function Base.iterate(
  heap :: TrackingHeap{K, V, N, O, S},
  last_ix :: K
) :: Union{Nothing, Tuple{Pair{K, V}, K}} where {K, V, N, O, S}
  curr_ix = onemore(last_ix)
  curr_ix > length(heap) && return nothing
  (heap.trcks[curr_ix] => heap.svals[curr_ix], curr_ix)
end

# BELOW ARE GENERAL COLLECTION INTERFACE METHODS: isempty, empty!, and length
# https://docs.julialang.org/en/v1/base/collections/#General-Collections-1

"""
    isempty(heap) :: Bool

Returns true if the heap has no values stored; false otherwise. `O(1)`.
"""
function Base.isempty(
  heap :: TrackingHeap
) :: Bool
  isempty(heap.svals)
end

"""
    empty!(heap) -> heap

Delete all stored values efficiently (consider internals structure unused).
`O(1)`.

CAUTION: this does not reset which trackers are considered "never used",
so a call to `next_tracker` before `empty!` and another after will return
the same.

Do not guarantee that previous `sizehint!` or any vector growth have their
effect negated.

See also: [`extract!`](@ref), [`delete!`](@ref), [`track!`](@ref)
"""
function Base.empty!(
  heap :: TrackingHeap{K, V, N, O, S}
) :: TrackingHeap{K, V, N, O, S} where {K, V, N, O, S}
  empty!(heap.svals)
  empty!(heap.trcks)
  heap
end

"""
    empty!!(heap) -> heap

Delete all stored values efficiently (consider internals structure unused).
`O(1)`.

CAUTION: this resets which trackers are considered "never used", so a call to
`next_tracker` or `track!` after it will return the first tracker a newly
created TrackingHeap gives (and so on).

Do not guarantee that previous `sizehint!` or any vector growth have their
effect negated.

See also: [`extract!`](@ref), [`delete!`](@ref)
[`track!`](@ref)
"""
function empty!!(
  heap :: TrackingHeap{K, V, N, O, S}
) :: TrackingHeap{K, V, N, O, S} where {K, V, N, O, S}
  empty!(heap.svals)
  empty!(heap.trcks)
  empty!(heap.tk2ix)
  heap
end

"""
    length(heap) :: K

Returns the number of values currently stored inside the heap. `O(1)`.
"""
function Base.length(
  heap :: TrackingHeap{K, V, N, O, S}
) :: K where {K, V, N, O, S}
  convert(K, length(heap.svals))
end

# NO ITERABLE COLLECTION INTERFACE METHODS IS IMPLEMENTED.
# The eltype, indexin, in, ∈, ∋, ∉, ∌, unique, allunique, reduce, fold{l,r},
# sum, prod, any, all, count, foreach, map*, first, and collect are default.
# Base.{maximum,minimum,extrema} are default and note they work over pairs, so
# they search for the maximum/minimum tracker key.
# Base.{argmin,argmax,findmin,findmax} are default and give the expected value
# considering the values (and not the keys). They are not optimized as top.
# The destructive version of many of the previous methods (basically anything
# that interacted with AbstractArray) was not implemented.

# BELOW ARE INDEXABLE COLLECTION INTERFACE METHODS: getindex, setindex!,
# firstindex, and lastindex.

"""
    getindex(heap, trck) -> sval

Return the value referred by `trck` in `heap`. O(1).

Only check if the tracker key exists if `S === SafeFromYourself`.

See also: [`setindex!`](@ref)
"""
function Base.getindex(
  heap :: TrackingHeap{K, V, N, O, S},
  trck :: K
) :: V where {K, V, N, O, S}
  S === SafeFromYourself && !haskey(heap, trck) && throw(KeyError(trck))
  return heap.svals[heap.tk2ix[trck]]
end

"""
    setindex!(heap, new_value, trck) -> heap

Begin tracking new_value using trck if such tracker is not yet in use,
update the value pointed by trck if such tracker was already in use.

Will rebalance the heap if necessary. `O(log m)`. Note: a vector with length
equal to the highest used `trck` is kept, so giving `new_value` and arbitrarily
large tracker can cause massive (and unnecessary) memory use.

See also: [`update!`](@ref)
"""
function Base.setindex!(
  heap :: TrackingHeap{K, V, N, O, S},
  new_value :: V,
  trck :: K
) :: TrackingHeap{K, V, N, O, S} where {K, V, N, O, S}
  S === SafeFromYourself && @assert length(heap.svals) == length(heap.trcks)
  if haskey(heap, trck)
    update!(heap, trck, new_value)
  else
    push!(heap.svals, new_value)
    push!(heap.trcks, trck)
    ix = convert(K, length(heap.trcks))

    max_trck = convert(K, length(heap.tk2ix))
    if trck > max_trck
      if trck == onemore(max_trck)
        push!(heap.tk2ix, ix)
      else
        resize!(heap.tk2ix, trck)
        heap.tk2ix[trck] = ix
      end
    else
      heap.tk2ix[trck] = ix
    end

    sift_up!(heap, ix)
  end

  heap
end

function Base.firstindex(
  heap :: TrackingHeap{K, V, N, O, S}
) :: K where {K, V, N, O, S}
  first(heap.trcks)
end

function Base.lastindex(
  heap :: TrackingHeap{K, V, N, O, S}
) :: K where {K, V, N, O, S}
  last(heap.trcks)
end

# BELOW ARE DICT INTERFACE METHODS: haskey, get, get!, getkey, delete!, pop!,
# keys, values, pairs, merge, merge!, sizehint!, keytype, and valtype.
# No get* method protect against type instability nor intend to do so.

function Base.haskey(
  heap :: TrackingHeap{K, V, N, O, S},
  trck :: K
) where {K, V, N, O, S}
  S === SafeFromYourself && ispositive(trck)
  trcks, tk2ix = heap.trcks, heap.tk2ix # abbreviate
  # The length of tk2ix need to allow trck and, as tk2ix can have uninitialized
  # values (from resize), the value of tk2ix[trck] needs to be inside the
  # expected bounds (the number of values currently inside the heap) and,
  # finally, trcks[tk2ix[trck]] needs to agree on which is its tracker (if it
  # does not, it can be be an expired tracker).
  length(tk2ix) >= trck && tk2ix[trck] in one(K):convert(K, length(trcks)) &&
    trcks[tk2ix[trck]] == trck
end

function Base.get(
  heap :: TrackingHeap{K, V, N, O, S},
  trck :: K,
  default
) where {K, V, N, O, S}
  haskey(heap, trck) ? heap[trck] : default
end

function Base.get(
  f :: Function,
  heap :: TrackingHeap{K, V, N, O, S},
  trck :: K
) where {K, V, N, O, S}
  haskey(heap, trck) ? heap[trck] : f()
end

function Base.get!(
  heap :: TrackingHeap{K, V, N, O, S},
  trck :: K,
  default :: V
) where {K, V, N, O, S}
  haskey(heap, trck) ? heap[trck] : heap[trck] = default
end

function Base.get!(
  f :: Function,
  heap :: TrackingHeap{K, V, N, O, S},
  trck :: K
) where {K, V, N, O, S}
  haskey(heap, trck) ? heap[trck] : heap[trck] = f()
end

function Base.getkey(
  heap :: TrackingHeap{K, V, N, O, S},
  trck :: K,
  default
) where {K, V, N, O, S}
  haskey(heap, trck) ? trck : default
end

# TODO: Maybe, in the future, implement a `merge!` method specific for heaps
# that is more optimized than the fallback provided in abstractdict.jl
# NOTE: both `merge!` variants are implemented by abstractdict.jl, but both
# `merge` variants needs to be defined for the specific type to avoid returning a
# `Dict` instead of a `TrackingHeap`.
function Base.merge(
  h :: TrackingHeap{K, V, N, O, S},
  others :: AbstractDict...
) where {K, V, N, O, S}
  pK = Union{K, map(keytype, others)...}
  pV = Union{V, map(valtype, others)...}
  new_h = TrackingHeap(pV; K = pK, N = N, O = O, S = S)
  merge!(new_h, h, others...)
  new_h
end

function Base.merge(
  combine :: Function,
  h :: TrackingHeap{K, V, N, O, S},
  others :: AbstractDict...
) where {K, V, N, O, S}
  pK = Union{K, map(keytype, others)...}
  pV = Union{V, map(valtype, others)...}
  new_h = TrackingHeap(pV; K = pK, N = N, O = O, S = S)
  merge!(combine, new_h, h, others...)
  new_h
end

function Base.sizehint!(
  heap :: TrackingHeap{K, V, N, O, S}, size :: K
) where {K, V, N, O, S}
  sizehint!(heap.svals, size)
  sizehint!(heap.trcks, size)
  # It can be that tk2ix is (or will be) larger than the other two (because
  # after some insertal and removal of values without reusing keys it will
  # grow larger than the other two), but there is no harm in passing sizehint!
  # to it.
  sizehint!(heap.tk2ix, size)
end

function Base.sizehint!(
  heap :: TrackingHeap{K, V, N, O, S}, size
) where {K, V, N, O, S}
  sizehint!(heap, convert(K, size))
end

"""
    delete!(heap, trck) -> heap

Delete the value referred by `trck` in the `heap`, and return the `heap`.

The `trck` can then be used to re-insert it, or to insert another value.

The heap may need to be rebalanced, so `O(log m)`.

See also: [`extract!`](@ref), [`pop!`](@ref)
"""
function Base.delete!(
  heap :: TrackingHeap{K, V, N, O, S},
  trck :: K
) :: TrackingHeap{K, V, N, O, S} where {K, V, N, O, S}
  S === SafeFromYourself && !haskey(heap, trck) && throw(KeyError(trck))
  extract!(heap, trck)
  heap
end

"""
    pop!(heap) -> Pair(sval, trck)

Delete the stored value at top of the heap, and return the `tracker => value`.

The `trck` can then be used to re-insert it, or to insert another value.

The heap may need to be rebalanced, so `O(log m)`.

See also: [`extract!`](@ref), [`delete!`](@ref)
"""
function Base.pop!(
  heap :: TrackingHeap{K, V, N, O, S}
) :: Pair{K, V} where {K, V, N, O, S}
  S === SafeFromYourself && @assert !isempty(heap)
  old_top = top(heap)

  extract!(heap, first(old_top))

  old_top
end

Base.keys(heap :: TrackingHeap) = heap.trcks
Base.values(heap :: TrackingHeap) = heap.svals
function Base.pairs(heap :: TrackingHeap)
  (Pair(k, v) for (k, v) in zip(keys(heap), values(heap)))
end

Base.keytype(heap :: TrackingHeap{K, V, N, O, S}) where {K, V, N, O, S} = K
Base.valtype(heap :: TrackingHeap{K, V, N, O, S}) where {K, V, N, O, S} = V

end # module


text:
@discord_object struct Connection
    id::String
    name::String
    type::String
    revoked::Bool
    integrations::Vector{Integration}
    verified::Bool
    friend_sync::Bool
    show_activity::Bool
    visibility::ConnectionVisibility.ConnectionVisibilityEnum
end


text:
M = [-1 1 0 0
    0 -1 1 0
    0 0 -1 1
    1 0 0 -1
]
using LinearAlgebra

Λ, V = eigen(M)
Λᵀ, Vᵀ = eigen(M')

max_stars_repo_name: sandreza/GradientLearning

max_issues_repo_name: sandreza/GradientLearning

max_forks_repo_name: sandreza/GradientLearning


text:
# This file is a part of SIS3316Digitizers.jl, licensed under the MIT License (MIT).

const ch_buffer_size = 64 * 1024 * 1024

@enum DataTransferCmd::UInt8 begin
    DataTransferCtlReset        = 0 # Reset Transfer FSM
    DataTransferCtlReset_Also   = 1 #
    DataTransferCtlRead         = 2 # Start Read Transfer
    DataTransferCtlWrite        = 3 # Start Write Transfer
end
export DataTransferCmd
export DataTransferCtlReset, DataTransferCtlReset_Also, DataTransferCtlRead, DataTransferCtlWrite

@enum MemSpaceSel::UInt8 begin
    MemSpaceMem1                = 0 # Memory 1 (Ch. 1 and 2, resp. 5 and 6, etc.)
    MemSpaceMem2                = 1 # Memory 2 (Ch. 3 and 4, resp. 7 and 8, etc.)
    MemSpaceStatCounter         = 3 # Statistic Counter (128 32-bit words)
end
export MemSpaceSel
export MemSpaceMem1, MemSpaceMem2, MemSpaceStatCounter

const data_transfer_ctl_register_layout = (
    cmd               = BitRange{RW}(30:31, DataTransferCmd), # Cmd
    mem_space_sel     = BitRange{RW}(28:29, MemSpaceSel), # Space select bits
    mem_addr          = BitRange{RW}( 0:27), # Memory 32-bit start address (128 Meg x 16 : only address bits 25-0 are used )
)
export data_transfer_ctl_register_layout

# ADC FPGA Data Transfer Control Register (only write complete register!)
data_transfer_ctl_register(addr::Unsigned) = Register{RW,RegVal}(
    MemAddr(addr),
    data_transfer_ctl_register_layout
)

# ADC FPGA Data Transfer Control Register, Ch. 1 to 4
const data_transfer_ch1_4_ctrl_reg = data_transfer_ctl_register(0x80)
export data_transfer_ch1_4_ctrl_reg

# ADC FPGA Data Transfer Control Register, Ch. 5 to 8
const data_transfer_ch5_8_ctrl_reg = data_transfer_ctl_register(0x84)
export data_transfer_ch5_8_ctrl_reg

# ADC FPGA Data Transfer Control Register, Ch. 9 to 12
const data_transfer_ch9_12_ctrl_reg = data_transfer_ctl_register(0x88)
export data_transfer_ch9_12_ctrl_reg

# ADC FPGA Data Transfer Control Register, Ch. 13 to 16
const data_transfer_ch13_16_ctrl_reg = data_transfer_ctl_register(0x8C)
export data_transfer_ch13_16_ctrl_reg

# ADC FPGA Data Transfer Control Registers
const data_transfer_ctrl_reg = (
    data_transfer_ch1_4_ctrl_reg,
    data_transfer_ch5_8_ctrl_reg,
    data_transfer_ch9_12_ctrl_reg,
    data_transfer_ch13_16_ctrl_reg,
)
export data_transfer_ctrl_reg

function data_transfer_ctrl_regval(
    cmd::DataTransferCmd,
    mem_space_sel::MemSpaceSel,
    mem_addr::Unsigned
)
    v::RegVal = zero(RegVal)
    v = setval(v, data_transfer_ctl_register_layout.cmd, cmd)
    v = setval(v, data_transfer_ctl_register_layout.mem_space_sel, mem_space_sel)
    v = setval(v, data_transfer_ctl_register_layout.mem_addr, mem_addr)
    v
end
export data_transfer_ctrl_regval

# ADC FPGA Data Transfer Status Register
data_transfer_status_register(addr::Unsigned) = Register{RO,RegVal}(
    MemAddr(addr),
    (
        busy              = Bit{RO}(31),  # Data Transfer Logic busy
        direction         = Bit{RO}(30),  # Data Transfer Direction (Write-Flag; 0: Memory -> VME FPGA; 1: VME FPGA -> Memory)
        fifo_near_full    = Bit{RO}(28),  # FIFO (read VME FIFO) Data AlmostFull Flag
        max_no_pend_read  = Bit{RO}(27),  # "max_nof_pending_read_requests"
        no_pend_read      = Bit{RO}(26),  # "no_pending_read_requests"
        int_addr_counter  = BitRange{RO}(0:25),  # Data Transfer internal 32-bit Address counter
    )
)

# ADC FPGA Data Transfer Status Register, Ch. 1 to 4
const data_transfer_adc1_4_status_reg = data_transfer_status_register(0x90)
export data_transfer_adc1_4_status_reg

# ADC FPGA Data Transfer Status Register, Ch. 5 to 8
const data_transfer_adc5_8_status_reg = data_transfer_status_register(0x94)
export data_transfer_adc5_8_status_reg

# ADC FPGA Data Transfer Status Register, Ch. 9 to 12
const data_transfer_adc9_12_status_reg = data_transfer_status_register(0x98)
export data_transfer_adc9_12_status_reg

# ADC FPGA Data Transfer Status Register, Ch. 13 to 16
const data_transfer_adc13_16_status_reg = data_transfer_status_register(0x9C)
export data_transfer_adc13_16_status_reg

# ADC FPGA Data Transfer Status Registers
const data_transfer_status_reg = (
    data_transfer_adc1_4_status_reg,
    data_transfer_adc5_8_status_reg,
    data_transfer_adc9_12_status_reg,
    data_transfer_adc13_16_status_reg,
)
export data_transfer_status_reg

const dls_fpga_bits = BitGroup{RW}(
    frame_err_cl    = Bit{RW}(7),   # Write: ADC FPGA: Clear Frame_error_latch; Read: ADC FPGA: Frame_error_latch
    soft_err_cl     = Bit{RW}(6),   # Write: ADC FPGA: Clear Soft_error_latch; Read: ADC FPGA: Soft_error_latch
    hard_err_cl     = Bit{RW}(5),   # Write: ADC FPGA: Clear Hard_error_latch; Read: ADC FPGA: Hard_error_latch
    lane_up         = Bit{RO}(4),   # Read: ADC FPGA: Lane_up_flag
    ch_up           = Bit{RO}(3),   # Read: ADC FPGA: Channel_up_flag
    frame_err       = Bit{RO}(2),   # Read: ADC FPGA: Frame_error_flag
    soft_err        = Bit{RO}(1),   # Read: ADC FPGA: Soft_error_flag
    hard_err        = Bit{RO}(0),   # Read: ADC FPGA: Hard_error_flag
)

"""
    vme_fpga_link_adc_prot_status_register(addr::Unsigned)

VME FPGA – ADC FPGA Data Link Status register (R/W, 0xA0)

Note: Only write complete Register.
"""
const vme_fpga_link_adc_prot_status_register = Register{RW,RegVal}(
    MemAddr(0xA0),
    (
        adc = BitGroup{RW}(
            fpga1 = dls_fpga_bits <<  0,  # FPGA 1 bits
            fpga2 = dls_fpga_bits <<  8,  # FPGA 2 bits
            fpga3 = dls_fpga_bits << 16,  # FPGA 3 bits
            fpga4 = dls_fpga_bits << 24,  # FPGA 4 bits
        ),
    )
)
export vme_fpga_link_adc_prot_status_register

fpga_num(dev_ch::Int) = div((dev_ch - 1), 4) + 1
fpga_ch(dev_ch::Int) = rem((dev_ch - 1), 4) + 1
export fpga_num, fpga_ch

const data_transfer_ctl_register_layout = (
    cmd               = BitRange{RW}(30:31, DataTransferCmd), # Cmd
    mem_space_sel     = BitRange{RW}(28:29, MemSpaceSel), # Space select bits
    mem_addr          = BitRange{RW}( 0:27), # Memory 32-bit start address (128 Meg x 16 : only address bits 25-0 are used )
)
export data_transfer_ctl_register_layout

const mem_data_region = 0x00100000:0x00500000
export mem_data_region

const mem_fifo_region = first(mem_data_region):first(mem_data_region) + 0x000FFFFD
export mem_fifo_region

const fifo_start_addr = IdDict(
    1 => 0x00100000,
    2 => 0x00200000,
    3 => 0x00300000,
    4 => 0x00400000,
)
export fifo_start_addr

const fpga_ch_offset = IdDict(
    1 => 0x00000000,
    2 => 0x02000000,
    3 => 0x00000000,
    4 => 0x02000000,
)
export fpga_ch_offset

const bank_offset = IdDict(
    1 => 0x00000000,
    2 => 0x01000000,
)
export bank_offset

function fpga_ch_fifo_addr_offset(fpga_ch::Integer, bank::Integer)
    1 <= fpga_ch <= 4 || throw(ArgumentError("ADC FPGA channel number must be between 1 and 4"))
    1 <= bank <= 2 || throw(ArgumentError("Bank must be 1 or 2"))

    fpga_ch_offset[fpga_ch] + bank_offset[bank]
end
export fpga_ch_fifo_addr_offset

const fpga_ch_mem_space_sel = IdDict(
    1 => MemSpaceMem1,
    2 => MemSpaceMem1,
    3 => MemSpaceMem2,
    4 => MemSpaceMem2
)
export fpga_ch_mem_space_sel

function jumbo_frames_enabled(mem::SIS3316Memory)
    # Called during fifo I/O, bypass R/W task:
    r = UInt32[0]
    SIS3316Digitizers.read_registers!(mem.gw, [getaddress(SIS3316Digitizers.udp_protocol_config_reg)], r)
    getval(first(r), getlayout(SIS3316Digitizers.udp_protocol_config_reg.jumbo_packet_enable))
end
export jumbo_frames_enabled

function jumbo_frames_enabled!(mem::SIS3316Memory, enabled::Bool)
    mem[SIS3316Digitizers.udp_protocol_config_reg.jumbo_packet_enable] = enabled
end
export jumbo_frames_enabled!

function start_fifo_read!(mem::SIS3316Memory, ch::Integer, bank::Integer, from::Unsigned)
    @debug "Starting FIFO data transfer for channel $ch, bank $bank, address $(repr(from))"
    @argcheck from % 4 == 0
    @argcheck from < 0x10000000

    from_word = div(from, 4) #!!!!????

    group = fpga_num(ch)
    grp_ch = fpga_ch(ch)

    cmd = DataTransferCtlRead
    mem_space_sel = fpga_ch_mem_space_sel[grp_ch]
    mem_addr = fpga_ch_fifo_addr_offset(grp_ch, bank) + from_word

    ctlreg_addr = getaddress(data_transfer_ctrl_reg[group])
    ctlreg_val = data_transfer_ctrl_regval(cmd, mem_space_sel, mem_addr)

    # @debug "Writing value $(repr(ctlreg_val)) to data transfer control register $(repr(ctlreg_addr))"
    #mem.sync() #!!!
    #mem[ctlreg_addr] = ctlreg_val
    # Bypass R/W task:
    write_registers!(mem.gw, [ctlreg_addr], [ctlreg_val])
    nothing
end
export start_fifo_read!

function reset_fifo!(mem::SIS3316Memory, ch::Integer)
    @argcheck 1 <= ch <= 16
    @debug "Resetting FIFO data transfer for channel $ch"
    group = fpga_num(ch)

    ctlreg_addr = getaddress(data_transfer_ctrl_reg[group])
    ctlreg_val = setval(zero(RegVal), data_transfer_ctl_register_layout.cmd, DataTransferCtlReset)

    # @debug "Writing value $(repr(ctlreg_val)) to data transfer control register $(repr(ctlreg_addr))"
    #mem.sync() #!!!
    #mem[ctlreg_addr] = ctlreg_val
    # Bypass R/W task:
    write_registers!(mem.gw, [ctlreg_addr], [ctlreg_val])
    nothing
end
export reset_fifo!

const max_udp_req_pkg_size = 1485
const max_udp_resp_pkg_size = 1485
const max_udp_jumbo_pkg_size = 8237

const udp_resp_header_size = 45 # According to SIS3316 docs - why not 28 (17 + 8 + 3)?
const max_n_resp_data_frames = 15

const max_n_bytes_per_req = div(min(
    max_n_resp_data_frames * (max_udp_resp_pkg_size - udp_resp_header_size),
    (typemax(UInt16) + 1) * sizeof(UInt32)
), sizeof(UInt32)) * sizeof(UInt32)

function read_adc_fifo_raw!(mem::SIS3316Memory, address::MemAddr, data::AbstractVector{T}) where {T<:Unsigned}
    use_jumbo_frames = jumbo_frames_enabled(mem) # Register read also ensures no pending reads/writes

    n_bytes = Int(length(eachindex(data)) * sizeof(T))
    @argcheck n_bytes % sizeof(UInt32) == 0
    @debug "Reading $n_bytes from FIFO memory, starting at address $(repr(address))"

    udp_resp_size_limit = use_jumbo_frames ? max_udp_jumbo_pkg_size : max_udp_resp_pkg_size

    max_n_bytes_per_req = div(min(
        max_n_resp_data_frames * (udp_resp_size_limit - udp_resp_header_size),
        (typemax(UInt16) + 1) * sizeof(UInt32)
    ), sizeof(UInt32)) * sizeof(UInt32)

    data_words = reinterpret(UInt32, data)
    max_n_words_per_req = div(max_n_bytes_per_req, sizeof(UInt32))

    for data_bytes_part in Base.Iterators.partition(data_words, max_n_words_per_req)
        # With the SIS3316 UDP Interface, address does not increase from
        # chunk to chunk for a FIFO read (different from SIS3153):
        from_addr = address

        read_bulk!(mem.gw, from_addr, data_bytes_part)
    end

    data
end
export read_adc_fifo_raw!

function _read_ch_bank_data_padded!(data::AbstractVector{UInt8}, mem::SIS3316Memory, ch::Int, bank::Int, padded_from::MemAddr)
    @argcheck padded_from + length(eachindex(data)) <= MemAddr(0x04000000)

    lock(mem.fifo_lock)
    try
        fifo_addr = fifo_start_addr[fpga_num(ch) ]
        reset_fifo!(mem, ch)
        start_fifo_read!(mem, ch, bank, MemAddr(padded_from))
        read_adc_fifo_raw!(mem, fifo_addr, data)
        reset_fifo!(mem, ch)
    finally
        unlock(mem.fifo_lock)
    end

    nothing
end

function read_ch_bank_data(mem::SIS3316Memory, ch::Integer, bank::Integer, from::Integer = 0, n_bytes::Integer = ch_buffer_size)
    @argcheck from >= 0
    @argcheck n_bytes >= 0

    unpadded_from = MemAddr(from)
    unpadded_until = MemAddr(unpadded_from + n_bytes)
    unpadded_nbytes = unpadded_until - unpadded_from

    @debug "Reading bank $bank data for channel $ch from $(repr(unpadded_from)) until $(repr(unpadded_until))"

    padded_from = MemAddr(div(unpadded_from, 8) * 8) # Must be a multiple of 64 bit
    padded_until = MemAddr(div(unpadded_until + 7, 8) * 8) # Must be a multiple of 64 bit
    padded_nbytes = padded_until - padded_from

    @assert padded_from <= unpadded_from
    @assert padded_until >= unpadded_until

    data = Vector{UInt8}(undef, padded_nbytes)
    multiple_tries(5, (EOFError,), "Reading padded bank data for channel $ch") do
        _read_ch_bank_data_padded!(data, mem, Int(ch), Int(bank), padded_from)
    end

    view_start = firstindex(data) + (unpadded_from - padded_from)
    view_end = view_start + unpadded_nbytes - 1
    view(data, view_start:view_end)
end
export read_ch_bank_data


text:
@testset "bayesianeva.jl" begin

    n = 100

    x = Variable("x",randn(n))

    θ = [1; 1; -.5; .3; .1]

    μ = θ[1] .+ θ[2]*x.value
    ϕ = θ[3] .+ θ[4]*x.value
    σ = exp.(ϕ)
    ξ = θ[5]

    pd = GeneralizedExtremeValue.(μ,σ, ξ)

    y = rand.(pd)

    model = BlockMaxima(Variable("y", y), locationcov = [x] ,logscalecov = [x])

    fm = Extremes.BayesianEVA(model, Mamba.Chains(collect(θ')))

    rl = returnlevel(fm, 100)

    ci = cint(rl)

    @testset "getdistribution(fittedmodel)" begin
        @test all(vec(Extremes.getdistribution(fm)) .== pd)
    end

    @testset "quantile(fm, p)" begin
        # p not in [0, 1] throws
        @test_throws AssertionError Extremes.quantile(fm, -1)

        # Test with known values
        @test vec(quantile(fm, .95)) ≈ quantile.(pd, .95)
    end

    @testset "returnlevel(fm, returnPeriod)" begin
        # returnPeriod < 0 throws
        @test_throws AssertionError returnlevel(fm, -100)

        # Test with known values
        @test vec(rl.value) ≈ quantile.(pd, 1-1/100)
    end

    @testset "cint(fm, returnPeriod, confidencelevel)" begin
        # confidencelevel not in [0, 1]
        @test_throws AssertionError cint(rl, -1)

        # Test with known values
        @test ci[1][1] ≈ rl.value[1,1]
        @test ci[5][1] ≈ rl.value[1,5]

    end

    n = 1000
    nobservation = 1000
    nobsperblock = 1

    threshold = 10.0

    x = Variable("x",randn(n))

    θ = [-.5; 1; 0.1]

    ϕ = θ[1] .+ θ[2]*x.value
    σ = exp.(ϕ)
    ξ = θ[3]

    pd = GeneralizedPareto.(σ, ξ)

    y = rand.(pd)

    model = ThresholdExceedance(Variable("y", y), logscalecov = [x])

    fm = Extremes.BayesianEVA(model, Mamba.Chains(collect(θ')))

    rl = returnlevel(fm, 0, nobservation, nobsperblock, 100)

    ci = cint(rl)

    @testset "getdistribution(fittedmodel)" begin
        @test all(vec(Extremes.getdistribution(fm)) .== pd)
    end

    @testset "quantile(fm, p)" begin
        # p not in [0, 1] throws
        @test_throws AssertionError Extremes.quantile(fm, -1)

        # Test with known values
        @test vec(quantile(fm, .95)) ≈ quantile.(pd, .95)
    end

    @testset "returnlevel(fm, returnPeriod)" begin
        # returnPeriod < 0 throws
        @test_throws AssertionError returnlevel(fm, 0, nobservation, nobsperblock, -100)

        # Test with known values
        p = 1-nobservation/(100 * nobsperblock * n)
        @test vec(rl.value) ≈ quantile.(pd, p)
    end

    @testset "cint(fm, returnPeriod, confidencelevel)" begin
        # confidencelevel not in [0, 1]
        @test_throws AssertionError cint(returnlevel(fm, 0, nobservation, nobsperblock, 100), -1)

        # Test with known values
        @test ci[1][1] ≈ rl.value[1,1]
        @test ci[5][1] ≈ rl.value[1,5]

    end

    @testset "showfittedEVA(io, obj, prefix)" begin
        # print does not throw
        buffer = IOBuffer()
        @test_logs Extremes.showfittedEVA(buffer, fm, prefix = "\t")
    end

    @testset "showChain(io, chain, prefix)" begin
        # print does not throw
        buffer = IOBuffer()
        @test_logs Extremes.showChain(buffer, fm.sim)
    end

    @testset "findposteriormode(fm::BayesianEVA)" begin

        x = Variable("x", randn(10))
        μ = 10 .+ x.value
        σ = 1.0
        ξ = .1
        pd = GeneralizedExtremeValue.(μ, σ, ξ)
        y = rand.(pd)
        fm = Extremes.BayesianEVA(Extremes.BlockMaxima(Variable("y", y), locationcov=[x]),
            Mamba.Chains([10.0 1.0 0.0 .1; -10.0 1.0 0.0 .1; 20.0 1.0 0.0 .1]))

        θ̂ = Extremes.findposteriormode(fm)
        @test θ̂ ≈ [10.0; 1.0; 0.0; .1]

    end

end


text:
using Test
using ProbNumDiffEq
using TimerOutputs
const to = TimerOutput()
macro timedtestset(name, code)
    return esc(:(println("Start testset: ", $name);
    @timeit to $name @testset $name $code;
    println("Done.")))
end

@testset "ProbNumDiffEq" begin
    @timedtestset "Correctness" begin
        include("correctness.jl")
    end

    @timedtestset "Filtering" begin
        include("filtering.jl")
    end

    @timedtestset "Convergence" begin
        include("convergence.jl")
    end

    @timedtestset "Priors" begin
        include("priors.jl")
    end

    @timedtestset "Solution" begin
        include("solution.jl")
    end

    @timedtestset "Diffusions" begin
        include("diffusions.jl")
    end

    @timedtestset "State Initialization" begin
        include("state_init.jl")
    end

    @timedtestset "Preconditioning" begin
        include("preconditioning.jl")
    end

    @timedtestset "Smoothing" begin
        include("smoothing.jl")
    end

    @timedtestset "Errors" begin
        include("errors.jl")
    end

    @timedtestset "IEKS" begin
        include("ieks.jl")
    end

    @timedtestset "Automatic Differentiation" begin
        include("autodiff.jl")
    end

    @timedtestset "Second order ODEs" begin
        include("secondorderode.jl")
    end

    @timedtestset "Specific Problems" begin
        include("specific_problems.jl")
    end

    @timedtestset "DiffEqDevTools.jl Compatibility" begin
        include("diffeqdevtools.jl")
    end
end

println("\nTimer Summary:")
display(to)
println("\n")


text:
asDouble(x) = float(x)
asDouble(x :: Complex) = asDouble(real(x))
asDouble(x :: AbstractArray{T}) where {T <: Complex} = [asDouble(xx) for xx in x]

asCharacter(x) = string(x)

asLogical(x :: Bool) = x
asLogical(x :: Number) = x != 0
asLogical(x :: AbstractString) = ifelse(x in ["T", "TRUE", "true", "True"], true, false)
asLogical(x :: AbstractArray{Bool}) = x
asLogical(x :: AbstractArray{T}) where {T} = [asLogical(xx) for xx in x]

isNumeric(x) = isa(x, Real)
isNumeric(x :: AbstractArray{T}) where {T} = T <: Real

max_stars_repo_name: KevCaz/JuliaCall

max_issues_repo_name: tinagui/JuliaCall

max_forks_repo_name: tinagui/JuliaCall


text:
"""
    IndexedBiDiGraph{T<:Integer} <: AbstractIndexedDiGraph{T}

A type representing a sparse directed graph with access to both outedges and inedges.

### FIELDS

- `A` -- square matrix filled with `NullNumber`s. `A[i,j]` corresponds to edge `j=>i`.
- `X` -- square matrix for efficient access by row. `X[j,i]` points to the index of element `A[i,j]` in `A.nzval`.
"""
struct IndexedBiDiGraph{T<:Integer} <: AbstractIndexedDiGraph{T}
    A :: SparseMatrixCSC{NullNumber, T}
    X :: SparseMatrixCSC{T, T}
end

"""
    IndexedBiDiGraph(A::AbstractMatrix)

Construct an `IndexedBiDiGraph` from the adjacency matrix `A`.

`IndexedBiDiGraph` internally stores the transpose of `A`. To avoid overhead due
to the transposition, use `IndexedBiDiGraph(transpose(At))` where `At` is the
transpose of `A`.
"""
function IndexedBiDiGraph(A::AbstractMatrix)
    _checksquare(A)
    _check_selfloops(A)
    B = convert(SparseMatrixCSC, transpose(A))
    X = permutedims( SparseMatrixCSC(B.m, B.n, B.colptr, B.rowval, collect(1:nnz(B))) )
    At = SparseMatrixCSC(B.m, B.n, B.colptr, B.rowval, fill(NullNumber(), nnz(B)))
    IndexedBiDiGraph(At, X)
end

"""
    IndexedBiDiGraph(A::AbstractSimpleGraph)

Construct an `IndexedBiDiGraph` from any `AbstractSimpleGraph` (Graphs.jl),
directed or otherwise.
"""
IndexedBiDiGraph(sg::AbstractSimpleGraph) = IndexedBiDiGraph(adjacency_matrix(sg))

inneighbors(g::IndexedBiDiGraph, i::Integer) = @inbounds @view g.X.rowval[nzrange(g.X,i)]

Base.zero(g::IndexedBiDiGraph) = IndexedBiDiGraph(zero(g.A))

"""
    inedges(g::AbstractIndexedBiDiGraph, i::Integer)

Return a lazy iterator to the edges ingoing to node `i` in `g`.
"""
inedges(g::IndexedBiDiGraph, i::Integer) = @inbounds (IndexedEdge{Int}(g.X.rowval[k], i, g.X.nzval[k]) for k in nzrange(g.X, i))

# Return a copy of the adjacency matrix with elements of type `T`
function adjacency_matrix(g::IndexedBiDiGraph, T::DataType=Int)
    SparseMatrixCSC(g.X.m, g.X.n, g.X.colptr, g.X.rowval, ones(T, nnz(g.X)))
end


text:
using OpenQuantumBase, Test

struct T_OPENSYS <: OpenQuantumBase.AbstractLiouvillian end
struct T_COUPLINGS <: OpenQuantumBase.AbstractCouplings end
struct T_BATH <: OpenQuantumBase.AbstractBath end

H = DenseHamiltonian([(x) -> x], [σz])
u0 = PauliVec[1][1]

annealing = Annealing(H, u0)
@test annealing.H == H
@test annealing.annealing_parameter(10, 5) == 0.5

evo = Evolution(H, u0)
@test evo.H == H
@test evo.annealing_parameter(10, 5) == 0.5

ode_params = ODEParams(T_OPENSYS(), 10, (tf, t)->t / tf)
@test typeof(ode_params.L) == T_OPENSYS
@test ode_params.tf == 10
@test ode_params(5) == 0.5

coupling = ConstantCouplings(["Z"])
inter = Interaction(coupling, T_BATH())
inter_set = InteractionSet(inter, inter)
@test inter_set[1] == inter
annealing = Annealing(H, u0, coupling = coupling, bath = T_BATH())
@test annealing.interactions[1].coupling == coupling
@test typeof(annealing.interactions[1].bath) <: T_BATH
@test_throws ArgumentError Annealing(H, u0, coupling = coupling, bath = T_BATH(), interactions = inter_set)


text:
# API requests data structures

# report request

@enum SamplingLevel begin
    SAMPLING_UNSPECIFIED
    DEFAULT
    SMALL
    LARGE
end

# generics

@enum Operator begin
    OPERATOR_UNSPECIFIED
end

@enum GeneralOperator <: Operator begin
    REGEXP
    BEGINS_WITH
    ENDS_WITH
    PARTIAL
    EXACT
    NUMERIC_EQUAL
    NUMERIC_GREATER_THAN
    NUMERIC_LESS_THAN
    IN_LIST
end

@enum DimensionOperator <: Operator begin
    OR
    AND
end

@enum MetricOperator <: Operator begin
    EQUAL = 2
    LESS_THAN = 3
    GREATER_THAN = 4
    IS_MISSING = 5
end

# general dimension struct
struct Dimension
    name::AbstractString
    histogram_buckets::Vector{Int64}
end

struct DimensionFilter
    dimension_name::AbstractString
    not::Bool
    operator::GeneralOperator
    expressions::Vector{AbstractString}
    case_sensitive::Bool
end

struct DimensionFilter
    operator::DimensionOperator
    filters::Vector{DimensionFilter}
end

# Metrics related

@enum MetricType begin
    METRIC_TYPE_UNSPECIFIED
    INTEGER
    FLOAT
    CURRENCY
    PERCENT
    TIME
end

struct Metric
    expression::AbstractString
    alias::AbstractString
    formatting_type::MetricType
end

struct MetricFilter
    metric_name::AbstractString
    not::Bool
    operator::MetricOperator
end

@enum OrderType begin
    ORDER_TYPE_UNSPECIFIED
    VALUE
    DELTA
    SMART
    HISTOGRAM_BUCKET
    DIMENSION_AS_INTEGER
end

struct OrderBy
    field_name::AbstractString
    order_type::OrderType
    sort_order::SortOrder
end

@enum SortOrder begin
end

struct Segment end

struct DynamicSegment end

struct SegmentDefinition end

struct SegmentFilter end

struct SimpleSegment end

struct SequenceSegment end

struct OrFiltersForSegment end

struct SegmentFilterClause end

struct SegmentDimensionFilter end

struct SegmentMetricFilterClause end

struct SegmentSequenceStep end

@enum Scope begin
end

@enum MatchType begin
end

struct Pivot end

# cohorts
struct CohortGroup end

struct Cohort end

@enum CohortType begin
end

struct ReportRequest
    viewid::AbstractString
    dateranges::DateRange
    sampling_level::SamplingLevel
    dimensions::Vector{Dimension}
    dimension_filters::DimensionFilter
    metrics::Vector{Metric}
    metric_filters::MetricFilter
    filter::AbstractString
    order_by::OrderBy
    segments::Segment
    pivots::Pivot
    cohort_group::CohortGroup
    page_token::AbstractString
    page_size::Integer
    include_null_rows::Bool
    hide_totals::Bool
    hide_value_ranges::Bool
end

struct Request
    report_requests::Vector{ReportRequest}
    use_resource_quotas::Bool
end