text:
import Base: +, getindex, isempty

export MinkowskiSum, ⊕,
       MinkowskiSum!,
       swap

"""
    MinkowskiSum{N<:Real, S1<:LazySet{N}, S2<:LazySet{N}} <: LazySet{N}

Type that represents the Minkowski sum of two convex sets.

### Fields

- `X` -- first convex set
- `Y` -- second convex set

### Notes

The `ZeroSet` is the neutral element and the `EmptySet` is the absorbing element
for `MinkowskiSum`.
"""
struct MinkowskiSum{N<:Real, S1<:LazySet{N}, S2<:LazySet{N}} <: LazySet{N}
    X::S1
    Y::S2

    # default constructor with dimension check
    function MinkowskiSum(X::S1, Y::S2) where {N<:Real, S1<:LazySet{N},
                                               S2<:LazySet{N}}
        @assert dim(X) == dim(Y) "sets in a Minkowski sum must have the " *
            "same dimension"
        return new{N, S1, S2}(X, Y)
    end
end

isoperationtype(::Type{<:MinkowskiSum}) = true
isconvextype(::Type{MinkowskiSum{N, S1, S2}}) where {N, S1, S2} = isconvextype(S1) && isconvextype(S2)

# ZeroSet is the neutral element for MinkowskiSum
@neutral(MinkowskiSum, ZeroSet)

# EmptySet and Universe are the absorbing elements for MinkowskiSum
@absorbing(MinkowskiSum, EmptySet)
# @absorbing(MinkowskiSum, Universe)  # TODO problematic

"""
    X + Y

Convenience constructor for Minkowski sum.

### Input

- `X` -- a convex set
- `Y` -- another convex set

### Output

The symbolic Minkowski sum of ``X`` and ``Y``.
"""
+(X::LazySet, Y::LazySet) = MinkowskiSum(X, Y)

"""
    ⊕(X::LazySet, Y::LazySet)

Unicode alias constructor ⊕ (`oplus`) for the lazy Minkowski sum operator.
"""
⊕(X::LazySet, Y::LazySet) = MinkowskiSum(X, Y)

"""
    swap(ms::MinkowskiSum)

Return a new `MinkowskiSum` object with the arguments swapped.

### Input

- `ms` -- Minkowski sum of two convex sets

### Output

A new `MinkowskiSum` object with the arguments swapped.
"""
function swap(ms::MinkowskiSum)
    return MinkowskiSum(ms.Y, ms.X)
end

"""
    dim(ms::MinkowskiSum)

Return the dimension of a Minkowski sum.

### Input

- `ms` -- Minkowski sum

### Output

The ambient dimension of the Minkowski sum.
"""
function dim(ms::MinkowskiSum)
    return dim(ms.X)
end

"""
    σ(d::AbstractVector{N}, ms::MinkowskiSum{N}) where {N<:Real}

Return the support vector of a Minkowski sum.

### Input

- `d`  -- direction
- `ms` -- Minkowski sum

### Output

The support vector in the given direction.
If the direction has norm zero, the result depends on the summand sets.

### Algorithm

The support vector in direction ``d`` of the Minkowski sum of two sets ``X``
and ``Y`` is the sum of the support vectors of ``X`` and ``Y`` in direction
``d``.
"""
function σ(d::AbstractVector{N}, ms::MinkowskiSum{N}) where {N<:Real}
    return σ(d, ms.X) + σ(d, ms.Y)
end

"""
    ρ(d::AbstractVector{N}, ms::MinkowskiSum{N}) where {N<:Real}

Return the support function of a Minkowski sum.

### Input

- `d`  -- direction
- `ms` -- Minkowski sum

### Output

The support function in the given direction.

### Algorithm

The support function in direction ``d`` of the Minkowski sum of two sets ``X``
and ``Y`` is the sum of the support functions of ``X`` and ``Y`` in direction
``d``.
"""
function ρ(d::AbstractVector{N}, ms::MinkowskiSum{N}) where {N<:Real}
    return ρ(d, ms.X) + ρ(d, ms.Y)
end

"""
    isbounded(ms::MinkowskiSum)

Determine whether a Minkowski sum is bounded.

### Input

- `ms` -- Minkowski sum

### Output

`true` iff both wrapped sets are bounded.
"""
function isbounded(ms::MinkowskiSum)
    return isbounded(ms.X) && isbounded(ms.Y)
end

"""
    isempty(ms::MinkowskiSum)

Return if a Minkowski sum is empty or not.

### Input

- `ms` -- Minkowski sum

### Output

`true` iff any of the wrapped sets are empty.
"""
function isempty(ms::MinkowskiSum)
    return isempty(ms.X) || isempty(ms.Y)
end

"""
    constraints_list(ms::MinkowskiSum)

Return the list of constraints of a lazy Minkowski sum of two polyhedral sets.

### Input

- `ms` -- Minkowski sum of two polyhedral sets

### Output

The list of constraints of the Minkowski sum.

### Algorithm

We compute a concrete set representation via `minkowski_sum` and call
`constraints_list` on the result.
"""
function constraints_list(ms::MinkowskiSum)
    return constraints_list(minkowski_sum(ms.X, ms.Y))
end

"""
    ∈(x::AbstractVector{N}, ms::MinkowskiSum{N, <:AbstractSingleton, <:LazySet}) where {N}

Check whether a given point is contained in the Minkowski sum of a singleton
and a set.

### Input

- `x`  -- point
- `ms` -- lazy Minkowski sum of a singleton and a set

### Output

`true` iff ``x ∈ ms``.

### Algorithm

Note that ``x ∈ (S ⊕ P)``, where ``S`` is a singleton set, ``S = \\{s\\}`` and
``P`` is a set, if and only if ``(x-s) ∈ P``.
"""
function ∈(x::AbstractVector{N}, ms::MinkowskiSum{N, S1, S2}) where {N, S1<:AbstractSingleton, S2<:LazySet}
    return _in_singleton_msum(x, ms.X, ms.Y)
end

# symmetric method
function ∈(x::AbstractVector{N}, ms::MinkowskiSum{N, <:LazySet, <:AbstractSingleton}) where {N}
    return _in_singleton_msum(x, ms.Y, ms.X)
end

# disambiguation
function ∈(x::AbstractVector{N}, ms::MinkowskiSum{N, <:AbstractSingleton, <:AbstractSingleton}) where {N}
    return _in_singleton_msum(x, ms.X, ms.Y)
end

@inline _in_singleton_msum(x, X, Y) = (x - element(X)) ∈ Y

function concretize(ms::MinkowskiSum)
    return minkowski_sum(concretize(ms.X), concretize(ms.Y))
end

# ================
# Helper functions
# ================

@inline function σ_helper(d::AbstractVector{N},
                          array::AbstractVector{<:LazySet}) where {N<:Real}
    svec = zeros(N, length(d))
    for sj in array
        svec += σ(d, sj)
    end
    return svec
end

"""
    vertices_list(ms::MinkowskiSum{N, Z1, Z2}) where {N<:Real, Z1<:AbstractZonotope{N}, Z2<:AbstractZonotope{N}}

Return the list of vertices for the Minkowski sum of two zonotopic sets.

### Input

- `ms` -- Minkowski sum of two zonotopic sets

### Output

The list of vertices of the Minkowski sum of two zonotopic sets.
"""
function vertices_list(ms::MinkowskiSum{N, Z1, Z2}) where {N<:Real, Z1<:AbstractZonotope{N}, Z2<:AbstractZonotope{N}}
    return vertices_list(minkowski_sum(ms.X, ms.Y))
end


text:
"""
        elastic_forces(rod::basic_rod, rod_props::rod_properties, ref_strains::rod_strains)

Calculate elastic forces acting on a rod. Output is a rod_delta type, containing
the (positive) gradient of elastic energy with respect to each degree of freedom.
"""
function elastic_forces(rod::basic_rod, rod_props::rod_properties, ref_strains::rod_strains)

    # unpack
    l0, κ0, τ0 =  ref_strains.l, ref_strains.κ, ref_strains.τ
    k, B, β = rod_props.k, rod_props.B, rod_props.β
    x, d1 = rod.x, rod.d

    # Tangents and directors
    T = eltype(x)
    nv = size(x,1)
    e = x[2:nv,:] .- x[1:nv-1,:]
    l = norm3(e)
    t = e./l
    d2 = cross3(t, d1)

    # Split into halves of the rod
    l1 = @view l[1:nv-2]
    l2 = @view l[2:nv-1]
    t1 = @view t[1:nv-2,:]
    t2 = @view t[2:nv-1,:]
    d1l = @view d1[1:nv-2,:]
    d1r = @view d1[2:nv-1,:]
    d2l = @view d2[1:nv-2,:]
    d2r = @view d2[2:nv-1,:]

    chi = 1 .+ dot3(t1, t2)
    ttild = (t1 .+ t2)./chi
    d1tild = (d1l .+ d1r)./chi
    d2tild = (d2l .+ d2r)./chi
    kb = 2 .*cross3(t1, t2)./chi

    κ1l =  dot3(kb, d2l)
    κ1r =  dot3(kb, d2r)
    κ2l = -dot3(kb, d1l)
    κ2r = -dot3(kb, d1r)
    κ1 = (κ1l + κ1r)./2
    κ2 = (κ2l + κ2r)./2
    κ = hcat(κ1, κ2)

    # Curvature Derivatives
    dκ1de1 = (-κ1.*ttild .+ cross3(t2, d2tild))./l1
    dκ1de2 = (-κ1.*ttild .- cross3(t1, d2tild))./l2

    dκ2de1 = (-κ2.*ttild .- cross3(t2, d1tild))./l1
    dκ2de2 = (-κ2.*ttild .+ cross3(t1, d1tild))./l2

    # Twists
    d1_transport = ptransport(d1l, t1, t2)
    cosτ = dot3(d1_transport, d1r)
    sinτ = dot3(d1_transport, d2r)
    τ = -atan.(sinτ, cosτ)

    # Twist Derivatives
    dτde1 = kb./(2 .*l1);
    dτde2 = kb./(2 .*l2);

    # Energy derivatives
    vl0 = [l0[1]; l0[2:end-1]/2] + [l0[2:end-1]/2; l0[end]]
    dEdl = k.*(l - l0)./l0
    dEdκ = B.*(κ .- κ0)./vl0
    dEdτ = β.*(τ .- τ0)./vl0

    dEsde = dEdl.*t
    dEde1 = dEdκ[:,1].*dκ1de1 .+ dEdκ[:,2].*dκ2de1 .+ dEdτ.*dτde1
    dEde2 = dEdκ[:,1].*dκ1de2 .+ dEdκ[:,2].*dκ2de2 .+ dEdτ.*dτde2

    dEdx = [-dEde1; zeros(T,2,3)] .+ [zeros(T,1,3); dEde1 .- dEde2; zeros(T,1,3)] .+ [zeros(T,2,3); dEde2]
    dEdx += [-dEsde; zeros(T,1,3)] .+ [zeros(T,1,3); dEsde]

    dEtdθ = [0; dEdτ] .- [dEdτ; 0]
    dEbdθ = ([(κ2l.*dEdκ[:,1] .- κ1l.*dEdκ[:,2]); 0]./2
          .+ [0; (κ2r.*dEdκ[:,1] .- κ1r.*dEdκ[:,2])]./2)
    dEdθ = dEtdθ .+ dEbdθ

    rod_delta(dEdx, dEdθ)
end

function gradient_update(x0, d0, Δx, Δθ)

    x1 = x0 + Δx
    e1 = edges(x1)
    l1 = norm3(e1)
    t1 = e1./l1
    t0 = tangents(x0)
    d1 = ptransport(d0, t0, t1)
    d1 = rotate_orthogonal_unit(d1, t1, Δθ)

    ζ = cross3(t1, t0)./(l1.*(1 .+ dot3(t1, t0)))

    x1, d1, ζ
end

function gradient_update(ref::basic_rod, dr::rod_delta)
    x, d, ζ = gradient_update(ref.x, ref.d, dr.Δx, dr.Δθ)
    basic_rod(x, d), ζ
end

# Cache should be ((ns + 1 x 7), (ns x 28))
function elastic_forces!(forces::rod_delta, rod::basic_rod,
    rod_props::rod_properties, ref_strains::rod_strains, caches)

    # unpack
    l0, κ0, τ0 =  ref_strains.l, ref_strains.κ, ref_strains.τ
    k, B, β = rod_props.k, rod_props.B, rod_props.β
    x, d1 = rod.x, rod.d
    dEdx, dEdθ = forces.Δx, forces.Δθ

    elastic_forces!(dEdx, dEdθ, x, d1, k, B, β, l0, κ0, τ0, caches)
end

function elastic_forces!(dEdx, dEdθ, x, d1, k, B, β, l0, κ0, τ0, caches)

    nv = size(x,1)
    ne = nv - 1

    # Unpack and divy up cache1
    cache1, cache2 = caches
    t = view(cache1, :, 1:3)
    l = view(cache1, :, 4)
    dEsde = view(cache1, :, 5:7) # Overwritten early
    d2 = view(cache1, :, 5:7)

    # Permanent cache2 space (never overwritten)
    vl0 = view(cache2, :, 1)
    kb = view(cache2, :, 2:4)
    dEde1 = view(cache2, :, 5:7)
    dEde2 = view(cache2, :, 8:10)

    # Bending cache2 space (overwritten for twist)
    chi = view(cache2, :, 11)
    ttild = view(cache2, :, 12:14)
    dtild = view(cache2, :, 15:17)
    κ1l = view(cache2, :, 18)
    κ1r = view(cache2, :, 19)
    κ2l = view(cache2, :, 20)
    κ2r = view(cache2, :, 21)
    κ1 = view(cache2, :, 22)
    κ2 = view(cache2, :, 23)
    dEdκ1 = view(cache2, :, 24)
    dEdκ2 = view(cache2, :, 25)
    dκde = view(cache2, :, 26:28)

    # Edges and lengths
    edges!(t, x)
    norm3!(l, t)
    t ./= l

    vl0[1] = l0[1] + l0[2]/2
    vl0[2:ne-2] .= (view(l0, 2:ne-2) .+ view(l0, 3:ne-1))./2
    vl0[ne-1] = l0[ne-1]/2 + l0[ne]

    # Stretching Forces
    dEsde .= (k.*(l .- l0)./l0).*t
    dEdx .= 0.0
    view(dEdx, 1:nv-1, :) .-= dEsde
    view(dEdx, 2:nv, :) .+= dEsde

    # Break tangents and directors into sets
    cross3!(d2, t, d1)
    l1 = @view l[1:nv-2]
    l2 = @view l[2:nv-1]
    t1 = @view t[1:nv-2,:]
    t2 = @view t[2:nv-1,:]
    d1l = @view d1[1:nv-2,:]
    d1r = @view d1[2:nv-1,:]
    d2l = @view d2[1:nv-2,:]
    d2r = @view d2[2:nv-1,:]

    # Curvature
    dot3!(chi, t1, t2)
    chi .+= 1.0
    cross3!(kb, t1, t2)
    kb .*= 2.0./chi

    dot3!(κ1l, kb, d2l)
    dot3!(κ1r, kb, d2r)
    dot3!(κ2l, kb, d1l)
    dot3!(κ2r, kb, d1r)
    κ2l .*= -1.0
    κ2r .*= -1.0
    κ1 .= (κ1l .+ κ1r)./2
    κ2 .= (κ2l .+ κ2r)./2
    dEdκ1 .= view(B, :, 1).*(κ1 .- view(κ0, :, 1))./vl0
    dEdκ2 .= view(B, :, 2).*(κ2 .- view(κ0, :, 2))./vl0

    # Curvature Derivatives
    ttild .= (t1 .+ t2)./chi
    dtild .= (d2l .+ d2r)./chi

    # dκ1de1
    dtild .= (d2l .+ d2r)./chi # d2tild
    cross3!(dκde, t2, dtild)
    dκde .-= κ1.*ttild
    dκde ./= l1
    dEde1 .= dEdκ1.*dκde

    # dκ1de2
    cross3!(dκde, t1, dtild)
    dκde .+= κ1.*ttild
    dκde ./= -1.0.*l2
    dEde2 .= dEdκ1.*dκde

    # dκ2de1
    dtild .= (d1l .+ d1r)./chi # d1tild
    cross3!(dκde, t2, dtild)
    dκde .+= κ2.*ttild
    dκde ./= -1.0.*l1
    dEde1 .+= dEdκ2.*dκde

    # dκ2de2
    cross3!(dκde, t1, dtild)
    dκde .-= κ2.*ttild
    dκde ./= l2
    dEde2 .+= dEdκ2.*dκde

    # Bending forces and moments
    dEdθ .= 0.0
    view(dEdθ, 1:nv-2) .+= (κ2l.*dEdκ1 .- κ1l.*dEdκ2)./2
    view(dEdθ, 2:nv-1) .+= (κ2r.*dEdκ1 .- κ1r.*dEdκ2)./2

    # Reallocate Cache Space for Twists
    pt_cache = view(cache2, :, 11:16)
    d1trans = view(cache2, :, 17:19)
    cosτ = view(cache2, :, 20)
    sinτ = view(cache2, :, 21)
    τ = view(cache2, :, 22)
    dEdτ = view(cache2, :, 23)

    # Twists
    ptransport!(d1trans, d1l, t1, t2, pt_cache)
    dot3!(cosτ, d1r, d1trans)
    dot3!(sinτ, d2r, d1trans)
    τ .= -1.0.*atan.(sinτ, cosτ)
    #
    # # Twisting Moments
    dEdτ .= β.*(τ .- view(τ0, :))./vl0
    view(dEdθ, 1:nv-2) .-= dEdτ
    view(dEdθ, 2:nv-1) .+= dEdτ
    dEde1 .+= dEdτ.*kb./(2 .*l1)
    dEde2 .+= dEdτ.*kb./(2 .*l2)

    # Assemble into output array
    view(dEdx, 1:nv-2, :) .-= dEde1
    view(dEdx, 2:nv-1, :) .+= dEde1 .- dEde2
    view(dEdx, 3:nv, :) .+= dEde2

    return nothing
end

allocate_cache(f::typeof(elastic_forces!), T::Type, ns) =
    (zeros(T, ns+1, 7), zeros(T, ns, 28))

function gradient_update!(x1, d1, ζ, x0, d0, Δx, Δθ, cache)
    t0 = view(cache, :, 1:3)
    t1 = view(cache, :, 4:6)
    l1 = view(cache, :, 7)
    t1_dot_t0 = view(cache, :, 8)
    cache_pt = view(cache, :, 8:13)
    cache_rt = view(cache, :, 8:12)

    tangents!(t0, x0)

    x1 .= x0 .+ Δx
    edges!(t1, x1)
    norm3!(l1, t1)
    t1 ./= l1

    ptransport!(d1, d0, t0, t1, cache_pt)
    rotate_orthogonal_unit!(d1, t1, Δθ, cache_pt)

    cross3!(ζ, t1, t0)
    dot3!(t1_dot_t0, t1, t0)
    ζ ./= (l1.*(1.0 .+ t1_dot_t0))
end

function gradient_update!(new::basic_rod, ζ, old::basic_rod, dr::rod_delta, cache)
    gradient_update!(new.x, new.d, ζ, old.x, old.d, dr.Δx, dr.Δθ, cache)
end

allocate_cache(f::typeof(gradient_update!), T::Type, ns) = zeros(T, ns+1, 13)

#--------------------------------------------------------
# Fusing the gradient update and forces step
#--------------------------------------------------------

function fused_update_gradient!(∇E, r, Δr, p, s, caches)

    cache_vertices, cache_edges, cache_interior = caches

    gu_cache = cache_edges
    ef_cache = (cache_edges, cache_interior)
    x = cache_vertices
    d = @view cache_edges[:, 14:16]
    ζ = @view cache_edges[:, 17:19]

    gradient_update!(x, d, ζ, r.x, r.d, Δr.Δx, Δr.Δθ, gu_cache)
    elastic_forces!(∇E.Δx, ∇E.Δθ, x, d, p.k, p.B, p.β, s.l, s.κ, s.τ, ef_cache)
    (@view ∇E.Δx[1:end-1, :]) .-= ∇E.Δθ .* ζ
    (@view ∇E.Δx[2:end,:]) .+= ∇E.Δθ .* ζ
end

allocate_cache(f::typeof(fused_update_gradient!), T::Type, ns) =
    (zeros(T, ns+2, 3), zeros(T, ns+1, 19), zeros(T, ns, 28))

ext: jl


text:
using DeIdentification
using CSV
using YAML
using DataFrames
using Dates
using Memento

using Test

# SET UP COMMON TEST VARIABLES AND ENVIRONMENT
test_file = "ehr_data.yml"
test_file_batch= "ehr_data_batch.yml"
logpath = joinpath(@__DIR__, "logs")
outputpath = joinpath(@__DIR__, "output")
basepath_not_created = joinpath(@__DIR__, "path")
logpath_not_created = joinpath(basepath_not_created, "to", "new", "logs")
outputpath_not_created = joinpath(basepath_not_created, "to", "new", "output")
default_max_days = 30
default_shift_years = 100
default_date_format = Dates.DateFormat("y-m-dTH:M:S.s")

try
    isdir(logpath) && rm(logpath, recursive=true, force=true)
catch
end

try
    isdir(outputpath) && rm(outputpath, recursive=true, force=true)
catch
end

try
    isdir(basepath_not_created) && rm(basepath_not_created, recursive=true, force=true)
catch
end

mkpath(logpath)
mkpath(outputpath)
# ----------------------------

@testset "config creation" begin
    cfg_raw = YAML.load_file(test_file)

    # nominally check YAML loading worked
    @test cfg_raw["project"] == "ehr"

    cfg = ProjectConfig(test_file)

    @test cfg_raw["project"] == cfg.name
    @test cfg_raw["datasets"][1]["name"] == cfg.file_configs[1].name
end

@testset "create config from csv" begin
    config_file = joinpath(@__DIR__, "example_config.yml")
    try
        build_config_from_csv("example_config", joinpath(@__DIR__, "data_worksheet_fields.csv"))

        config_raw = YAML.load_file(config_file)
        @test config_raw["project"] == "example_config"
        @test length(config_raw["datasets"]) == 11
    finally
        if isfile(config_file)
            try
                rm(config_file)
            catch
            end
        end
    end
end

@testset "rid generation" begin
    ids1 = [4, 6, 7, 3, 3, 5, 7]
    ids2 = [6, 5, 3, 4, 5]

    # Check hashing and research ID generation
    dicts = DeIdDicts(default_max_days, default_shift_years, default_date_format)
    hash1 = map( x-> DeIdentification.getoutput(dicts, DeIdentification.Hash, x, 0), ids1)
    hash2 = map( x-> DeIdentification.getoutput(dicts, DeIdentification.Hash, x, 0), ids2)

    rid1 = map( x-> DeIdentification.setrid(x, dicts), hash1)
    rid2 = map( x-> DeIdentification.setrid(x, dicts), hash2)

    @test rid1 == [1, 2, 3, 4, 4, 5, 3]
    @test rid2 == [2, 5, 4, 1, 5]
end

@testset "integration tests" begin
    proj_config = ProjectConfig(test_file)
    proj_config_batch = ProjectConfig(test_file_batch)
    deid = deidentify(proj_config)

    @test typeof(deid) == DeIdDicts

    @test isfile(joinpath(logpath,"ehr.log.0001"))

    dx = false
    salts = false
    df = DataFrame()
    df_pat = DataFrame()
    for (root, dirs, files) in walkdir(outputpath)
        for file in files
            if occursin(r"^deid_dx_.*csv", file)
                dx = true
                df = CSV.read(joinpath(root,file))
            elseif occursin(r"deid_pat_.*csv", file)
                df_pat = CSV.read(joinpath(root,file))
            elseif occursin(r"salts_.*json", file)
                salts = true
            end
        end
    end

    dfo = CSV.read(joinpath(@__DIR__, "data", "dx.csv"))

    # test column name change
    @test in(:EncounterBrownCSN, getfield(getfield(dfo, :colindex),:names))
    @test in(:CSN, getfield(getfield(df, :colindex),:names))

    # test that hash column was hashed
    @test length(df[1, :PatientPrimaryMRN]) == 64

    # test that dropped column was dropped
    @test in(:DiagnosisTerminologyType, getfield(getfield(dfo, :colindex),:names))
    @test !in(:DiagnosisTerminologyType, getfield(getfield(df, :colindex),:names))

    # test that dateshifted column was dateshifted
    @test df[1,:ArrivalDateandTime] != dfo[1,:ArrivalDateandTime]
    @test Dates.days(abs(df[1,:ArrivalDateandTime] - dfo[1,:ArrivalDateandTime])) <= proj_config.maxdays + Dates.days(Dates.Year(proj_config.shiftyears))
    @test isapprox(Dates.days(abs(df[1,:ArrivalDateandTime] - dfo[1,:ArrivalDateandTime]))/365, proj_config.shiftyears, atol = 1.0)

    # test the transforms work
    @test length(string(df_pat[1,:PatientBirthDate])) == 4

    # test that files were created as expected
    @test dx == true
    @test salts == true

    # test that when config has glob pattern in file name, ProjectConfig gets all files in directory.
    @test length(proj_config_batch.file_configs) == 9
    @test split(proj_config_batch.file_configs[3].filename, Sys.iswindows() ? "\\" : "/")[end] == "dx_2.csv"
end

@testset "primary identifiter" begin
    cfg = ProjectConfig("ehr_data_bad_pk.yml")

    @test_throws AssertionError deidentify(cfg)
end

@testset "hash seeding" begin
    cfg1 = ProjectConfig(test_file)
    cfg2 = ProjectConfig("ehr_data_alt_seed.yml")

    deid1 = deidentify(cfg1)
    deid1a = deidentify(cfg1)
    deid2 = deidentify(cfg2)

    @test deid1.salt == deid1a.salt
    @test deid1.salt != deid2.salt
end

@testset "create output directories" begin
    cfg = ProjectConfig("ehr_data_alt_paths.yml")
    deid = deidentify(cfg)

    @test isfile(joinpath(logpath_not_created,"ehr.log.0001"))
    @test isdir(joinpath(outputpath_not_created))
end

# TEAR DOWN
# this is necessary to ensure directories are deleted on Windows
GC.gc()

try
    rm(logpath, recursive=true, force=true)
catch
end

try
    rm(outputpath, recursive=true, force=true)
catch
end

try
    rm(joinpath(@__DIR__, "path"), recursive=true, force=true)
catch
end
# --------------------------


text:
# Wiktorsson method
# First published by Wiktorsson, 2001

struct Wiktorsson <: AbstractIteratedIntegralAlgorithm end

convorder(::Wiktorsson) = 1//1
errcoeff(m, h, ::Wiktorsson, ::MaxL2) = √(5*m)*h/(√12*π)
norv(m, n, ::Wiktorsson) = 2*m*n+(m^2-m)÷2

"""
    levyarea(W, n, alg::Wiktorsson)

Simulates an approximation of the iterated Itô-integrals ``\\int_0^1W_s\\otimes dW_s``
of the given ``m``-dimensional increment of a Wiener process with step size 1.
The parameter ``n`` specifies the number of terms in the approximation and thus determines the accuracy.
This is an efficient implementation of the algorithm proposed in [Wiktorsson, 2001](@ref wiktorsson2001).
It is based on the Fourier method from Milstein but incorporates an additional tail sum approximation.
The algorithm needs approximately ``2\\cdot m^2+2\\cdot m\\cdot n+m`` Float's.
The time complexity is ``\\mathcal{O}(m^2\\cdot n)``.
"""
function levyarea(W::AbstractVector{T}, n::Integer, alg::Wiktorsson) where {T}
    rng = default_rng()
    m = length(W)
    # Preallocate
    A = similar(W,m,m) # allocates m*m Floats
    G = similar(W,m,m) # allocates m*m Floats
    # 1. Simulate Xₖ and Yₖ and approximate stochastic area integral
    X = randn(rng, T, n, m) # allocates m*n Floats
    Y = randn(rng, T, m, n) # allocates m*n Floats
    Y .= (Y .- √(T(2)).*W) ./ (1:n)'
    mul!(A,Y,X)
    # 2.a Simulate Gₙ
    a = T(√(2*trigamma(n+1)))
    for j=1:m
        @inbounds G[j,j] = zero(T)
        for i=j+1:m
            g = a * randn(rng, T)
            @inbounds G[i,j] = g
            @inbounds G[j,i] = -g
            @inbounds A[i,j] += g
        end
    end
    # 2.b and add the tail-sum approximation
    A .+= inv(1+√(1+W'*W)) .* (G*W) .* W' # allocates m Floats
    # 3. Calculate the iterated integrals
    G .= inv(2*T(π)).*(A .- A') # reuse G to save allocations
    return G
end

max_stars_repo_name: mrbauff/MasterThesis

max_issues_repo_name: mrbauff/MasterThesis


text:
using Revise
using Soss
using Random

# Build the model

m = @model Prior,x,σ begin
    β ~ Prior
    yhat = β .* x
    y ~ For(eachindex(x)) do j
        Normal(yhat[j], σ)
    end
end

# Sample from generative model

Random.seed!(42);

x = randn(3);
Prior = Normal()
args = (Prior=Prior, x=x, σ=2.0);
truth = rand(m(args));

pairs(truth)

# Joint Distributions

pairs(args)

m(args)

# Evaluate log density

pairs(truth)

logpdf(m(args), truth)

# Sample from posterior

pairs(args)

post = dynamicHMC(m(args), (y=truth.y,));
particles(post)

# Determine predictive distribution

pred = predictive(m, :β)

# Sample from posterior predictive distribution

particles(post)

argspost = merge(args, particles(post));
pairs(argspost)

postpred = pred(argspost) |> rand;
pairs(postpred)

q = @model λ begin
β ~ Normal(λ,1)
end

# imp(m(args),q(λ=1.0),truth)
# importanceSample(m(args),q(λ=1.0),truth)

# particles(m(args))


text:
# https://fenicsproject.org/docs/dolfin/latest/python/demos/biharmonic/demo_biharmonic.py.html

module demo_biharmonic

using FenicsPy

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = true
parameters["form_compiler"]["optimize"] = true

# Make mesh ghosted for evaluation of DG terms
parameters["ghost_mode"] = "shared_facet"

# Create mesh and define function space
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "CG", 2)

# Define Dirichlet boundary
@pydef mutable struct DirichletBoundary <: dolfin.SubDomain
    function inside(self, x, on_boundary)
        return on_boundary
    end
end

py"""
from math import sin, pi
from dolfin import UserExpression
class Source(UserExpression):
    def __init__(self, **kwargs):
        super().__init__(**kwargs)
    def eval(self, values, x):
        values[0] = 4.0*pi**4*sin(pi*x[0])*sin(pi*x[1])
    def value_shape(self):
        return ()
"""

# Define boundary condition
u0 = Constant(0.0)
bc = DirichletBC(V, u0, SubDomain(DirichletBoundary()))

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)

# Define normal component, mesh size and right-hand side
h = CellDiameter(mesh)
h_avg = (h("+") + h("-"))/2.0
n = FacetNormal(mesh)
f = Expression(py"Source"(degree=2))

# Penalty parameter
alpha = Constant(8.0)

# Define bilinear form
a = inner(div(grad(u)), div(grad(v)))*dx -
   inner(avg(div(grad(u))), jump(grad(v), n))*dS -
   inner(jump(grad(u), n), avg(div(grad(v))))*dS +
   alpha/h_avg*inner(jump(grad(u),n), jump(grad(v),n))*dS

# Define linear form
L = f*v*dx

# Solve variational problem
u = FeFunction(V)
solve(a == L, u, bc)

# Save solution to file
file = File("biharmonic/biharmonic.pvd")
file << u

# Plot solution
#plot(u)

end # module demo_biharmonic


text:
using Gaugefields
function test()
    NX = 4
    NY = 4
    NZ = 4
    NT = 4
    L = [NX,NY,NZ,NT]

    nn = CovNeuralnet()
    ρ = [0.1]
    layername = ["plaquette"]
    st = STOUT_Layer(layername,ρ,L)
    push!(nn,st)

    ρ = [0.3]
    layername = ["plaquette"]
    st = STOUT_Layer(layername,ρ,L)
    push!(nn,st)
    show(nn)

    Dim = 4
    Nwing = 1
    NC = 3
    u1 = RandomGauges(NC,Nwing,L...)
    U = Array{typeof(u1),1}(undef,Dim)
    U[1] = u1
    for μ=2:Dim
        U[μ] = RandomGauges(NC,Nwing,L...)
    end

    @time Uout,Uout_multi,_ = calc_smearedU(U,nn)
end
test()


text:
module FEOperatorsTests

##
using Test
using Gridap

import Gridap: ∇
##

# Define manufactured functions
ufun(x) = x[1] + x[2]
ufun_grad(x) = VectorValue(1.0,1.0)
∇(::typeof(ufun)) = ufun_grad
bfun(x) = 0.0

# Construct the discrete model
model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0), partition=(4,4))

order = 1

include("FEOperatorsTestsMixin.jl")

include("FEOperatorsTestsRobin.jl")

order = 3

include("FEOperatorsTestsMixin.jl")

include("FEOperatorsTestsRobin.jl")

model = simplexify(model)

order = 1

include("FEOperatorsTestsMixin.jl")

include("FEOperatorsTestsRobin.jl")

order = 3

include("FEOperatorsTestsMixin.jl")

include("FEOperatorsTestsRobin.jl")

# Further tests
op = LinearFEOperator(V,U,t_Ω,t_ΓN,t_ΓR,t_ΓS)
op = LinearFEOperator(V,U,t_Ω)

end # module


text:
# Base Classes
@enum(Dynamics, dynamicFreeInitial, fixedInitial, steadyStateInitial, steadyState) # this works but it can't find it in the the "Model"

# Volumes
UnitVolume = Model(
    nParallel=1,
    V=1.0u"m^3",
    d=1000.0u"kg/m^3",
    cp=1.0u"J/(kg*K)",
    Q_gen=0.0u"W",
    T_reference  = 273.15u"K",
    U=Var(start=6000.0u"J"),
    port = HeatPort,
    equations = :[
        m = d*V
        U = m*cp*(T-T_reference)
        Qb_flow = port.Q_flow / nParallel + Q_gen
        der(U) = Qb_flow
        port.T = T
    ]
)
#    U=Var(start=6000.0u"J"),
#    T = Var(start=(20+273.15)u"K"),


text:
import SampledSignals
import SampledSignals: SampleBuf, domain

using DSP, DSP.Filters
using PaddedViews
using ProgressMeter

export signal, @rate, @samerateas, domain, analytic, isanalytic, samples
export padded, slide, toframe, rowview

SignalBase.nframes(x::SampleBuf) = SampledSignals.nframes(x)
SignalBase.framerate(x::SampleBuf) = SampledSignals.samplerate(x)
SignalBase.nchannels(x::SampleBuf) = SampledSignals.nchannels(x)
SignalBase.sampletype(x::SampleBuf) = eltype(x)

SignalBase.nframes(x::AbstractArray) = size(x,1)
SignalBase.framerate(x::AbstractArray) = 1.0
SignalBase.nchannels(x::AbstractArray) = size(x,2)
SignalBase.sampletype(x::AbstractArray) = eltype(x)

"""
$(SIGNATURES)
Creates a signal with frame rate `fs`.
"""
signal(x::AbstractArray, fs) = SampleBuf(x, inHz(fs))

"""
$(SIGNATURES)
Creates a signal with frame rate `fs`. If the original signal's frame rate is
the same as `fs`, this method simply returns the original signal. Otherwise, it
creates a new signal with the specified frame rate and data from the original
signal. Do note that this method does not resample the signal.
"""
signal(x::SampleBuf, fs) = fs == framerate(x) ? x : SampleBuf(x.data, inHz(fs))

"""
    @rate fs expr
Creates a signal from `expr` with frame rate `fs`. It provides syntactic sugar
on the [`signal`](@ref) method.

# Example:
```julia-repl
julia> x = @rate 44.1kHz randn(44100)
44100-frame, 1-channel SampleBuf{Float64, 1}
1.0s sampled at 44100.0Hz
```
"""
macro rate(fs, expr)
  :( signal($(esc(expr)), $(esc(fs))) )
end

"""
    @samerateas x expr
Creates a signal from `expr` with the same frame rate as signal `x`. This is
useful to preserve frame rate metadata across functions that do not return a
signal.

# Examples:
```julia-repl
julia> x = @rate 44.1kHz randn(44100)
44100-frame, 1-channel SampleBuf{Float64, 1}
1.0s sampled at 44100.0Hz

julia> using DSP
julia> y = filt([1.0,0.5], x)     # frame rate stripped by DSP.filt
44100-element Array{Float64,1}:
   ⋮

julia> y = @samerateas x filt([1.0,0.5], x)
44100-frame, 1-channel SampleBuf{Float64, 1}
1.0s sampled at 44100.0Hz
```
"""
macro samerateas(x, expr)
  :( signal($(esc(expr)), framerate($(esc(x)))))
end

"""
    @samerateas n x expr
Creates a signal from `expr` with the a frame rate `n` times that of signal `x`.

# Examples:
```julia-repl
julia> x = @rate 44.1kHz randn(44100)
44100-frame, 1-channel SampleBuf{Float64, 1}
1.0s sampled at 44100.0Hz

julia> y = @samerateas 1//2 x x[1:2:end]
22050-frame, 1-channel SampleBuf{Float64, 1}
1.0s sampled at 22050.0Hz

julia> using DSP
julia> y = @samerateas 2//3 x resample(x, 2//3)
29400-frame, 1-channel SampleBuf{Float64, 1}
1.0s sampled at 29400.0Hz
```
"""
macro samerateas(n, x, expr)
  :( signal($(esc(expr)), $(esc(n))*framerate($(esc(x)))))
end

"""
$(SIGNATURES)
Converts a signal to analytic representation.
"""
analytic(s::SampleBuf) = isanalytic(s) ? s : signal(hilbert(s)/√2.0, framerate(s))
analytic(s) = isanalytic(s) ? s : hilbert(s)/√2.0

"""
$(SIGNATURES)
Checks if a signal is analytic.
"""
isanalytic(s) = eltype(s) <: Complex

"""
$(SIGNATURES)
Gets the underlying samples in the signal.
"""
samples(s::SampleBuf) = s.data
samples(s) = s

"""
$(SIGNATURES)
Generates a padded view of a signal with optional delay/advance.
"""
function padded(s::AbstractVector{T}, padding; delay=0, fill=zero(T)) where {T, N}
  if length(padding) == 1
    left = padding
    right = padding
  else
    left = padding[1]
    right = padding[2]
  end
  PaddedView(fill, s, (1-left:length(s)+right,), (1+delay:delay+length(s),))
end

"""
    x[a:b,∘]
Generates a view of specified row range `a:b` of a vector or matrix `x`.
"""
Base.getindex(a::AbstractVector, i, ::typeof(∘)) = @view a[i]
Base.getindex(a::AbstractMatrix, i, ::typeof(∘)) = @view a[i,:]

"""
    slide(f::Function, s::AbstractVector, nframes, overlap=0, args...; showprogress=true)

Slides a window over a signal, processing each window. If the total number of frames
in the signal is not an integral multiple of `nframes`, the last incomplete
block of samples remains unprocessed.

The processing function receives a view on the original signal, and therefore
may modify the signal if desired.

# Examples:
```julia-repl
julia> x = signal(ones(1000), 8kHz);
julia> slide(x, 250) do x1, blknum, firstframe
         println(size(x1), ", ", blknum, ", ", firstframe)
       end
(250,), 1, 1
(250,), 2, 251
(250,), 3, 501
(250,), 4, 751

julia> slide(x, 250) do x1, blknum, firstframe
         x1 .= blknum
       end

julia> x[1], x[251], x[501], x[751]
(1.0, 2.0, 3.0, 4.0)
```
"""
function slide(f::Function, s::AbstractVecOrMat, nframes, overlap=0, args...; showprogress=true)
  @assert overlap < nframes "overlap must be less than nframes"
  n = size(s,1)
  m = nframes - overlap
  mmax = (n-nframes)÷m
  showprogress && (p = Progress(mmax+1, 1, "Processing: "))
  for j = 0:mmax
    s1 = s[j*m+1:j*m+nframes,∘]
    f(s1, j+1, j*m+1, args...)
    showprogress && next!(p)
  end
end

"""
    slide(f::Function, ::Type{T}, s::AbstractVector, nframes, overlap=0, args...; showprogress=true) where T

Slides a window over a signal, processing each window, and collecting the results of type `T`.
If the total number of frames in the signal is not an integral multiple of
`nframes`, the last incomplete block of samples remains unprocessed.

# Examples:
```julia-repl
julia> x = signal(ones(1000), 8kHz);
julia> slide(Float32, x, 250) do x1, blknum, firstframe
         sum(x1)*blknum
       end
4-element Array{Float32,1}:
  250.0
  500.0
  750.0
 1000.0

julia> slide(Tuple{Int,Float64}, x, 250) do x1, blknum, firstframe
          (blknum, sum(x1)*blknum)
        end
4-element Array{Tuple{Int64,Float64},1}:
  (1, 250.0)
  (2, 500.0)
  (3, 750.0)
  (4, 1000.0)
```
"""
function slide(f::Function, ::Type{T}, s::AbstractVecOrMat, nframes, overlap=0, args...; showprogress=true) where {T}
  @assert overlap < nframes "overlap must be less than nframes"
  n = size(s,1)
  m = nframes - overlap
  mmax = (n-nframes)÷m
  out = Array{T,1}(undef, 1+mmax)
  showprogress && (p = Progress(mmax+1, 1, "Processing: "))
  for j = 0:mmax
    s1 = s[j*m+1:j*m+nframes,∘]
    out[j+1] = f(s1, j+1, j*m+1, args...)
    showprogress && next!(p)
  end
  return out
end

"""
    slide(f::Function, ::Type{Array{T}}, noutput::Int, s::AbstractVector, nframes, overlap=0, args...; showprogress=true) where T

Slides a window over a signal, processing each window, and collecting the results
of type `Array{T}` of length `noutput`. If the total number of frames in the signal is
not an integral multiple of `nframes`, the last incomplete block of samples remains unprocessed.

# Examples:
```julia-repl
julia> x = signal(ones(1000), 8kHz);
julia> slide(Array{Float64}, 2, x, 250) do x1, blknum, firstframe
          [sum(x1), prod(x1)]
       end
4×2 Array{Float64,2}:
 250.0  1.0
 250.0  1.0
 250.0  1.0
 250.0  1.0
```
"""
function slide(f::Function, ::Type{Array{T}}, noutput::Int, s::AbstractVecOrMat, nframes, overlap=0, args...; showprogress=true) where {T}
  @assert overlap < nframes "overlap must be less than nframes"
  n = size(s,1)
  m = nframes - overlap
  mmax = (n-nframes)÷m
  out = Array{T,2}(undef, (1+mmax, noutput))
  showprogress && (p = Progress(mmax+1, 1, "Processing: "))
  for j = 0:mmax
    s1 = s[j*m+1:j*m+nframes,∘]
    out[j+1,:] = f(s1, j+1, j*m+1, args...)
    showprogress && next!(p)
  end
  return out
end

"""
$(SIGNATURES)
Converts time to signal frame number.

# Examples:
```julia-repl
julia> x = signal(randn(2000), 8kHz);
julia> toframe(0.2s, x)
1601

julia> toframe([0.2s, 201ms], x)
2-element Array{Int64,1}:
 1601
 1609

julia> toframe(0.2:0.01:0.3, x)
 11-element Array{Int64,1}:
  1601
  1681
  1761
   ⋮
```
"""
toframe(t, s::SampleBuf) = 1 .+ round.(Int, inseconds.(t)*framerate(s))


text:
### Not altering the tree length ###

"""
    slide!(root::T) where T<:AbstractNode

This function performs a slide move on an intermediate node. The node is moved
upwards or downwards on the path specified by its mother and one of its
daughters.

* `root` : root Node of tree.
"""
function slide!(root::T) where {T<:AbstractNode}

    available = [node.num for node in post_order(root)]
    n = rand(available)
    target::T = find_num(root, n)

    while target.nchild == 0 || target.root
        n = rand(available)
        target = find_num(root, n)
    end

    # proportion of slide move is randomly selected
    proportion::Float64 = rand()
    # pick a random child
    child::T = rand(target.children)

    # calculate and set new values
    move!(target, child, proportion)
    nothing

end # function slide!

"""
    slide(root::T)::T where T<:AbstractNode

This function performs a slide move on an intermediate node. The node is moved
upwards or downwards on the path specified by its mother and one of its
daughters.

Returns root Node of new tree.

* `root` : root Node of tree.
"""
function slide(root::T)::T where {T<:AbstractNode}
    new_root = deepcopy(root)
    slide!(new_root)
    return new_root
end

"""
    swing!(root::T) where T<:AbstractNode

This function performs a swing node. A random non-leave node is selected and
moved along the path specified by its two children.

* `root` : root Node of tree.
"""
function swing!(root::T) where {T<:AbstractNode}

    available = [node.num for node in post_order(root)]
    n = rand(available)
    target::T = find_num(root, n)

    while target.nchild < 2
        n = rand(available)
        target = find_num(root, n)
    end
    proportion::Float64 = rand()

    child1 = target.children[1]
    child2 = target.children[2]

    # calculate and set new values
    move!(child1, child2, proportion)
    nothing
end # function swing!

"""
    swing(root::T)::T where T<:AbstractNode

This function performs a swing node. A random non-leave node is selected and
moved along the path specified by its two children.

Returns root Node of new tree.

* `root` : root Node of tree.
"""
function swing(root::T)::T where {T<:AbstractNode}
    new_root = deepcopy(root)
    swing!(new_root)
    return new_root
end

"""
    move!(node1::T, node2::T, proportion::Float64) where T <:AbstractNode

Change the incoming length of node1 and node2 while keeping their combined length
constant.

* `node1` : Node whose inc\\_length will be modified; this node's inc\\_length will be the total inc\\_length of both nodes, times proportion.

* `node2` : Node whose inc\\_length will be modified; this node's inc\\_length will be the remainder of total - the new inc\\_length value of node1.

* `proportion` : Float64, determines proportion of the inc\\_length of both nodes assigned to node1.
"""
function move!(node1::T, node2::T, proportion::Float64) where {T<:AbstractNode}
    total::Float64 = node1.inc_length + node2.inc_length
    fp::Float64 = total * proportion
    sp::Float64 = total - fp
    node1.inc_length = fp
    node2.inc_length = sp
end # function move!

### Altering the length of the tree ###

"""
    change_edge_length!(root::T) where T <:AbstractNode

Pick a random node and increase or decrease its length randomly.

* `root` : root node of tree.
"""
function change_edge_length!(root::T) where {T<:AbstractNode}
    available = [node.num for node in post_order(root)]
    n = rand(available)
    target::T = find_num(root, n)
    while target.root
        n = rand(available)
        target = find_num(root, n)
    end
    factor = abs(randn())
    target.inc_length *= factor
end


text:
using Random
using NeuralVerification;

"""
    make_random_network(layer_sizes::Vector{Int}, [min_weight = -1.0], [max_weight = 1.0], [min_bias = -1.0], [max_bias = 1.0], [rng = 1.0])
    read_layer(output_dim::Int, f::IOStream, [act = ReLU()])

Generate a network with random weights and bias. The first layer is treated as the input.
The values for the weights and bias will be uniformly drawn from the range between min_weight
and max_weight and min_bias and max_bias respectively. The last layer will have an ID()
activation function and the rest will have ReLU() activation functions. Allow a random number
generator(rng) to be passed in. This allows for seeded random network generation.
"""
function make_random_network(layer_sizes::Vector{Int}, min_weight = -1.0, max_weight = 1.0, min_bias = -1.0, max_bias = 1.0, rng=MersenneTwister(0))
    # Create each layer based on the layer_size
    layers = []
    for index in 1:(length(layer_sizes)-1)
        cur_size = layer_sizes[index]
        next_size = layer_sizes[index+1]
        # Use Id activation for the last layer - otherwise use ReLU activation
        if index == (length(layer_sizes)-1)
            cur_activation = NeuralVerification.Id()
        else
            cur_activation = NeuralVerification.ReLU()
        end

        # Dimension: num_out x num_in
        cur_weights = min_weight .+ (max_weight - min_weight) * rand(rng, Float64, (next_size, cur_size))
        cur_weights = reshape(cur_weights, (next_size, cur_size)) # for edge case where 1 dimension is equal to 1 this keeps it from being a 1-d vector

        # Dimension: num_out x 1
        cur_bias = min_bias .+ (max_bias - min_bias) * rand(rng, Float64, (next_size))
        push!(layers, NeuralVerification.Layer(cur_weights, cur_bias, cur_activation))
    end

    return Network(layers)
end


text:
using TypeClasses

# MonoidAlternative
# =================

TypeClasses.neutral(::Type{Writer{Accumulator, Value}}) where {Accumulator, Value} = neutral(Accumulator) => neutral(Value)
TypeClasses.combine(w1::Writer, w2::Writer) = Writer(combine(w1.acc, w2.acc), combine(w1.value, w2.value))

# we don't implement orelse, as it is commonly meant on container level, but there is no obvious failure semantics here

# FunctorApplicativeMonad
# =======================

TypeClasses.foreach(f, p::Writer) = f(p.value); nothing
TypeClasses.map(f, p::Writer) = Writer(p.acc, f(p.value))

# pure needs Neutral on First
function TypeClasses.pure(::Type{<:Writer{Acc}}, a) where {Acc}
  Writer(neutral(Acc), a)
end
# We fall back to assume that the general writer uses Option values in order to have a neutral value
# this should be okay, as it is canonical extension for any Semigroup to an Monoid
function TypeClasses.pure(::Type{<:Writer}, a)
  Writer(neutral, a)
end

# Writer always defines `combine` on `acc`
function TypeClasses.ap(f::Writer, a::Writer)
  Writer(combine(f.acc, a.acc), f.value(a.value))
end

# we cannot overload this generically, because `Base.map(f, ::Vector...)` would get overwritten as well (even without warning surprisingly)
TypeClasses.map(f, a::Writer, b::Writer, more::Writer...) = mapn(f, a, b, more...)

function TypeClasses.flatmap(f, a::Writer)
  # we intentionally only convert the container to Writer, and leave the responsibility for the accumulator to the `combine` function
  # (one reason is that all monads are only converted on the container side, another is that it seems quite difficult to convert MyType{Option} correctly, because it is a Union type.)
  nested_writer = convert(Writer, f(a.value))
  Writer(combine(a.acc, nested_writer.acc), nested_writer.value)
end

# flip_types
# ==========

function TypeClasses.flip_types(a::Writer)
  TypeClasses.map(x -> Writer(a.acc, x), a.value)
end


text:
@testset "MaskedSetMosaic" begin
    using OptEnrichedSetCover: MaskOverlap, MaskedSetMosaic

    @testset "empty" begin
        sm = SetMosaic(Set{Symbol}[])
        msm = mask(sm, [Set{Symbol}()])
        @test msm isa MaskedSetMosaic
        @test originalmosaic(msm) == sm
        @test nelements(msm) == 0
        @test nsets(msm) == 0
        @test nmasks(msm) == 1
        @test msm.ix2experiment == [1]
    end

    @testset "empty but with elements" begin
        sm = SetMosaic(Set{Symbol}[], Set([:a, :b]))
        msm = mask(sm, [Set([:a])])
        @test originalmosaic(msm) == sm
        @test nelements(msm) == 2
        @test nsets(msm) == 0
        @test nmasks(msm) == 1
        @test nmasked(msm, 1) == 1
        @test nunmasked(msm, 1) == 1
        @test msm.ix2experiment == [1]

        # named
        nmsm = mask(sm, Dict(:X => Set([:a])))
        @test originalmosaic(nmsm) == sm
        @test nelements(nmsm) == 2
        @test nsets(nmsm) == 0
        @test nmasks(nmsm) == 1
        @test nmasked(nmsm, :X) == 1
        @test nunmasked(nmsm, :X) == 1
        @test nmsm.ix2experiment == [:X]
    end

    @testset "[:a :b] [:c :d] [:a :b :c :d], mask=[:a :b]" begin
        sm = SetMosaic([Set([:a, :b]), Set([:c, :d]), Set([:a, :b, :c, :d])])

        msm = mask(sm, [Set([:a, :b])])

        @test nelements(msm) == 4
        @test nmasked(msm, 1) == 2
        @test nunmasked(msm, 1) == 2
        @test nsets(msm) == 2
        @test nmasks(msm) == 1
        @test overlap(msm, 1, 1) == MaskOverlap(2, 0, true)
        @test overlap(msm, 1, 3) == MaskOverlap(2, 2, true)
        @test overlap(msm, 1, 2) === missing

        msm_copy = copy(msm)
        @test nelements(msm_copy) == nelements(msm)
        @test nmasks(msm_copy) == nmasks(msm)
        @test msm_copy.original === msm.original
        @test msm_copy.elmasks !== msm.elmasks
        @test msm_copy.elmasks == msm.elmasks
        @test msm_copy.total_masked !== msm.total_masked
        @test msm_copy.total_masked == msm.total_masked
        @test msm_copy.setXexp_olaps !== msm.setXexp_olaps
        @test msm_copy.setXexp_olaps == msm.setXexp_olaps

        # mask with nonexisting element
        msm2 = mask(sm, [Set([:a, :b, :g])])
        @test nelements(msm2) == 4
        @test nmasked(msm2, 1) == 2
        @test nunmasked(msm2, 1) == 2
        @test nsets(msm2) == 2
        @test overlap(msm2, 1, 1) == MaskOverlap(2, 0, true)
        @test overlap(msm2, 1, 3) == MaskOverlap(2, 2, true)

        # mask with max_overlap_logpvalue, [:a :b :c :d] is excluded
        msm3 = mask(sm, [Set([:a, :b])], max_overlap_logpvalue=-0.1)
        @test nelements(msm3) == 4
        @test nmasked(msm3, 1) == 2
        @test nunmasked(msm3, 1) == 2
        @test nsets(msm3) == 1
        @test overlap(msm3, 1, 1) == MaskOverlap(2, 0, true)
        @test overlap(msm3, 1, 3) === missing
    end

    @testset "A=[:a :b] B=[:c :d] C=[:a :b :c :d], mask=[:a :b]" begin
        sm = SetMosaic(Dict(:A=>Set([:a, :b]), :B=>Set([:c, :d]), :C=>Set([:a, :b, :c, :d])))
        msm = mask(sm, [Set([:a, :b])])

        @test nelements(msm) == 4
        @test nsets(msm) == 2
        @test nmasks(msm) == 1
        @test overlap(msm, 1, :A) == MaskOverlap(2, 0, true)
        @test overlap(msm, 1, :B) === missing
        @test overlap(msm, 1, :C) == MaskOverlap(2, 2, true)
    end

    @testset "multimask: A=[:a :b] B=[:c :d] C=[:a :b :c :d], mask=[[:a :b] [:a :b :c] [:a :b :d]]" begin
        sm = SetMosaic(Dict(:A=>Set([:a, :b]), :B=>Set([:c, :d]), :C=>Set([:a, :b, :c, :d])))
        msm = mask(sm, [Set([:a, :b]), Set([:a, :b, :c]), Set([:a, :b, :d])], min_nmasked=1)

        @test nelements(msm) == 4
        @test nsets(msm) == 3
        @test nmasks(msm) == 3
        @test_throws KeyError overlap(msm, 4, :A)
        @test_throws KeyError overlap(msm, 1, :D)
        @test overlap(msm, 1, :A) == MaskOverlap(2, 0, true)
        @test overlap(msm, 2, :A) == MaskOverlap(2, 0, true)
        @test overlap(msm, 3, :A) == MaskOverlap(2, 0, true)
        @test overlap(msm, 1, :B) === MaskOverlap(0, 2, true)
        @test overlap(msm, 2, :B) == MaskOverlap(1, 1, true)
        @test overlap(msm, 3, :B) == MaskOverlap(1, 1, true)
        @test overlap(msm, 1, :C) == MaskOverlap(2, 2, true)
        @test overlap(msm, 2, :C) == MaskOverlap(3, 1, true)
        @test overlap(msm, 3, :C) == MaskOverlap(3, 1, true)

        msm2 = mask(sm, [Set([:a, :b]), Set([:a, :b, :c]), Set([:a, :b, :d])], min_nmasked=2)
        @test nsets(msm2) == 2
        @test overlap(msm2, 1, :A) == MaskOverlap(2, 0, true)
        @test overlap(msm2, 2, :A) == MaskOverlap(2, 0, true)
        @test overlap(msm2, 3, :A) == MaskOverlap(2, 0, true)
        @test overlap(msm2, 1, :B) === missing
        @test overlap(msm2, 2, :B) === missing
        @test overlap(msm2, 3, :B) === missing
        @test overlap(msm2, 1, :C) == MaskOverlap(2, 2, true)
        @test overlap(msm2, 2, :C) == MaskOverlap(3, 1, true)
        @test overlap(msm2, 3, :C) == MaskOverlap(3, 1, true)

        nmsm2 = mask(sm, Dict(:X => Set([:a, :b]), :Y => Set([:a, :b, :c]), :Z => Set([:a, :b, :d])), min_nmasked=2)
        @test nsets(nmsm2) == 2
        @test_throws KeyError overlap(nmsm2, :U, :A)
        @test_throws KeyError overlap(nmsm2, :X, :D)
        @test overlap(nmsm2, :X, :A) == MaskOverlap(2, 0, true)
        @test overlap(nmsm2, :Y, :A) == MaskOverlap(2, 0, true)
        @test overlap(nmsm2, :Z, :A) == MaskOverlap(2, 0, true)
        @test overlap(nmsm2, :X, :B) === missing
        @test overlap(nmsm2, :Y, :B) === missing
        @test overlap(nmsm2, :Z, :B) === missing
        @test overlap(nmsm2, :X, :C) == MaskOverlap(2, 2, true)
        @test overlap(nmsm2, :Y, :C) == MaskOverlap(3, 1, true)
        @test overlap(nmsm2, :Z, :C) == MaskOverlap(3, 1, true)
    end
end


text:
using naumer_Dimensionality_2022
using SafeTestsets

# Break up the testing into a per-file setup

@safetestset "core.jl" begin
    include("core_test.jl")
end

@safetestset "systems.jl" begin
    include("systems_test.jl")
end


text:
using Unitful
using DFTK

lattice = load_lattice(joinpath(@__DIR__, "Cr19.json"))
atoms   = load_atoms(  joinpath(@__DIR__, "Cr19.json"))
atoms   = attach_pseudos(atoms, Cr="hgh/pbe/cr-q14.hgh")
magnetic_moments = build_magnetic_moments(atoms; Cr=5.0)

cases = standard_cases()
cases["standard_02"] = (damping=DFTK.FixedDamping(0.2), tol=1e-12, )
cases["standard_03"] = (damping=DFTK.FixedDamping(0.3), tol=1e-12, )
cases["standard_04"] = (damping=DFTK.FixedDamping(0.4), tol=1e-12, )

config[@__DIR__] = (
    directory        =  @__DIR__,
    lattice          = lattice,
    atoms            = atoms,
    magnetic_moments = magnetic_moments,
    #
    kwargs_model     = (smearing=Smearing.Gaussian(), temperature=0.01),
    Ecut             = 30,
    kgrid            = kgrid_size_from_minimal_spacing(lattice, 2π * 0.04 / u"Å"),
    supersampling    = 1.5,
    potential        = nothing,
    kwargs_scf       = (; ),
    mixing           = KerkerMixing(),
    cases            = cases,
)

max_stars_repo_name: mfherbst/supporting-adaptive-damping

max_issues_repo_name: mfherbst/supporting-adaptive-damping

max_forks_repo_name: mfherbst/supporting-adaptive-damping


text:
mutable struct Callbacks
    cbs::Vector
    runner::CallbackRunner
    graph::SimpleDiGraph
    initialized::Bool
end
Callbacks(cbs, runner=LinearRunner()) = Callbacks(cbs, runner, callbackgraph(cbs), false)

init!(cbs::Callbacks, learner) = foreach(cb -> init!(cb, learner), cbs.cbs)

mutable struct Learner
    model
    data::PropDict
    optimizer
    lossfn
    params
    step::PropDict
    callbacks::Callbacks
    cbstate::PropDict
end

"""
    Learner(model, data, optimizer, lossfn, [callbacks...; kwargs...])

Holds and coordinates all state of the training. `model` is trained by
optimizing `lossfn` with `optimizer` on `data`.

## Arguments

- `model`: A Flux.jl model or a `NamedTuple` of models.
- `data`: Data iterators. A 2-tuple will be treated as `(trainingdataiter, validdataiter)`.
    You can also pass in an empty tuple `()` and use the [`epoch!`](#) method with a
    `dataiter` as third argument.

    A data iterator is an iterable over batches. For regular supervised training,
    each batch should be a tuple `(xs, ys)`.
- `lossfn`: Function with signature `lossfn(model(x), y) -> Number`
- `optimizer`
- `callbacks...`: Any other unnamed arguments are callbacks

## Keyword arguments

- `usedefaultcallbacks = true`: Whether to add some basic callbacks. Included
    are [`Metrics`](#), [`Recorder`](#), [`ProgressPrinter`](#),
    [`StopOnNaNLoss`](#), and [`MetricsPrinter`](#).
- `cbrunner = LinearRunner()`: Callback runner to use.

## Fields

*(Use this as a reference when implementing callbacks)*

- `model`, `optimizer`, and `lossfn` are stored as passed in
- `data` is a `PropDict` of data iterators, usually `:training` and `:validation`.
- `params`: An instance of `model`'s parameters of type `Flux.Params`. If `model` is
    a `NamedTuple`, then `params` is a `NamedTuple` as well.
- `step::`[`PropDict`](#): State of the last step. Contents depend on the last run
    [`Phase`](#).
- `cbstate::`[`PropDict`](#): Special state container that callbacks can
    save state to for other callbacks. Its keys depend on what callbacks
    are being used. See the [custom callbacks guide](/documents/docs/callbacks/custom.md)
    for more info.
"""
function Learner(
        model, data, optimizer, lossfn, callbacks::Vararg{<:Callback};
        usedefaultcallbacks=true, cbrunner=LinearRunner()
    )
    callbacks = collect(Callback, callbacks)

    if usedefaultcallbacks
        for cb in defaultcallbacks()
            if !any(typeof(cb) .== typeof.(callbacks))
                push!(callbacks, cb)
            end
        end
    end

    cbs = Callbacks(callbacks, cbrunner)

    learner = Learner(
        model,
        _dataiters(data),
        optimizer,
        lossfn,
        paramsrec(model),
        PropDict(),
        cbs,
        PropDict())
    init!(cbs, learner)
    return learner
end

Base.show(io::IO, learner::Learner) = print(io, "Learner()")

defaultcallbacks()::Vector{AbstractCallback} = [
    ProgressPrinter(),
    MetricsPrinter(),
    StopOnNaNLoss(),
    Recorder(),
    Metrics(),
]

#  Callback handling
handle(event, learner, phase) = handle(learner.callbacks.runner, event, phase, learner)

# Other

phasedataiter(::AbstractTrainingPhase) = :training
phasedataiter(::AbstractValidationPhase) = :validation

function model!(learner, model)
    learner.model = model
    learner.params = paramsrec(model)
end

_dataiters(d::PropDict) = d
_dataiters(t::NamedTuple) = PropDict(pairs(t))
function _dataiters(t::Tuple)
    if length(t) == 0
        return PropDict(Dict{Symbol,Any}())
    elseif length(t) == 1
        return _dataiters((training = t[1]))
    elseif length(t) == 2
        return _dataiters((training = t[1], validation = t[2]))
    else
        error("Please pass a `NamedTuple` or `PropDict` as `data`.")
    end
end

paramsrec(m) = Flux.params(m)
paramsrec(t::Union{Tuple,NamedTuple}) = map(paramsrec, t)

# Callback utilities


text:
"seek h by Laplacian approximation to L/q, L: log of true_density, q: current_approx.
Laplacian: μ = argmax (log L - Log q), Σ = (- Hessian of log L at μ)^(-1)
The new compoent is N(μ, Σ/2).

Univariate version."
function get_new_component_laplacian(ll::Function, current_approx::UnivariateDistribution;
                                    init_approx = false, lower = -100., upper = 100.,
                                    show_trace=false, verbose=false)
    # get init point
    opt_init = rand(current_approx)
    if init_approx
        # get the 1st component: scaled Laplacian approximation
        f(x) = -ll(x)
    else
        # maximize log L - log q
        # Julia optimization routine uses "Brent's" for univariate optimization
        # lower and upper bound must be provided
        f = stable_negative_log_residual(ll, current_approx)  # obj to minimize
    end

    opt_result = optimize(f, lower, upper, method = GoldenSection())
    # check convergence
    if verbose
        println(opt_result)
    end
    converged(opt_result) || error("Optimization not converged", opt_result)
    # get the Hessian
    η = opt_result.minimizer
    # Finite Differencing
    H = second_derivative(f, η)
    println(H)
    new_component = laplacian_fit(get_component_distribution(current_approx), η, H)
    if init_approx
        println(@sprintf("Initial approximation: N(μ=%f, σ=%f)", new_component.μ, new_component.σ))
    end
    return new_component
end

"seek h by Laplacian approximation to L/q, L: log of true_density, q: current_approx.
Laplacian: μ = argmax (log L - Log q), Σ = (- Hessian of log L at μ)^(-1)
The new compoent is N(μ, Σ/2).

Multivariate version, random initialization."
function get_new_component_laplacian(ll::Function, current_approx::MultivariateDistribution;
                                    init_approx = false, diagonal = false,
                                    max_iter=1000, show_trace=false, verbose=false)

    opt_init = rand(current_approx)
    return get_new_component_laplacian(ll, current_approx, opt_init;
                                       init_approx = init_approx,
                                       diagonal = diagonal,
                                       max_iter=max_iter,
                                       show_trace=show_trace, verbose=verbose)
end

"with supplied initialization."
function get_new_component_laplacian{T <: Real}(ll::Function,
                                                current_approx::MultivariateDistribution,
                                                x_init::AbstractVector{T};
                                                init_approx = false,
                                                diagonal = false,
                                                max_iter=1000, show_trace=false, verbose=false)
    d = length(current_approx)
    if init_approx
        # get the 1st component: scaled Laplacian approximation
        f(x) = -ll(x)
    else
        # maximize log L - log q
        f = stable_negative_log_residual(ll, current_approx)  # obj to minimize
    end
    opt_result = optimize(f, x_init, LBFGS(),
                Options(iterations = max_iter,
                                    show_trace = show_trace))
    # check convergence
    converged(opt_result) || error("Optimization not converged", opt_result)
    if verbose
        println(opt_result)
    end
    # get the Hessian
    η = opt_result.minimizer
    # Finite Differencing, H should be positive definite
    if d > 20
        println("(.. estimating hessian ..)")
    end
    H = hessian(f, η)
    H = (H + H') / 2.
    if diagonal
        H = diagm(diag(H))
    end
    if d > 20
        if diagonal
            println("(.. diagonal hessian approximated ..)")
        else
            println("(.. cholesky on hessian ..)")
        end
    end
    H = PDMat(H)
    new_component = laplacian_fit(get_component_distribution(current_approx), η, H)
    if init_approx
        println("Initial approximation:N(μ = ", new_component.μ, ", Σ = ", new_component.Σ)
    end
    return new_component
end

"""
Laplacian fit for multivariate Gaussian

* η: argmax log L - log q
* H: Hessian of (-log L + log q) at η

return h = MvNormal(η, inv(H)/2)
"""
function laplacian_fit{D <: MvNormal, T <: Real}(::Type{D}, η::Vector{T}, H::PDMat)
    @assert length(η) == size(H, 1)
    Σ = inv(H) / 2.0
    # Σ = inv(H)
    all(isfinite.(η)) || error("Invalid number encountered in Laplacian fit to Gaussian")
    if length(η) > 20
        println("(.. new Gaussian component determined ..)")
    end
    return MvNormal(η, Σ) # optimal scaling
end

"Laplacian fit for univariate Gaussian"
function laplacian_fit{D <: Normal, T <: Real}(::Type{D}, η::T, H::T)
    # σ2 = 1.0 / H / 2.0
    σ2 = 1.0 / H / 2.0
    if !(isfinite.(σ2) && σ2 > 0)
        warn("ill H = ", H)
        return Normal(η, 1e-4)
    end
    return Normal(η, sqrt(σ2)) # optimal scaling, 2nd argument is sd
end

max_stars_repo_name: richardkwo/boosting-VB

max_issues_repo_name: richardkwo/boosting-VB

max_forks_repo_name: richardkwo/boosting-VB


text:
#####
##### Gradient AD implementation using ReverseDiff
#####

import .ReverseDiff

struct ReverseDiffLogDensity{L,C} <: ADGradientWrapper
    ℓ::L
    compiledtape::C
end

"""
    ADgradient(:ReverseDiff, ℓ)
    ADgradient(Val(:ReverseDiff), ℓ)

Gradient using algorithmic/automatic differentiation via ReverseDiff.
"""
ADgradient(::Val{:ReverseDiff}, ℓ) = begin
    f = _logdensity_closure(ℓ)
    x = rand(dimension(ℓ)) #init random parameters
    tape = ReverseDiff.GradientTape(f, x)
    compiledtape = ReverseDiff.compile(tape)
    ReverseDiffLogDensity(ℓ, compiledtape)
end

Base.show(io::IO, ∇ℓ::ReverseDiffLogDensity) = print(io, "ReverseDiff AD wrapper for ", ∇ℓ.ℓ)

function logdensity_and_gradient(∇ℓ::ReverseDiffLogDensity, x::AbstractVector{T}) where {T}
    @unpack ℓ, compiledtape = ∇ℓ
    buffer = _diffresults_buffer(ℓ, x)
    result = ReverseDiff.gradient!(buffer, compiledtape, x)
    _diffresults_extract(result)
end


text:
"""
    get_quantile_scores(scores, quantiles = 0.0:0.01:1.0)

return the quantiles of the given N dimensional anomaly `scores` cube. `quantiles` (default: `quantiles = 0.0:0.01:1.0`) is a Float range of quantiles.
Any score being greater or equal `quantiles[i]` and beeing smaller than `quantiles[i+1]` is assigned to the respective quantile `quantiles[i]`.

# Examples
```
julia> scores1 = rand(10, 2)
julia> quantile_scores1 = get_quantile_scores(scores1)
```
"""
function get_quantile_scores(scores::AbstractArray{tp,N}, quantiles::StepRangeLen{Float64} = 0.0:0.01:1.0) where {tp,N}
  quantile_scores = zeros(Float64, size(scores))
  get_quantile_scores!(quantile_scores, scores, quantiles)
  return quantile_scores
end

"""
    get_quantile_scores!{tp,N}(quantile_scores::AbstractArray{Float64, N}, scores::AbstractArray{tp,N}, quantiles::StepRangeLen{Float64} = 0.0:0.01:1.0)

return the quantiles of the given N dimensional `scores` array into a preallocated `quantile_scores` array, see `get_quantile_scores()`.
"""
function get_quantile_scores!(quantile_scores::AbstractArray{Float64, N}, scores::AbstractArray{tp,N}, quantiles::StepRangeLen{Float64} = 0.0:0.01:1.0; return_thresholds::Bool = false) where {tp,N}
  LENGTH = length(scores)
  thresholds = quantile(reshape(scores, prod(size(scores))), collect(quantiles))
  n_quants = size(thresholds, 1)
  for j = 1:LENGTH
      if(scores[j] <= thresholds[1])
        quantile_scores[j] = quantiles[1]
      end
      for i = 2:n_quants
        if(scores[j] > thresholds[i-1] && scores[j] <= thresholds[i])
          quantile_scores[j] = quantiles[i]
        end
      end
  end
  if return_thresholds == false
    return quantile_scores
  else
    return thresholds
  end
end

"""
    compute_ensemble(m1_scores, m2_scores[, m3_scores, m4_scores]; ensemble = "mean")

compute the mean (`ensemble = "mean"`), minimum (`ensemble = "min"`), maximum (`ensemble = "max"`) or median (`ensemble = "median"`) of the given anomaly scores.
Supports between 2 and 4 scores input arrays (`m1_scores, ..., m4_scores`). The scores of the different anomaly detection algorithms should be somehow comparable,
e.g., by using `get_quantile_scores()` before.

# Examples
```
julia> using MultivariateAnomalies
julia> scores1 = rand(10, 2)
julia> scores2 = rand(10, 2)
julia> quantile_scores1 = get_quantile_scores(scores1)
julia> quantile_scores2 = get_quantile_scores(scores2)
julia> compute_ensemble(quantile_scores1, quantile_scores2, ensemble = "max")
```
"""
function compute_ensemble(m1_scores::Array{Tp, N}, m2_scores::Array{Tp, N}; ensemble = "mean") where {Tp, N}
    @assert any(ensemble .== ["mean","min","max","median"])

    scores = cat(m1_scores, m2_scores, dims = N+1);

    if ensemble == "mean"
      ensemble_scores = dropdims(mean(scores, dims = N+1), dims = N+1)
    end

    if ensemble == "median"
      ensemble_scores = dropdims(median(scores, dims = N+1), dims = N+1)
    end

    if ensemble == "max"
      ensemble_scores = dropdims(maximum(scores, dims = N+1), dims = N+1)
    end

    if ensemble == "min"
      ensemble_scores = dropdims(minimum(scores, dims = N+1), dims = N+1)
    end
  return ensemble_scores
end

function compute_ensemble(m1_scores::Array{Tp, N}, m2_scores::Array{Tp, N}, m3_scores::Array{Tp, N}; ensemble = "mean") where {Tp, N}
    @assert any(ensemble .== ["mean","min","max","median"])

    scores = cat(m1_scores, m2_scores, m3_scores, dims = N+1);

    if ensemble == "mean"
      ensemble_scores = dropdims(mean(scores, dims = N+1), dims = N+1)
    end

    if ensemble == "median"
      ensemble_scores = dropdims(median(scores, dims = N+1), dims = N+1)
    end

    if ensemble == "max"
      ensemble_scores = dropdims(maximum(scores, dims = N+1), dims = N+1)
    end

    if ensemble == "min"
      ensemble_scores = dropdims(minimum(scores, dims = N+1), dims = N+1)
    end
  return ensemble_scores
end

function compute_ensemble(m1_scores::Array{Tp, N}, m2_scores::Array{Tp, N}, m3_scores::Array{Tp, N}, m4_scores::Array{Tp, N}; ensemble = "mean") where {Tp, N}
    @assert any(ensemble .== ["mean","min","max","median"])

    scores = cat(m1_scores, m2_scores, m3_scores, m4_scores, dims = N+1);

    if ensemble == "mean"
      ensemble_scores = dropdims(mean(scores, dims = N+1), dims = N+1)
    end

    if ensemble == "median"
      ensemble_scores = dropdims(median(scores, dims = N+1), dims = N+1)
    end

    if ensemble == "max"
      ensemble_scores = dropdims(maximum(scores, dims = N+1), dims = N+1)
    end

    if ensemble == "min"
      ensemble_scores = dropdims(minimum(scores, dims = N+1), dims = N+1)
    end
  return ensemble_scores
end

###################################
# end


text:
# External packages
using Revise, DelimitedFiles, SymPy

# Add Exasim to Julia search path
cdir = pwd(); ii = findlast("Exasim", cdir);
include(cdir[1:ii[end]] * "/Installation/setpath.jl");

# Exasim packages
using Preprocessing, Mesh, Gencode, Postprocessing

# create pde structure and mesh structure
pde, mesh = Preprocessing.initializeexasim();

# Define PDE model: governing equations, initial solutions, and boundary conditions
pde.model = "ModelC";            # ModelC, ModelD, ModelW
include("pdemodel.jl");          # include the PDE model file

# Set discretization parameters, physical parameters, and solver parameters
pde.porder = 3;          # polynomial degree
pde.torder = 2;          # time-stepping order of accuracy
pde.nstage = 2;          # time-stepping number of stages
pde.dt = 0.02*ones(200);   # time step sizes
pde.soltime = collect(10:10:length(pde.dt)); # steps at which solution are collected
pde.visdt = 1; # visualization timestep size

gam = 1e4;                      # gravity
pde.physicsparam = [gam 0.0];
pde.tau = [1.0];          # DG stabilization parameter
pde.GMRESrestart=15;            # number of GMRES restarts
pde.linearsolvertol=1e-12;      # GMRES tolerance
pde.linearsolveriter=16;        # number of GMRES iterations
pde.precMatrixType=2;           # preconditioning type
pde.NLtol = 1e-12;              # Newton tolerance
pde.NLiter=2;                   # Newton iterations

# Choose computing platform and set number of processors
#pde.platform = "gpu";           # choose this option if NVIDIA GPUs are available
pde.mpiprocs = 1;                # number of MPI processors

# create a linear mesh for a square domain
mesh.p, mesh.t = Mesh.SquareMesh(64,64,1); # a mesh of 8 by 8 quadrilaterals
mesh.p = (4*pi)*mesh.p .- 2*pi;
# expressions for disjoint boundaries
mesh.boundaryexpr = [p -> (p[2,:] .< -2*pi+1e-3), p -> (p[1,:] .> 2*pi-1e-3), p -> (p[2,:] .> 2*pi-1e-3), p -> (p[1,:] .< -2*pi+1e-3)];
mesh.boundarycondition = [1 1 1 1]; # Set boundary condition for each disjoint boundary
mesh.periodicexpr = [2 p->p[2,:] 4 p->p[2,:]; 1 p->p[1,:] 3 p->p[1,:]];

# call exasim to generate and run C++ code to solve the PDE model
#sol, pde, mesh,~,~,~,~  = Postprocessing.exasim(pde,mesh);

# search compilers and set options
pde = Gencode.setcompilers(pde);

# generate input files and store them in datain folder
pde, mesh, master, dmd = Preprocessing.preprocessing(pde,mesh);

# generate source codes and store them in app folder
Gencode.gencode(pde);

# compile source codes to build an executable file and store it in app folder
compilerstr = Gencode.compilecode(pde);

# run executable file to compute solution and store it in dataout folder
runstr = Gencode.runcode(pde);

# get solution from output files in dataout folder
sol = Postprocessing.fetchsolution(pde,master,dmd);


text:
"""
    lcparray

Kasai's algorithm for linear-time construction of LCP array from Suffix Array
"""
function lcparray(sa::Vector{Int}, data::Vector{T}) where T<:Integer
    n = length(sa)
    lcps = similar(sa)
    rank = similar(sa)
    for i = 1:n
        rank[sa[i]] = i
    end

    lcp = 0
    for i = 1:n
        rank[i] > 1 || continue
        j = sa[rank[i]-1]
        while i+lcp <= n && j+lcp <= n && data[i+lcp] == data[j+lcp]
            lcp += 1
        end
        lcps[rank[i]] = lcp
        lcp > 0 && (lcp -= 1)
    end
    lcps[1] = 0
    lcps
end


text:
include("runtests_preamble.jl")

using LESbrary.TurbulenceStatistics:
    horizontally_averaged_tracers,
    velocity_covariances,
    tracer_covariances,
    third_order_velocity_statistics,
    third_order_tracer_statistics,
    first_order_statistics,
    second_order_statistics,
    third_order_statistics,
    first_through_second_order,
    first_through_third_order,
    subfilter_momentum_fluxes,
    subfilter_tracer_fluxes,
    turbulent_kinetic_energy_budget

using Oceananigans.Fields: XYReducedField

function output_works(simulation, output, output_name="")
    model = simulation.model
    model.clock.time = 0
    model.clock.iteration = 0
    simulation.stop_iteration = 1

    simulation.output_writers[:test] = JLD2OutputWriter(model, output,
                                                        schedule = IterationInterval(1),
                                                        prefix = "test",
                                                        dir = ".")

    success = try
        run!(simulation)
        true
    catch err
        @warn "Output test for $output_name failed with " * sprint(showerror, err)
        false
    finally
        rm("test.jld2")
        pop!(simulation.output_writers, :test)
    end

    return success
end

for arch in architectures
    @testset "Turbulence Statistics [$(typeof(arch))]" begin
        @info "Testing turbulence statistics [$(typeof(arch))]..."

        model = NonhydrostaticModel(grid = RectilinearGrid(arch, size=(1, 1, 1), extent=(1, 1, 1)),
                                    tracers = :b,
                                    buoyancy = BuoyancyTracer(),
                                    closure = AnisotropicMinimumDissipation())

        simulation = Simulation(model, Δt=1.0, stop_iteration=1)

        C  = horizontally_averaged_tracers(model)
        u² = velocity_covariances(model)
        c² = tracer_covariances(model)
        u³ = third_order_velocity_statistics(model)
        u³ = third_order_tracer_statistics(model)

        ψ¹ = first_order_statistics(model)
        ψ² = second_order_statistics(model)
        ψ³ = third_order_statistics(model)

        ψ¹_ψ² = first_through_second_order(model)
        ψ¹_ψ³ = first_through_third_order(model)

        Qᵘ = subfilter_momentum_fluxes(model)
        Qᶜ = subfilter_tracer_fluxes(model)

        tke_budget = turbulent_kinetic_energy_budget(model)

        @test all(ϕ isa XYReducedField for ϕ in values( C     ))
        @test all(ϕ isa XYReducedField for ϕ in values( u²    ))
        @test all(ϕ isa XYReducedField for ϕ in values( c²    ))
        @test all(ϕ isa XYReducedField for ϕ in values( u³    ))
        @test all(ϕ isa XYReducedField for ϕ in values( u³    ))

        @test all(ϕ isa XYReducedField for ϕ in values( ψ¹    ))
        @test all(ϕ isa XYReducedField for ϕ in values( ψ²    ))
        @test all(ϕ isa XYReducedField for ϕ in values( ψ³    ))

        @test all(ϕ isa XYReducedField for ϕ in values( ψ¹_ψ² ))
        @test all(ϕ isa XYReducedField for ϕ in values( ψ¹_ψ³ ))

        @test all(ϕ isa XYReducedField for ϕ in values( Qᵘ    ))
        @test all(ϕ isa XYReducedField for ϕ in values( Qᶜ    ))

        @test output_works(simulation, C, "Horizontally averaged tracers")
        @test output_works(simulation, u², "Velocity covariances")
        @test output_works(simulation, c², "Tracer covariances")
        @test output_works(simulation, ψ¹, "First-order statistics")

        # Individually test outputs for better failure messages
        for (name, output) in ψ³
            @test output_works(simulation, Dict(name => output), name)
        end

        for (name, output) in tke_budget
            @test output_works(simulation, Dict(name => output), name)
        end

        for (name, output) in Qᵘ
            @test output_works(simulation, Dict(name => output), name)
        end

        for (name, output) in Qᶜ
            @test output_works(simulation, Dict(name => output), name)
        end
    end
end

max_stars_repo_name: ali-ramadhan/LESbrary

max_issues_repo_name: ali-ramadhan/LESbrary

max_forks_repo_name: ali-ramadhan/LESbrary


text:
module TSne

using LinearAlgebra, Statistics, Distances, ProgressMeter
using Printf: @sprintf

export tsne

"""
Compute the point perplexities `P` given its squared distances to the other points `D`
and the precision of Gaussian distribution `beta`.
"""
function Hbeta!(P::AbstractVector, D::AbstractVector, beta::Number)
    @inbounds P .= exp.(-beta .* D)
    sumP = sum(P)
    @assert (isfinite(sumP) && sumP > 0.0) "Degenerated P: sum=$sumP, beta=$beta"
    H = log(sumP) + beta * dot(D, P) / sumP
    @assert isfinite(H) "Degenerated H"
    @inbounds P .*= 1/sumP
    return H
end

"""
    perplexities(D::AbstractMatrix, tol::Number = 1e-5, perplexity::Number = 30.0;
                 [keyword arguments])

Convert `n×n` squared distances matrix `D` into `n×n` perplexities matrix `P`.
Performs a binary search to get P-values in such a way that each conditional
Gaussian has the same perplexity.
"""
function perplexities(D::AbstractMatrix{T}, tol::Number = 1e-5, perplexity::Number = 30.0;
                      max_iter::Integer = 50,
                      verbose::Bool=false, progress::Bool=true) where T<:Number
    (issymmetric(D) && all(x -> x >= 0, D)) ||
        throw(ArgumentError("Distance matrix D must be symmetric and positive"))

    # initialize
    n = size(D, 1)
    P = fill(zero(T), n, n) # perplexities matrix
    beta = fill(one(T), n)  # vector of Normal distribution precisions for each point
    Htarget = log(perplexity) # the expected entropy
    Di = fill(zero(T), n)
    Pcol = fill(zero(T), n)

    # Loop over all datapoints
    progress && (pb = Progress(n, "Computing point perplexities"))
    for i in 1:n
        progress && update!(pb, i)

        # Compute the Gaussian kernel and entropy for the current precision
        betai = 1.0
        betamin = 0.0
        betamax = Inf

        copyto!(Di, view(D, :, i))
        Di[i] = prevfloat(Inf) # exclude D[i,i] from minimum(), yet make it finite and exp(-D[i,i])==0.0
        minD = minimum(Di) # distance of i-th point to its closest neighbour
        @inbounds Di .-= minD # entropy is invariant to offsetting Di, which helps to avoid overflow

        H = Hbeta!(Pcol, Di, betai)
        Hdiff = H - Htarget

        # Evaluate whether the perplexity is within tolerance
        tries = 0
        while abs(Hdiff) > tol && tries < max_iter
            # If not, increase or decrease precision
            if Hdiff > 0.0
                betamin = betai
                betai = isfinite(betamax) ? (betai + betamax)/2 : betai*2
            else
                betamax = betai
                betai = (betai + betamin)/2
            end

            # Recompute the values
            H = Hbeta!(Pcol, Di, betai)
            Hdiff = H - Htarget
            tries += 1
        end
        verbose && abs(Hdiff) > tol && warn("P[$i]: perplexity error is above tolerance: $(Hdiff)")
        # Set the final column of P
        @assert Pcol[i] == 0.0 "Diagonal probability P[$i,$i]=$(Pcol[i]) not zero"
        @inbounds P[:, i] .= Pcol
        beta[i] = betai
    end
    progress && finish!(pb)
    # Return final P-matrix
    verbose && @info(@sprintf("Mean σ=%.4f", mean(sqrt.(1 ./ beta))))
    return P, beta
end

"""
    pca(X::Matrix, ncols::Integer = 50)

Run PCA on `X` to reduce the number of its dimensions to `ndims`.

FIXME use PCA routine from JuliaStats?
"""
function pca(X::AbstractMatrix, ndims::Integer = 50)
    (n, d) = size(X)
    (d <= ndims) && return X
    Y = X .- mean(X, dims=1)
    C = Symmetric((Y' * Y) ./ (n-1))
    Ceig = eigen(C, (d-ndims+1):d) # take eigvects for top ndims largest eigvals
    return Y * reverse(Ceig.vectors, dims=2)
end

# K-L divergence element
kldivel(p, q) = ifelse(p > zero(p) && q > zero(q), p*log(p/q), zero(p))

# pairwise squared distance
# if X is the matrix of objects, then the distance between its rows
pairwisesqdist(X::AbstractMatrix, dist::Bool) =
    dist ? X.^2 : pairwise(SqEuclidean(), X')

pairwisesqdist(X::AbstractVector, dist::Function) =
    [dist(x, y)^2 for x in X, y in X] # note: some redundant calc since dist should be symmetric

pairwisesqdist(X::AbstractMatrix, dist::Function) =
    [dist(view(X, i, :), view(X, j, :))^2 for i in 1:size(X, 1), j in 1:size(X, 1)] # note: some redundant calc since dist should be symmetric

pairwisesqdist(X::AbstractMatrix, dist::SemiMetric) =
    pairwise(dist, X').^2 # use Distances

"""
    tsne(X::Union{AbstractMatrix, AbstractVector}, ndims::Integer=2, reduce_dims::Integer=0,
         max_iter::Integer=1000, perplexity::Number=30.0; [keyword arguments])

Apply t-SNE (t-Distributed Stochastic Neighbor Embedding) to `X`,
i.e. embed its points (rows) into `ndims` dimensions preserving close neighbours.

Returns the points×`ndims` matrix of calculated embedded coordinates.

Different from the orginal implementation,
the default is not to use PCA for initialization.

### Arguments
* `distance` if `true`, specifies that `X` is a distance matrix,
  if of type `Function` or `Distances.SemiMetric`, specifies the function to
  use for calculating the distances between the rows
  (or elements, if `X` is a vector) of `X`
* `reduce_dims` the number of the first dimensions of `X` PCA to use for t-SNE,
  if 0, all available dimension are used
* `pca_init` whether to use the first `ndims` of `X` PCA as the initial t-SNE layout,
  if `false` (the default), the method is initialized with the random layout
* `max_iter` how many iterations of t-SNE to do
* `perplexity` the number of "effective neighbours" of a datapoint,
  typical values are from 5 to 50, the default is 30
* `verbose` output informational and diagnostic messages
* `progress` display progress meter during t-SNE optimization
* `min_gain`, `eta`, `initial_momentum`, `final_momentum`, `momentum_switch_iter`,
  `stop_cheat_iter`, `cheat_scale` low-level parameters of t-SNE optimization
* `extended_output` if `true`, returns a tuple of embedded coordinates matrix,
  point perplexities and final Kullback-Leibler divergence

See also [Original t-SNE implementation](https://lvdmaaten.github.io/tsne).
"""
function tsne(X::Union{AbstractMatrix, AbstractVector}, ndims::Integer = 2, reduce_dims::Integer = 0,
              max_iter::Integer = 1000, perplexity::Number = 30.0;
              distance::Union{Bool, Function, SemiMetric} = false,
              min_gain::Number = 0.01, eta::Number = 200.0, pca_init::Bool = false,
              initial_momentum::Number = 0.5, final_momentum::Number = 0.8, momentum_switch_iter::Integer = 250,
              stop_cheat_iter::Integer = 250, cheat_scale::Number = 12.0,
              verbose::Bool = false, progress::Bool=true,
              extended_output = false)
    # preprocess X
    ini_Y_with_X = false
    if isa(X, AbstractMatrix) && (distance !== true)
        verbose && @info("Initial X shape is $(size(X))")
        ndims < size(X, 2) || throw(DimensionMismatch("X has fewer dimensions ($(size(X,2))) than ndims=$ndims"))

        ini_Y_with_X = true
        X = X * (1.0/std(X)::eltype(X)) # note that X is copied
        if 0<reduce_dims<size(X, 2)
            reduce_dims = max(reduce_dims, ndims)
            verbose && @info("Preprocessing the data using PCA...")
            X = pca(X, reduce_dims)
        end
    end
    n = size(X, 1)
    # Initialize embedding
    if pca_init && ini_Y_with_X
        verbose && @info("Using the first $ndims components of the data PCA as the initial layout...")
        if reduce_dims >= ndims
            Y = X[:, 1:ndims] # reuse X PCA
        else
            @assert reduce_dims <= 0 # no X PCA
            Y = pca(X, ndims)
        end
    else
        verbose && @info("Starting with random layout...")
        Y = randn(n, ndims)
    end

    dY = fill!(similar(Y), 0)     # gradient vector
    iY = fill!(similar(Y), 0)     # momentum vector
    gains = fill!(similar(Y), 1)  # how much momentum is affected by gradient

    # Compute P-values
    verbose && (distance !== true) && @info("Computing pairwise distances...")
    D = pairwisesqdist(X, distance)
    P, beta = perplexities(D, 1e-5, perplexity,
                           verbose=verbose, progress=progress)
    P .+= P' # make P symmetric
    P .*= cheat_scale/sum(P) # normalize + early exaggeration
    sum_P = cheat_scale

    # Run iterations
    progress && (pb = Progress(max_iter, "Computing t-SNE"))
    Q = fill!(similar(P), 0)     # temp upper-tri matrix with 1/(1 + (Y[i]-Y[j])²)
    Ymean = similar(Y, 1, ndims) # average for each embedded dimension
    sum_YY = similar(Y, n, 1)    # square norms of embedded points
    L = fill!(similar(P), 0)     # temp upper-tri matrix for KLdiv gradient calculation
    Lcolsums = similar(L, n, 1)  # sum(Symmetric(L), 2)
    last_kldiv = NaN
    for iter in 1:max_iter
        # Compute pairwise affinities
        BLAS.syrk!('U', 'N', 1.0, Y, 0.0, Q) # Q=YY', updates only the upper tri of Q
        @inbounds for i in 1:size(Q, 2)
            sum_YY[i] = Q[i, i]
        end
        sum_Q = 0.0
        @inbounds for j in 1:size(Q, 2)
            sum_YYj_p1 = 1.0 + sum_YY[j]
            Qj = view(Q, :, j)
            Qj[j] = 0.0
            for i in 1:(j-1)
                sqdist_p1 = sum_YYj_p1 - 2.0 * Qj[i] + sum_YY[i]
                @fastmath Qj[i] = ifelse(sqdist_p1 > 1.0, 1.0 / sqdist_p1, 1.0)
                sum_Q += Qj[i]
            end
        end
        sum_Q *= 2 # the diagonal and lower-tri part of Q is zero

        # Compute the gradient
        inv_sum_Q = 1.0 / sum_Q
        fill!(Lcolsums, 0.0) # column sums
        # fill the upper triangle of L (gradient)
        @inbounds for j in 1:size(L, 2)
            Lj = view(L, :, j)
            Pj = view(P, :, j)
            Qj = view(Q, :, j)
            Lsumj = 0.0
            for i in 1:(j-1)
                @fastmath Lj[i] = l = (Pj[i] - Qj[i]*inv_sum_Q) * Qj[i]
                Lcolsums[i] += l
                Lsumj += l
            end
            Lcolsums[j] += Lsumj
        end
        @inbounds for (i, ldiag) in enumerate(Lcolsums)
            L[i, i] = -ldiag
        end
        # dY = -4LY
        BLAS.symm!('L', 'U', -4.0, L, Y, 0.0, dY)

        # Perform the update
        momentum = iter <= momentum_switch_iter ? initial_momentum : final_momentum
        @inbounds for i in eachindex(gains)
            gains[i] = max(ifelse((dY[i] > 0) == (iY[i] > 0),
                                  gains[i] * 0.8,
                                  gains[i] + 0.2),
                           min_gain)
            iY[i] = momentum * iY[i] - eta * (gains[i] * dY[i])
            Y[i] += iY[i]
        end
        @inbounds Y .-= mean!(Ymean, Y)

        # stop cheating with P-values
        if sum_P != 1.0 && iter >= min(max_iter, stop_cheat_iter)
            P .*= 1/sum_P
            sum_P = 1.0
        end
        # Compute the current value of cost function
        if !isfinite(last_kldiv) || iter == max_iter ||
            (progress && mod(iter, max(max_iter÷20, 10)) == 0)
            local kldiv = 0.0
            @inbounds for j in 1:size(P, 2)
                Pj = view(P, :, j)
                Qj = view(Q, :, j)
                kldiv_j = 0.0
                for i in 1:(j-1)
                    # P and Q are symmetric (only the upper triangle used)
                    @fastmath kldiv_j += kldivel(Pj[i], Qj[i])
                end
                kldiv += 2*kldiv_j + kldivel(Pj[j], Q[j])
            end
            last_kldiv = kldiv/sum_P + log(sum_Q/sum_P) # adjust wrt P and Q scales
        end
        progress && update!(pb, iter,
                            showvalues = Dict(:KL_divergence => @sprintf("%.4f%s", last_kldiv,
                                                                         iter <= stop_cheat_iter ? " (warmup)" : "")))
    end
    progress && (finish!(pb))
    verbose && @info(@sprintf("Final t-SNE KL-divergence=%.4f", last_kldiv))

    # Return solution
    if !extended_output
        return Y
    else
        return Y, beta, last_kldiv
    end
end

end


text:
@testset "abstractbq" begin
    @test_throws ArgumentError BQ.check_kernel_parameters(-1.0, 2)
    @test_throws ArgumentError BQ.check_kernel_parameters([1.0, -1.0], 2)
    @test_throws ArgumentError BQ.check_kernel_parameters(rand(4, 4), 2)
    @test_throws ArgumentError BQ.check_kernel_parameters(1.0, -1)
    @test_nowarn BQ.check_kernel_parameters(1.0, 1.0)

    @test BQ.get_kernel_params(SqExponentialKernel(); a=1, b=2) ==
        (SqExponentialKernel(), (a=1, b=2))
    @test BQ.get_kernel_params(3.0 * SqExponentialKernel(); l=1, σ=2) ==
        (SqExponentialKernel(), (l=1, σ=3.0))
    @test BQ.get_kernel_params(
        3.0 * with_lengthscale(SqExponentialKernel(), 2.0); l=1, σ=2
    ) == (SqExponentialKernel(), (l=2.0, σ=3.0))

    @test_throws ArgumentError BQ.check_transform(SelectTransform([1, 3]))
    @test_nowarn BQ.check_transform(ScaleTransform(2.0))

    k = 3.0 * with_lengthscale(SqExponentialKernel(), 2.0)
    bquad = BayesQuad(k)
    x = rand(10)
    y = rand(10)
    @test BQ.kernel(bquad)(x, y) == k(x, y) # Kernels can not be compared directly
end


text:
function writemime(io::IO, mime::MIME"text/plain", entries::Vector{(Any,Entry)})
    for (m, ent) in entries
        println(io, ">>>")
        println(io, AnsiColor.colorize(:blue, "\n • $(m)\n"; mode = "bold"))
        writemime(io, mime, ent)
    end
end

function writemime(io, ::MIME"text/plain", entry::Entry)
    println(io, entry.docs)
    # print metadata if any
    if !isempty(entry.meta)
        println(io, AnsiColor.colorize(:green, " • Details:\n"))
    end
    for (k, v) in entry.meta
        if isa(v, Vector)
            println(io, "\t ∘ ", k, ":")
            for line in v
                if isa(line, NTuple)
                    println(io, "\t\t", AnsiColor.colorize(:cyan, string(line[1])), ": ", line[2])
                else
                    println(io, "\t\t", string(line))
                end
            end
        else
            println(io, "\t ∘ ", k, ": ", v)
        end
        println(io)
    end
end


text:
using LinearAlgebra, Test
using DynamicPolynomials
using SetProg, SetProg.Sets
const MOI = SetProg.JuMP.MOI

@testset "Spaces" begin
    B = Symmetric([1.0 0.0; 0.0 1.0])
    b = [0.0, 0.0]
    β = 1.0
    h = [0.0, 0.0]
    @polyvar x y z
    dual = Sets.PerspectiveInteriorCone(B, b, β, [z, x, y], h)
    ◯ = Sets.perspective_dual(dual)
    mock = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
    model = JuMP.direct_model(mock);
    Q = [1.0, 0.0, 1.0]
    t = 1.0
    # The difference will be zero hence we put `zeros(6)`
    mock_optimize!(mock) = MOI.Utilities.mock_optimize!(mock, [Q; β; b; zeros(6); t])
    MOI.Utilities.set_mock_optimize!(mock, mock_optimize!)
    @variable(model, ◯◯, Ellipsoid(point=SetProg.InteriorPoint(h)))
    cref = @constraint(model, ◯◯ ⊆ ◯)
    @objective(model, Max, nth_root(volume(◯◯)))

    SetProg.optimize!(model)
    @test JuMP.termination_status(model) == MOI.OPTIMAL
    @test JuMP.objective_sense(model) == MOI.MAX_SENSE
    @test JuMP.objective_value(model) == t
    @test JuMP.value(◯◯).p == ◯.p
end


text:
# Autogenerated wrapper script for Pangolin_jll for x86_64-linux-gnu-cxx11
export libpangolin

using GLEW_jll
using libpng_jll
using JpegTurbo_jll
using GLU_jll
using Xorg_libX11_jll
JLLWrappers.@generate_wrapper_header("Pangolin")
JLLWrappers.@declare_library_product(libpangolin, "libpangolin.so")
function __init__()
    JLLWrappers.@generate_init_header(GLEW_jll, libpng_jll, JpegTurbo_jll, GLU_jll, Xorg_libX11_jll)
    JLLWrappers.@init_library_product(
        libpangolin,
        "lib/libpangolin.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Autogenerated wrapper script for Qt5Location_jll for armv7l-linux-gnueabihf-cxx11
export libqt5location, libqt5positioning, libqt5positioningquick

using Qt5Declarative_jll
JLLWrappers.@generate_wrapper_header("Qt5Location")
JLLWrappers.@declare_library_product(libqt5location, "libQt5Location.so.5")
JLLWrappers.@declare_library_product(libqt5positioning, "libQt5Positioning.so.5")
JLLWrappers.@declare_library_product(libqt5positioningquick, "libQt5PositioningQuick.so.5")
function __init__()
    JLLWrappers.@generate_init_header(Qt5Declarative_jll)
    JLLWrappers.@init_library_product(
        libqt5location,
        "lib/libQt5Location.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5positioning,
        "lib/libQt5Positioning.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5positioningquick,
        "lib/libQt5PositioningQuick.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()