text:
@testset "matching.jl" begin
    @testset "randomMatching(::Graph{Int}, ::Int)::Matching{Int}" begin
        bipartite = Graph(10)
        for i = 1:5
            add_edge!(bipartite, i, i+5)
        end
        matching = SeaPearl.randomMatching(bipartite, 5)
        edgeset = [Edge(m[1], m[2]) for m in matching.matches]

        @test matching.size == 5
        @test length(intersect(edgeset, edges(bipartite))) == 5
    end
    @testset "matchingFromDigraph(::DiGraph{Int}, ::Int)::Matching{Int}" begin
        bipartite = DiGraph(6)
        add_edge!(bipartite, 1, 4)
        add_edge!(bipartite, 5, 1)
        add_edge!(bipartite, 2, 5)
        add_edge!(bipartite, 6, 3)
        add_edge!(bipartite, 6, 2)
        matching = SeaPearl.matchingFromDigraph(bipartite, 3)

        @test matching.size == 2
        @test Pair(1, 4) in matching.matches
        @test Pair(2, 5) in matching.matches
    end
    @testset "augmentMatching!(::DiGraph, ::Int, ::Vector{Int})" begin
        bipartite = DiGraph(6)
        add_edge!(bipartite, 1, 4)
        add_edge!(bipartite, 5, 1)
        add_edge!(bipartite, 2, 5)
        add_edge!(bipartite, 5, 3)
        add_edge!(bipartite, 6, 2)
        free = Set([3])
        res = SeaPearl.augmentMatching!(bipartite, 6, free)

        @test !isnothing(res)
        @test res == Pair(3, 6)
    end
    @testset "buildDigraph!(::DiGraph{Int}, ::Graph{Int}, ::Matching{Int})" begin
        graph = Graph(7)
        digraph = DiGraph(7)
        add_edge!(graph, 1, 4)
        add_edge!(graph, 1, 5)
        add_edge!(graph, 2, 4)
        add_edge!(graph, 2, 7)
        add_edge!(graph, 3, 4)
        add_edge!(graph, 3, 6)
        add_edge!(graph, 3, 7)
        matching = SeaPearl.Matching(2, [Pair(1, 4), Pair(2, 7)])
        SeaPearl.buildDigraph!(digraph, graph, matching)

        @test Edge(1, 4) in edges(digraph)
        @test Edge(5, 1) in edges(digraph)
        @test Edge(4, 2) in edges(digraph)
        @test Edge(2, 7) in edges(digraph)
        @test Edge(4, 3) in edges(digraph)
        @test Edge(6, 3) in edges(digraph)
        @test Edge(7, 3) in edges(digraph)
    end
    @testset "maximizeMatching!(::Graph, ::DiGraph, ::Int)::Matching" begin
    #Replays the example from the paper
        graph = Graph(12)
        digraph = DiGraph(12)
        add_edge!(graph, 1, 7)
        add_edge!(graph, 1, 8)
        add_edge!(graph, 2, 8)
        add_edge!(graph, 2, 9)
        add_edge!(graph, 3, 9)
        add_edge!(graph, 4, 8)
        add_edge!(graph, 4, 10)
        add_edge!(graph, 5, 9)
        add_edge!(graph, 5, 10)
        add_edge!(graph, 5, 11)
        add_edge!(graph, 5, 12)
        add_edge!(graph, 6, 12)
        matching = SeaPearl.Matching(3, [Pair(1, 7), Pair(4, 8), Pair(5, 10)])
        SeaPearl.buildDigraph!(digraph, graph, matching)

        SeaPearl.maximizeMatching!(digraph, 6)
        target = [Edge(1, 7), Edge(2, 8), Edge(3, 9), Edge(4, 10), Edge(5, 11), Edge(6, 12)]
        @test all([e in edges(digraph) for e in target])
        @test all([outdegree(digraph, v) == 1 for v in 1:6])
    end
    @testset "maximumMatching!(::Graph{Int}, ::DiGraph{Int}, ::Int)::Matching{Int}" begin
    # This function and all its dependencies have been tested with an external library
    # Thus the code written at that time is certified to work based on 3000 randomly generated
    # matching cases with optimal matching detection.

    # Warning the testsets aren't as powerful as this external control and in case of a code
    # modification I recommend to run external tests again.

        graph = Graph(13)
        digraph = DiGraph(13)
        add_edge!(graph, 1, 7)
        add_edge!(graph, 1, 8)
        add_edge!(graph, 2, 8)
        add_edge!(graph, 2, 9)
        add_edge!(graph, 3, 7)
        add_edge!(graph, 3, 9)
        add_edge!(graph, 4, 8)
        add_edge!(graph, 4, 10)
        add_edge!(graph, 5, 9)
        add_edge!(graph, 5, 10)
        add_edge!(graph, 5, 11)
        add_edge!(graph, 6, 11)
        add_edge!(graph, 6, 12)
        matching = SeaPearl.maximumMatching!(graph, digraph, 6)

        @test matching.size == 6
        @test all(map(pair -> Edge(pair[1], pair[2]) in edges(graph), matching.matches))
        @test all(map(pair -> Edge(pair[1], pair[2]) in edges(digraph), matching.matches))
        @test all(map(e -> e in edges(graph), edges(digraph)))
        @test ne(graph) == ne(digraph)
    end
end


text:
import GeometryBasics: AbstractPoint, Point, Point2, Point3

export PressureTap, xcoord, ycoord, zcoord, coordinates

struct PressureTap{T} <: AbstractPoint{3,T}
    "Coordinate of the pressure tap"
    point::Point{3,T}
    "Face to which the the pressure tap belongs"
    face::Int
    "Integer id of the pressure tap"
    id::Int
    "Is the pressure tap on the inside surface of the face?"
    isint::Bool
    "Does the pressure tap present problems?"
    hasprob::Bool
end

PressureTap(x,y,z; face=1, id=1, isint=false, hasprob=false) =
    PressureTap(Point3{Float64}(x,y,z),face,id,isint, hasprob)

PressureTap(p::Point{3,T}; face=1, id=1, isint=false, hasprob=false) where {T} =
    PressureTap(p,face,id,isint, hasprob)

#import Base.getindex
#getindex(p::PressureTap, i) = p.position[i]

xcoord(p::PressureTap) = p.point[1]
ycoord(p::PressureTap) = p.point[2]
zcoord(p::PressureTap) = p.point[3]
coordinates(p::PressureTap) = p.point

import Base.convert
convert(::Type{Point{3,T}}, p::PressureTap{T}) where {T} = p.point

ext: jl


text:
module Streams

using Dates
using Base: RefValue

export Stream, TimeStream, lift

unixnow() = datetime2unix(now())

struct Event{T}
    time::Float64
    value::T
end

Event(v) = Event(unixnow(), v)

nullfunction() = nothing

# interace: current, current!, update!
abstract type AbstractStream{T} end

struct Stream{T, F, I} <: AbstractStream{T}
    f::F
    current::RefValue{Event{T}}
    inputs::I
end

Stream(f, e::Event, inputs=()) = Stream(f, Ref(e), inputs)
Stream(e::Event) = Stream(nullfunction, e)
Stream(val) = Stream(Event(val))

current(s::Stream) = s.current[]
current!(s::Stream, val) = s.current[] = val

Base.show(io::IO, s::Stream) = print(io, "Stream with value $(s[])")

function update!(s::Stream)
    current_time = current(s).time
    inputs = map(update!, s.inputs)
    inputs_time = mapreduce(max, inputs, init=-Inf) do str
        current(str).time
    end
    if inputs_time > current_time
        res = apply(s.f, inputs...)
        current!(s, Event(inputs_time, res))
    end
    return s
end

function Base.getindex(s::AbstractStream)
    update!(s)
    return current(s).value
end

Base.setindex!(s::AbstractStream, val) = current!(s, Event(val))

function apply(f, streams::AbstractStream...)
    args = map(getindex, streams)
    return f(args...)
end

function lift(f, streams::AbstractStream...)
    res = apply(f, streams...)
    return Stream(f, Event(res), streams)
end

struct TimeStream{T, F} <: AbstractStream{T}
    f::F
    TimeStream{T}(f::F) where {T, F} = new{T, F}(f)
end

function current(t::TimeStream)
    time = unixnow()
    return Event(time, t.f(time))
end
update!(t::TimeStream) = t

end


text:
# include("gp.jl")
# include("distances.jl")

function get_me_true_check(ℳ, 𝒟::ProfileData)
    # mean error on true check for a single value of γ

    total_error = 0.0
    gpr_prediction = predict(ℳ, 𝒟; postprocessed=true)
    # println("Nt-1$(𝒟.Nt-1)")
    # println("pred$(length(gpr_prediction))")
    # println("vavg$(length(𝒟.vavg))")

    # n = 𝒟.Nt-1
    n = length(gpr_prediction)
    for i in 1:n
        exact    = 𝒟.vavg[i]
        predi    = gpr_prediction[i]
        total_error += euclidean_distance(exact, predi) # euclidean distance
    end

    return total_error / n
end

function get_me_greedy_check(ℳ, 𝒟::ProfileData)
    # mean error on greedy check
    # compares the direct model output to the target for all of the timesteps in the validation set.

    total_error = 0.0
    validation_set = 𝒟.validation_set
    n = length(validation_set)
    for j in 1:n
        test_index = validation_set[j]
        y_prediction = model_output(𝒟.x[test_index], test_index, ℳ, 𝒟)
        error = euclidean_distance(y_prediction, 𝒟.y[test_index])
        total_error += error
    end
    return total_error / n
end

function disparity_vector(ℳ, 𝒟::ProfileData)

    total_error = 0.0
    gpr_prediction = predict(ℳ, 𝒟; postprocessed=true)

    total_error = zeros(length(𝒟.vavg[1])) # D
    n = length(gpr_prediction)
    for i in 1:n
        exact = 𝒟.vavg[i]
        predi = gpr_prediction[i]
        total_error .+= exact .- predi # euclidean distance
    end

    return abs.(total_error / n)
end

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection


text:
struct FakeScalar
end

@testset "polynomials" begin
    @polyvar x y z

    @testset "variables" begin
        @test @inferred(x * x) == @inferred(x^2)
        @test @inferred(x^1) == x
        @test @inferred(x^0) == 1
        @test (@wrappedallocs x * x) == 0
        @test (@wrappedallocs x^2) == 0
        @test (@wrappedallocs x^1) == 0
        @test (@wrappedallocs x^0) == 0
        @test @inferred(zero(x)) == @inferred(zero(typeof(x)))
        @test @inferred(one(x)) == 1
        @test @inferred(zero(x)) == @inferred(0 * x)
        @test nvariables(x) == 1
        @test nvariables(y) == 1
    end

    @testset "orderings" begin
        @test Monomial(y) < Monomial(x)
        @test Monomial(x) > Monomial(y)
        @test Monomial(y) < x
        @test y < Monomial(x)
        @test y < x
        @test x > y
        @test x > Monomial(y)
        @test Monomial(x) > y

        @test x == x
        @test x == Monomial(x)
        @test x == 1x
        @test x != y
        @test x != Monomial(y)
        @test x != 1y
        @test x != nothing
        @test x != 0

        @test x > y
        @test x^1 > y
        @test x^1 > y^1
        @test 1x^1 > 1y^1
        @test x^1 > 1y^1
        @test 1x^1 > y^1
        @test y < x
        @test y < x^1
        @test y^1 < x
        @test 1y^1 < x
        @test 1y^1 < 1x^1
        @test 1x < 2x
        @test 2x > 1x
    end

    @testset "monomials" begin
        m = x^2
        @test isa(m, Monomial)
        @test variables(m) == (x,)
        @test exponents(m) == (2,)
        @test degree(m, x) == 2
        @test degree(m, y) == 0
        @test degree(m, z) == 0
        @test nvariables(m) == 1

        m = x * y
        @test isa(m, Monomial)
        @test variables(m) == (x, y)
        @test exponents(m) == (1, 1)
        @test degree(m, x) == 1
        @test degree(m, y) == 1
        @test degree(m, z) == 0
        @test nvariables(m) == 2

        m = y * x
        @test isa(m, Monomial)
        @test variables(m) == (x, y)
        @test exponents(m) == (1, 1)
        @test degree(m, x) == 1
        @test degree(m, y) == 1
        @test degree(m, z) == 0

        m = y * x^2 * y^3
        @test isa(m, Monomial)
        @test variables(m) == (x, y)
        @test exponents(m) == (2, 4)
        @test degree(m, x) == 2
        @test degree(m, y) == 4

        @test x * y < x^2

        @test @wrappedallocs((x * y) < x^2) == 0

        @test @inferred(zero(m)) == @inferred(zero(typeof(m)))
        @test zero(m) == 0 * m

        @test @inferred(Monomial{tuple(), 0}()) == @inferred(x^0)

        @test @inferred(Monomial{(x, y)}([1, 2])) == Monomial{(x, y)}((1, 2))
    end

    @testset "terms" begin
        t = @inferred 3x
        @test isa(t, Term)
        @test variables(t) == (x,)
        @test coefficient(t) == 3
        @test monomial(t) == Monomial(x)

        t = @inferred x * 3
        @test isa(t, Term)
        @test variables(t) == (x,)
        @test coefficient(t) == 3
        @test monomial(t) == Monomial(x)

        t = @inferred (2x) * y * x
        @test isa(t, Term)
        @test variables(t) == (x, y)
        @test coefficient(t) == 2
        @test exponents(t) == (2, 1)
        @test degree(t, x) == 2
        @test degree(t, y) == 1

        @test @inferred(zero(t)) == @inferred(zero(typeof(t)))
        @test zero(t) == 0 * t
        @test @inferred(one(t)) == @inferred(one(typeof(t)))
        @test 1 == one(t) == 1 * t^0

        m1 = x * y * z^3 * x^2 * y
        m2 = x^3 * y^3 * z^2
        t1 = 2m1
        t2 = 2m2
        @test (@wrappedallocs t1 < t2) == 0

        @test -(2x^2 * y ) == -1 * (2y * x^2)

        @test TypedPolynomials.monomialtype(5 * x * y) == Monomial{(x, y), 2}
    end

    @testset "polynomials" begin
        p = @inferred x + 1
        @test isa(p, Polynomial)
        @test p.terms[1] == Term(x)
        @test variables(p) == (x,)

        @test @inferred(zero(p)) == @inferred(zero(typeof(p)))
        @test zero(p) == 0 * p
        @test @inferred(one(p)) == @inferred(one(typeof(p)))
        @test one(p) == 1

        @test mindegree(p) == 0
        @test maxdegree(p) == 1
        @test extdegree(p) == (0, 1)
        @test nvariables(p) == 1

        p = @inferred x^2 + x * x
        @test isa(p, Polynomial)
        @test length(terms(p)) == 1
        @test coefficient(terms(p)[1]) == 2

        @test mindegree(p) == 2
        @test maxdegree(p) == 2
        @test extdegree(p) == (2, 2)
        @test nvariables(p) == 1

        p = @inferred (1 + x) + (y + 2)
        @test coefficient(terms(p)[3]) == 3
        @test exponents(terms(p)[3]) == (0, 0)
        @test variables(p) == (x, y)
        @test coefficient(terms(p)[2]) == 1
        @test exponents(terms(p)[2]) == (0, 1)
        @test coefficient(terms(p)[1]) == 1
        @test exponents(terms(p)[1]) == (1, 0)

        @test mindegree(p) == 0
        @test maxdegree(p) == 1
        @test extdegree(p) == (0, 1)
        @test nvariables(p) == 2

        p = @inferred x^2 + y + x * x + 3 * x * y + x * y^2
        @test length(terms(p)) == 4
        @test coefficient.(terms(p)) == [1, 2, 3, 1]
        @test exponents.(terms(p)) == [(1, 2), (2, 0), (1, 1), (0, 1)]

        @test mindegree(p) == 1
        @test maxdegree(p) == 3
        @test extdegree(p) == (1, 3)
        @test nvariables(p) == 2

        @test (@wrappedallocs x^2 + y + x * x + 3 * x * y + x * y) <= 1376
        @test (@wrappedallocs x^2 + 1) <= 384

        @test (1 + x) * (x + 3) == 3 + 4x + x^2
        @test (2.0 + x) * (y + 1) == 2 + 2y + x + x * y
        @test (x + 1) - 1 == x

        @test @inferred(Polynomial(x) * x^2) == x^3
        @test @inferred((x + y) * x^2) == @inferred(x^3 + x^2 * y)

        @test @inferred(Polynomial([x, y])) == y + x
    end

    @testset "equality" begin
        @test 1 + x == 1 + x
        @test 1 + x == x + 1
        @test 1 + (1 + 1e-16) * x ≈ 1 + x
        @test 2 + 3 + x != 3 + x
        @test 2 + 3 + x == 5 + x
        @test x * y + 2 != y
        @test x * y + 2 == x * y + 2 * x^0
        @test x != nothing
        @test x^2 != nothing
        @test 3x != nothing
        @test x + 1 != nothing
        @test nothing != x
        @test nothing != x^2
        @test nothing != 3x
        @test nothing != x + 1

        @test Polynomial([1, 0 * x, y]) == Polynomial([1, y])
        @test Polynomial([1, 1 * x, y]) != Polynomial([1, y])
        @test Polynomial([1, y]) == Polynomial([1, 0 * x, y])
        @test Polynomial([1, y]) != Polynomial([1, 1 * x, y])
    end

    @testset "ordering" begin
        @test x > y
        @test Monomial(x) > Monomial(y)
        @test y < x
        @test Monomial(y) < Monomial(x)
        @test x^2 > x
        @test y^2 > x
        @test y^2 < x^2
        @test 3x < 5x
        @test 3x^2 > 5x
        @test 5x < x^2
        @test 1x^3 > 10y^3
    end

    @testset "linear algebra" begin
        @test @inferred([x, y]' * [1 2; 3 4] * [x, y]) == x^2 + 5 * x * y + 4 * y^2
        @test @inferred([x, y]' * [1 2; 3 4]) == [(x + 3y) (2x + 4y)]
        @test @inferred([x, y]' * [-1, 3]) == 3y - x
    end

    function testmonomials(var, degree)
        [var^i for i in 0:degree]
    end

    @testset "monomial vector" begin
        ms = @inferred testmonomials(x, 3)
        @test eltype(ms) == Monomial{(x,), 1}
    end

    @testset "operators" begin
        @test !iszero(x)
        @test !iszero(x^2)
        @test !iszero(5x)
        @test !iszero(5x + 1)
        @test iszero(0 * x)
        @test iszero(0 * x + 0)

        @test (x)' == x
        @test (x^2)' == x^2
        @test (Term([1 2; 3 4], x))' == Term([1 2; 3 4]', x)

        @test @inferred(dot(1, x)) == 1x
        @test @inferred(dot(x, 5)) == 5x
        @test @inferred(dot(Term([1, 2], x), Term([3, 4], y))) == [1, 2]' * [3, 4] * x * y
    end

    @testset "linear algebra with tuples" begin
        @polyvar x[1:2] y z

        @test x .+ x == 2 .* x
        @test x .- y == (x[1] - y, x[2] - y)

        @test dot(x, x) == x[1]^2 + x[2]^2
        @test dot(x, [y, z]) == x[1]*y + x[2]*z
        @test dot([y, z], x) == x[1]*y + x[2]*z

        xv = @inferred(vec(x))
        @test @inferred(dot(xv, xv)) == x[1]^2 + x[2]^2
        @test @inferred(dot(xv, [y, z])) == x[1]*y + x[2]*z
        @test @inferred(dot([y, z], xv)) == x[1]*y + x[2]*z

        @test @inferred([1 2; 3 4] * xv) == [x[1] + 2x[2], 3x[1] + 4x[2]]
        @test @inferred(dot(xv, [1 2; 3 4] * xv)) == x[1]^2 + 5x[1]*x[2] + 4x[2]^2
    end

    @testset "Term construction shortcut" begin
        @polyvar x y z

        # Verify that our shortcut that (*)(x, m::MonomialLike) = Term(x, m) works
        @test typeof(@inferred(FakeScalar() * x)) == Term{FakeScalar, Monomial{(x,), 1}}
        @test typeof(@inferred(FakeScalar() * x^2)) == Term{FakeScalar, Monomial{(x,), 1}}
        @test typeof(@inferred(x * FakeScalar())) == Term{FakeScalar, Monomial{(x,), 1}}
        @test typeof(@inferred(x^2 * FakeScalar())) == Term{FakeScalar, Monomial{(x,), 1}}
    end
end


text:
module TestUniversalFallback

using Test

import VecMathOptInterface

const MOI = VecMathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

###
### Utilities for running the tests
###

struct UnknownOptimizerAttribute <: MOI.AbstractOptimizerAttribute end

# A few objective/constraint types are supported to test both the fallback and the
# delegation to the internal model
MOI.Utilities.@model(
    ModelForUniversalFallback,
    (),
    (MOI.LessThan,),
    (),
    (),
    (),
    (MOI.ScalarAffineFunction,),
    (),
    ()
)

function MOI.supports(
    ::ModelForUniversalFallback{T},
    ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}},
) where {T}
    return false
end

function MOI.supports_constraint(
    ::ModelForUniversalFallback{T},
    ::Type{MOI.VariableIndex},
    ::Type{
        <:Union{
            MOI.EqualTo{T},
            MOI.GreaterThan{T},
            MOI.Interval{T},
            MOI.Integer,
            MOI.ZeroOne,
        },
    },
) where {T}
    return false
end

###
### The tests
###

function test_MOI_Test()
    inner = ModelForUniversalFallback{Float64}()
    model = MOI.Utilities.UniversalFallback(inner)
    MOI.Test.runtests(
        model,
        MOI.Test.Config(exclude = Any[MOI.optimize!]),
        exclude = String[
            # UniversalFallback fails all these tests because it supports
            # everything...
            "test_attribute_",
            "test_model_supports_constraint_",
            "test_model_copy_to_Unsupported",
            # Bugs in UniversalFallback
            "test_model_LowerBoundAlreadySet",
            "test_model_UpperBoundAlreadySet",
        ],
    )
    return
end

function _test_Optimizer_Model_attributes(uf, model, attr, listattr)
    @test !MOI.supports(model, attr)
    @test MOI.supports(uf, attr)
    @test isempty(MOI.get(uf, listattr))
    MOI.set(uf, attr, 0)
    @test MOI.get(uf, attr) == 0
    @test MOI.get(uf, listattr) == [attr]
    return
end

function _test_Variable_Constraint_attributes(
    uf,
    model,
    attr,
    listattr,
    I::Type{<:MOI.Index},
    addfun,
    x,
    y,
    z,
)
    @test !MOI.supports(model, attr, I)
    @test MOI.supports(uf, attr, I)
    @test isempty(MOI.get(uf, listattr))
    MOI.set(uf, attr, [x, y], [2, 0])
    @test MOI.get(uf, attr, z) === nothing
    @test !MOI.is_empty(uf)
    @test MOI.get(uf, listattr) == [attr]
    MOI.set(uf, attr, z, 5)
    @test MOI.get(uf, attr, y) == 0
    @test MOI.get(uf, attr, [z, x]) == [5, 2]
    @test MOI.get(uf, listattr) == [attr]
    u = addfun(uf)
    @test MOI.get(uf, attr, u) === nothing
    @test MOI.get(uf, listattr) == [attr]
    MOI.set(uf, attr, u, 8)
    @test MOI.get(uf, listattr) == [attr]
    w = addfun(uf)
    @test MOI.get(uf, listattr) == [attr]
    @test MOI.get(uf, attr, w) === nothing
    @test MOI.is_valid(uf, u)
    MOI.delete(uf, u)
    @test !MOI.is_valid(uf, u)
    @test_throws MOI.InvalidIndex{typeof(u)} MOI.delete(uf, u)
    @test MOI.get(uf, listattr) == [attr]
    MOI.set(uf, attr, [w, z], [9, 4])
    @test MOI.get(uf, listattr) == [attr]
    @test MOI.get(uf, attr, w) == 9
    @test MOI.get(uf, attr, x) == 2
    @test MOI.get(uf, attr, z) == 4
    @test MOI.get(uf, attr, y) == 0
    return
end

function test_optimizer_attributes()
    model = ModelForUniversalFallback{Float64}()
    uf = MOI.Utilities.UniversalFallback(model)
    _test_Optimizer_Model_attributes(
        uf,
        model,
        UnknownOptimizerAttribute(),
        MOI.ListOfOptimizerAttributesSet(),
    )
    return
end

function test_model_attributes()
    model = ModelForUniversalFallback{Float64}()
    uf = MOI.Utilities.UniversalFallback(model)
    _test_Optimizer_Model_attributes(
        uf,
        model,
        MOI.Test.UnknownModelAttribute(),
        MOI.ListOfModelAttributesSet(),
    )
    # Test that emptying the uf get's rid of the model attributes!
    @test !MOI.is_empty(uf)
    MOI.empty!(uf)
    @test MOI.is_empty(uf)
    return
end

function test_variable_attributes()
    model = ModelForUniversalFallback{Float64}()
    uf = MOI.Utilities.UniversalFallback(model)
    _test_Variable_Constraint_attributes(
        uf,
        model,
        MOI.Test.UnknownVariableAttribute(),
        MOI.ListOfVariableAttributesSet(),
        MOI.VariableIndex,
        MOI.add_variable,
        MOI.add_variables(uf, 3)...,
    )
    return
end

function test_supported_objective_attributes()
    model = ModelForUniversalFallback{Float64}()
    uf = MOI.Utilities.UniversalFallback(model)
    x, y = MOI.add_variables(uf, 2)
    attr = MOI.ObjectiveFunction{MOI.VariableIndex}()
    @test MOI.supports(model, attr)
    MOI.set(uf, attr, x)
    @test MOI.get(uf, attr) ≈ x
    MOI.set(uf, attr, y)
    @test MOI.get(uf, attr) ≈ y
    return
end

function test_unsupported_objective_attributes()
    model = ModelForUniversalFallback{Float64}()
    uf = MOI.Utilities.UniversalFallback(model)
    x, y = MOI.add_variables(uf, 2)
    attr = MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}()
    @test !MOI.supports(model, attr)
    @test MOI.supports(uf, attr)
    f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0)
    MOI.set(uf, attr, f)
    @test MOI.get(uf, attr) ≈ f
    MOI.modify(uf, attr, MOI.ScalarCoefficientChange(y, 1.0))
    new_obj = MOI.ScalarAffineFunction(
        [MOI.ScalarAffineTerm(1.0, x), MOI.ScalarAffineTerm(1.0, y)],
        0.0,
    )
    @test MOI.get(uf, attr) ≈ new_obj
    @test attr in MOI.get(uf, MOI.ListOfModelAttributesSet())
    @test MOI.get(uf, MOI.ObjectiveFunctionType()) ==
          MOI.ScalarAffineFunction{Float64}
    return
end

function test_supported_constraint_attributes()
    model = ModelForUniversalFallback{Float64}()
    uf = MOI.Utilities.UniversalFallback(model)
    x, y, z = MOI.add_variables(uf, 3)
    function _add_constraint(uf, x::MOI.VariableIndex, ub::Float64)
        return MOI.add_constraint(
            uf,
            MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0),
            MOI.LessThan(ub),
        )
    end
    cx = _add_constraint(uf, x, 0.0)
    cy = _add_constraint(uf, y, 1.0)
    cz = _add_constraint(uf, z, 2.0)
    _test_Variable_Constraint_attributes(
        uf,
        model,
        MOI.Test.UnknownConstraintAttribute(),
        MOI.ListOfConstraintAttributesSet{
            MOI.ScalarAffineFunction{Float64},
            MOI.LessThan{Float64},
        }(),
        MOI.ConstraintIndex{
            MOI.ScalarAffineFunction{Float64},
            MOI.LessThan{Float64},
        },
        uf -> _add_constraint(uf, x, 0.0),
        cx,
        cy,
        cz,
    )
    fy = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, y)], 0.0)
    MOI.set(uf, MOI.ConstraintFunction(), cx, fy)
    @test MOI.get(uf, MOI.ConstraintFunction(), cx) ≈ fy
    @test MOI.get(uf, MOI.CanonicalConstraintFunction(), cx) ≈ fy
    @test MOI.Utilities.is_canonical(
        MOI.get(uf, MOI.CanonicalConstraintFunction(), cx),
    )
    @test MOI.supports(uf, MOI.ConstraintName(), typeof(cx))
    MOI.set(uf, MOI.ConstraintName(), cx, "LessThan")
    @test MOI.get(uf, MOI.ConstraintName(), cx) == "LessThan"
    @test MOI.get(uf, typeof(cx), "LessThan") == cx
    MOI.delete(uf, cx)
    @test MOI.get(uf, typeof(cx), "LessThan") === nothing
    return
end

function test_unsupported_constraint_attributes()
    model = ModelForUniversalFallback{Float64}()
    uf = MOI.Utilities.UniversalFallback(model)
    x, y, z = MOI.add_variables(uf, 3)
    function _add_constraint(uf, x::MOI.VariableIndex, ub::Float64)
        return MOI.add_constraint(
            uf,
            MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0),
            MOI.EqualTo(ub),
        )
    end
    cx = _add_constraint(uf, x, 0.0)
    cy = _add_constraint(uf, y, 1.0)
    cz = _add_constraint(uf, z, 2.0)
    _test_Variable_Constraint_attributes(
        uf,
        model,
        MOI.Test.UnknownConstraintAttribute(),
        MOI.ListOfConstraintAttributesSet{
            MOI.ScalarAffineFunction{Float64},
            MOI.EqualTo{Float64},
        }(),
        MOI.ConstraintIndex{
            MOI.ScalarAffineFunction{Float64},
            MOI.EqualTo{Float64},
        },
        uf -> _add_constraint(uf, x, 0.0),
        cx,
        cy,
        cz,
    )
    fy = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, y)], 0.0)
    MOI.set(uf, MOI.ConstraintFunction(), cx, fy)
    @test MOI.get(uf, MOI.ConstraintFunction(), cx) ≈ fy
    @test MOI.get(uf, MOI.CanonicalConstraintFunction(), cx) ≈ fy
    @test MOI.Utilities.is_canonical(
        MOI.get(uf, MOI.CanonicalConstraintFunction(), cx),
    )
    @test MOI.supports(uf, MOI.ConstraintName(), typeof(cx))
    MOI.set(uf, MOI.ConstraintName(), cx, "EqualTo")
    @test MOI.get(uf, MOI.ConstraintName(), cx) == "EqualTo"
    @test MOI.get(uf, typeof(cx), "EqualTo") == cx
    MOI.delete(uf, cx)
    @test MOI.get(uf, typeof(cx), "EqualTo") === nothing
    return
end

function test_duplicate_names()
    model = ModelForUniversalFallback{Float64}()
    uf = MOI.Utilities.UniversalFallback(model)
    x = MOI.add_variable(uf)
    # These are both intercepted by the fallback.
    f = MOI.ScalarAffineFunction{Float64}([MOI.ScalarAffineTerm(1.0, x)], 0.0)
    x_eq = MOI.add_constraint(uf, f, MOI.EqualTo(1.0))
    x_gt = MOI.add_constraint(uf, f, MOI.GreaterThan(1.0))
    MOI.set(uf, MOI.ConstraintName(), x_eq, "a name")
    MOI.set(uf, MOI.ConstraintName(), x_gt, "a name")
    @test_throws Exception MOI.get(
        uf,
        MOI.ConstraintIndex{
            MOI.ScalarAffineFunction{Float64},
            MOI.EqualTo{Float64},
        },
        "a name",
    )
    @test_throws Exception MOI.get(uf, MOI.ConstraintIndex, "a name")
    return
end

function test_deterministic_constraint_ordering()
    model = ModelForUniversalFallback{Float64}()
    uf = MOI.Utilities.UniversalFallback(model)
    set1 = MOI.EqualTo(1.0)
    set2 = MOI.GreaterThan(1.0)
    F = MOI.ScalarAffineFunction{Float64}
    for sets in [[set1, set2], [set2, set1]]
        model = ModelForUniversalFallback{Float64}()
        uf = MOI.Utilities.UniversalFallback(model)
        x = MOI.add_variable(uf)
        y = MOI.add_variable(uf)
        fx = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0)
        fy = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, y)], 0.0)
        MOI.add_constraint(uf, fx, sets[1])
        MOI.add_constraint(uf, fx, sets[2])
        MOI.add_constraint(uf, fy, sets[1])
        MOI.add_constraint(uf, fy, sets[2])
        # check that the constraint types are in the order they were added in
        @test MOI.get(uf, MOI.ListOfConstraintTypesPresent()) ==
              [(F, typeof(sets[1])), (F, typeof(sets[2]))]
        # check that the constraints given the constraint type are in the order
        # they were added in
        for set in sets
            @test MOI.get(uf, MOI.ListOfConstraintIndices{F,typeof(set)}()) == [
                MOI.ConstraintIndex{F,typeof(set)}(1),
                MOI.ConstraintIndex{F,typeof(set)}(2),
            ]
        end
    end
    return
end

function test_show()
    model = ModelForUniversalFallback{Float64}()
    uf = MOI.Utilities.UniversalFallback(model)
    @test sprint(show, uf) == MOI.Utilities.replace_acronym("""
    $(MOI.Utilities.UniversalFallback{ModelForUniversalFallback{Float64}})
    fallback for $(ModelForUniversalFallback{Float64})""")
    return
end

function test_missing_attribute()
    model = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}())
    @test MOI.get(model, MOI.Test.UnknownModelAttribute()) === nothing
    @test MOI.get(model, MOI.TimeLimitSec()) === nothing
    return
end

end  # module

TestUniversalFallback.runtests()


text:
# ---
# title: 1006. Clumsy Factorial
# id: problem1006
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math
# link: <https://leetcode.com/problems/clumsy-factorial/description/>
# hidden: true
# ---
#
# Normally, the factorial of a positive integer `n` is the product of all
# positive integers less than or equal to `n`.  For example, `factorial(10) = 10
# * 9 * 8 * 7 * 6 * 5 * 4 * 3 * 2 * 1`.
#
# We instead make a _clumsy factorial:_  using the integers in decreasing order,
# we swap out the multiply operations for a fixed rotation of operations:
# multiply (*), divide (/), add (+) and subtract (-) in this order.
#
# For example, `clumsy(10) = 10 * 9 / 8 + 7 - 6 * 5 / 4 + 3 - 2 * 1`.  However,
# these operations are still applied using the usual order of operations of
# arithmetic: we do all multiplication and division steps before any addition or
# subtraction steps, and multiplication and division steps are processed left to
# right.
#
# Additionally, the division that we use is _floor division_  such that `10 * 9
# / 8` equals `11`.  This guarantees the result is an integer.
#
# `Implement the clumsy` function as defined above: given an integer `N`, it
# returns the clumsy factorial of `N`.
#
#
#
# **Example 1:**
#
#
#
#     Input: 4
#     Output:  7
#     Explanation: 7 = 4 * 3 / 2 + 1
#
#
# **Example 2:**
#
#
#
#     Input: 10
#     Output: 12
#     Explanation: 12 = 10 * 9 / 8 + 7 - 6 * 5 / 4 + 3 - 2 * 1
#
#
#
#
# **Note:**
#
#   1. `1 <= N <= 10000`
#   2. `-2^31 <= answer <= 2^31 - 1`  (The answer is guaranteed to fit within a 32-bit integer.)
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
####### SIMULATION OF MARGINAL ###########################################

### +++ INITIALIZATION +++

### import packages
using Random, LinearAlgebra, RDatasets, MixedModels, Distributions, Optim, CSV, Dates

### define functions
vcov2 = function(x, xbar, σ, τ, nj, mj, γ)
    m = length(xbar)
    a = sum(γ / τ^2 .* mj)  # 1-1 entry
    d = (sum(x.^2) - sum(rep(γ .* nj, mj) .* xbar.^2)) / σ^2 # 2-2 entry
    b = sum(rep(γ, mj) .* xbar) / τ^2 # 2-1 entry
    [d -b; -b a] / (a * d - b^2) # return inverse (X^t V^-1 X)^-1
end

beta = function(vcov1, x, xbar, y, ybar, σ, τ, nj, mj, γ)
    c1 = sum(rep(γ, mj) .* ybar) / τ^2 # = 1^t V^-1 y
    c2 = (x' * y - sum(rep(nj .* γ, mj) .* ybar .* xbar)) / σ^2 # = x^t V^-1 y
    vcov1 * [c1 c2]'  # (X^tV^-1X)^-1 X^tV^-1y
end

rep = function(x, y)
    n = length(y)
    out = []
    for i in 1:n
        append!(out, fill(x[i], y[i]))
    end
    out
end

max1 = function(x, y)
    out = zeros(length(x))
    for i in 1:length(x)
        out[i] = maximum([x[i], y])
    end
    out
end

τlist = repeat([4, 6, 8].^0.5, inner = 3)
σlist = repeat([8, 6, 4].^0.5, inner = 3)

mjlist = repeat([[10, 30, 10],[20, 20, 10],[30, 10, 10]], 3)

for j in 1:length(mjlist)

    ###     SELECT PARAMETERS
    mj = mjlist[j]
    τ0 = τlist[j]
    σ0 = σlist[j]

    ### fix parameters
    p = 2
    nsim = 10000
    Δ = [0:0.1:2;]
    Random.seed!(999)

    ### +++ PER SPECIFICATION SETTING +++

    ### to compile
    if !@isdefined(τ0) τ0 = 1 end
    if !@isdefined(σ0) σ0 = 1 end
    if !@isdefined(mj) mj = 1 end

    ### variable parameters
    nj = [1, 10, 50]
    m = sum(mj)
    N = sum(nj .* mj)

    β0 = [-4.9, 0.03] # rand(Uniform(0, 1), p)
    x = 10 * randn(N)
    xbar = by(DataFrame(x = x, v = rep(1:m, rep(nj, mj))), :v, :x => mean)[:,2]

    ### conditional evaluation
    v = τ0 * randn(m)
    u = β0[1] .+ β0[2] * x + rep(v, rep(nj, mj))
    μ0 = β0[1] .+ β0[2] .* xbar + v

    ### obtain conditional variance
    σ1 = zeros(10000, 1)
    τ1 = zeros(10000, 1)
    for i in 1:10000
        data = DataFrame(x = x, v = rep(v, rep(nj, mj)), y = u + σ0 * randn(N))
        model = fit(MixedModel, @formula(y ~ x + (1 | v)), data)
        τ1[i] = VarCorr(model).σρ.v.σ[1]
        σ1[i] = VarCorr(model).s
    end
    σ1 = mean(σ1) # should be (close to) σ0, but theory for conditional convergence is missing
    τ1 = mean(τ1) # should be (close to) τ0, but theory for conditional convergence is missing

    ### +++ KNOWN DELTA +++

    ###     MARGINAL COVARIANCE MATRIX

    ### helpful parameters
    γ0 = τ1^2 ./ (τ1^2 .+ σ1^2 ./ nj)
    vcov0 = vcov2(x, xbar, σ1, τ1, nj, mj, γ0) # = (X^tV^-1X)^-1
    d0 = rep(1 .- γ0, mj) .* hcat(fill(1, m), xbar) # d' = l' - b'X
    dd0 = d0' * d0

    ### assemble covariance matrix parts
    g1 = rep(γ0 .* σ1^2 ./ nj, mj)

    ### assemble inverse of covariance matrix with Woodbury
    Σ0diag = convert(Array{Float64,1}, g1)
    Σ0i = (I - (d0 ./ Σ0diag) *
                inv(inv(vcov0) +  d0' * (d0 ./ Σ0diag)) * d0') * Diagonal(Σ0diag.^-1)

    ###     CONDITIONAL COVARIANCE MATRIX
    l1 = rep(γ0.^2 ./ nj * σ1^2, mj)

    ### assemble inverse of covariance matrix with Woodbury
    cΣ0diag = convert(Array{Float64,1}, l1)
    cΣ0 = Diagonal(cΣ0diag) + d0 * (vcov0 + vcov0 * [1 sum(xbar); sum(xbar) sum(xbar.^2)] *
                vcov0 .* sum(nj .* γ0 .* (1 .- γ0) .* mj) ./ σ1^2) * d0' + (rep(γ0, mj) .* d0) *
                vcov0 * d0' + ((rep(γ0, mj) .* d0) * vcov0 * d0')'
    cΣ0i = inv(cΣ0)

    ###     NONCENTRALITY PARAMETER
    AZv = rep(γ0 .- 1, mj) .* v + rep(nj, mj) .* d0 * vcov0 * d0' * v ./ τ1^2
    λ0 = AZv' * cΣ0i * AZv

    ###     QUANTILES
    qc0 = cquantile(NoncentralChisq(m, λ0), 0.05) # quantile of alpha = 0.05
    qm = cquantile(Chisq(m), 0.05) # quantile of alpha = 0.05

    ### simulation
    Coverage = zeros(5, length(Δ))
    λ = zeros(nsim, 2)
    Volume = ones(nsim, 5)
    for i in 1:nsim

        ### +++ EMPIRICAL +++

        ### variable data
        y = u + σ0 * randn(N) # conditional
        ybar = by(DataFrame(y = y, v = rep(1:m, rep(nj, mj))), :v, :y => mean)[:,2]
        data = DataFrame(x = x, v = rep(v, rep(nj, mj)), y = u + σ0 * randn(N))

        ### fit model
        model = fit(MixedModel, @formula(y ~ x + (1 | v)), data)

        ### make sure the variance parameters are positive
        τ = maximum([VarCorr(model).σρ.v.σ[1], 0.001])
        σ = maximum([VarCorr(model).s, 0.001])

        ###     MARGINAL

        ### helpful parameters
        γ = τ^2 ./ (τ^2 .+ σ^2 ./ nj)
        vcov1 = vcov2(x, xbar, σ, τ, nj, mj, γ) # = (X^tV^-1X)^-1
        d = rep(1 .- γ, mj)  .* hcat(fill(1, m), xbar) # d' = l' - b'X
        #dd = d' * d
        aj = nj .* τ^2 .+ σ^2

        Ivv = sum(rep(nj.^2 ./ aj.^2, mj)) / 2
        Ive = sum(rep(nj ./ aj.^2, mj)) / 2
        Iee = sum(rep((nj .- 1) / σ^2 .+ aj.^-2, mj)) / 2
        Fisher = [Ivv Ive; Ive Iee]
        Vbar = inv(Fisher)

        ### assemble diagonal matrix parts
        g1 = γ .* σ^2 ./ nj
        g3 = zeros(length(nj))
        for l in 1:length(nj)
            g3[l] += sum(Diagonal(nj[l] / aj[l]^4 .*
            [τ^2 * σ^4 .* nj[l] .+ σ^6  -(τ^4 * σ^2 .* nj[l] .+ τ^2 * σ^4);
            -(τ^4 * σ^2 .* nj[l] .+ τ^2 * σ^4) τ^6 .* nj[l] .+ τ^4 * σ^2] * Vbar))
        end

        ### assemble inverse of covariance matrix with Woodbury
        Σdiag = convert(Array{Float64,1}, rep(g1 + 2 * g3, mj))
        Σi = (I - (d ./ Σdiag) * inv(inv(vcov1) + d' * (d ./ Σdiag)) * d') * Diagonal(Σdiag.^-1)

        ###     CONDITIONAL

        ### assemble diagonal matrix parts
        l1 = γ.^2 ./ nj .* σ^2

        ### helpful new parameters
        b = γ ./ nj
        Db = [σ^2; τ^2] / aj.^2

        bvv = - 2 .* nj .* σ^2 ./ aj.^3
        bve = (nj .* τ^2 .- σ^2) ./ aj.^3
        bee = 2 .* τ^2 ./ aj.^3
        D2b = Vector{Vector{Float64}}[ [bvv, bve], [bve, bee]]

        DvIvv = -sum(rep(nj.^3 ./ aj.^3, mj))
        DvIve = sum(rep(nj.^2 ./ aj.^3, mj))
        DvIee = sum(rep(nj / aj.^3, mj))
        DeIvv = -sum(rep(nj.^2 ./ aj.^3, mj))
        DeIve = sum(rep(nj ./ aj.^3, mj))
        DeIee = sum(rep((nj .- 1) / σ^6 + aj.^-3, mj))
        DFisher = Matrix{Float64}[[DvIvv DvIve; DvIve DvIee], [DeIvv DeIve; DeIve DeIee]]

        l3_1 = zeros(length(nj))
        mid = [sum(γ.^2 .* mj) sum(γ.^2 ./ nj .* mj) ] / τ^4
        for l in 1:length(nj)
            for e in 1:2
                l3_1[l] += mid[e] * σ^4 * b[l] * Vbar[e,:]' * Db[:,l]
            end
        end

        l3_2 = zeros(length(nj))
        for l in 1:length(nj)
            for e in 1:2
                for d in 1:2
                    for f in 1:2
                        for g in 1:2
                            l3_2[l] += Vbar[e,f] * Vbar[f,g] * Fisher[e,d] *  b[l] * D2b[e][d][l]
                        end
                    end
                end
            end
        end

        l3_3 = zeros(length(nj))
        for l in 1:length(nj)
            for e in 1:2
                for d in 1:2
                    for f in 1:2
                        l3_3[l] -= 2 * nj[l] * σ^2 * Vbar[e,f] * Fisher[e,d] * b[l] * Vbar[d,:]' * DFisher[e] * Vbar * Db[:,l]
                    end
                end
            end
        end

        l3_4 = zeros(length(nj))
        for l in 1:length(nj)
            for e in 1:2
                for d in 1:2
                    for g in 1:2
                        l3_4[l] += 2 * nj[l] * σ^2 * DFisher[g][e,d] * Vbar[e,d] * b[l] * Vbar[e,:]' * Db[:,l]
                    end
                end
            end
        end
        l3 = l3_1 .+ l3_2 .+ l3_3 .+ l3_4

        l4 = zeros(length(nj))
        for l in 1:length(nj)
            l4[l] += sum(Diagonal(nj[l] / aj[l]^4 .*
            [σ^6  -τ^2 * σ^4; -τ^2 * σ^4 τ^4 * σ^2] * Vbar))
        end

        l5 = zeros(length(nj))
        for l in 1:length(nj)
            DLvv = σ^4
            DLve = - τ^2 * σ^2
            DLee = τ^4
            l5[l] += nj[l] * (σ^2 - 2 * nj[l] * τ^2) * (nj[l] * τ^2 + σ^2)^-4 * sum(Diagonal([DLvv DLve; DLve DLee] * Vbar)) / 2
        end

        ### assemble inverse of covariance matrix with Woodbury
        cΣdiag = convert(Array{Float64,1}, rep(max1(l1 .+ l3 .+ l4 .- l5, 1/100000), mj))
        cΣ = Diagonal(cΣdiag) + d * (vcov1 + vcov1 * [1 sum(xbar); sum(xbar) sum(xbar.^2)] *
                    vcov1 .* sum(nj .* γ .* (1 .- γ) .* mj) ./ σ1^2) * d' + (rep(γ, mj) .* d) *
                    vcov1 * d' + ((rep(γ, mj) .* d) * vcov1 * d')'
        cΣi = inv(cΣ)

        ###     COMPUTE MEAN
        β1 = beta(vcov0, x, xbar, y, ybar, σ1, τ1, nj, mj, γ0)
        μ1 = rep(1 .- γ0, mj) .* (β1[1] .+ β1[2] .* xbar) + rep(γ0, mj) .* ybar

        β2 = beta(vcov1, x, xbar, y, ybar, σ, τ, nj, mj,  γ) # = fixef(model)
        μ2 = rep(1 .- γ, mj) .* (β2[1] .+ β2[2] .* xbar) + rep(γ, mj) .* ybar

        ###     NONCENTRALITY PARAMETER
        ay = d0 * β1 .- rep(1 .- γ0, mj) .* ybar
        λ[i, 1] = maximum([(ay' * cΣ0i * ay)[1] - tr(rep(γ0.^2 ./ nj.^3, mj) .* Diagonal(cΣ0i)) / τ1^4 * σ1^6, 0])

        ay = d * β1 .- rep(1 .- γ, mj) .* ybar
        λ[i, 2] = maximum([(ay' * cΣi * ay)[1] - tr(rep(γ.^2 ./ nj.^3 ./ aj.^2, mj) .* Diagonal(cΣi)) * σ^6 , 0])

        ###     QUANTILES
        qc1 = cquantile(NoncentralChisq(m, λ[i, 1]), 0.05) # quantile of alpha = 0.05
        qc2 = cquantile(NoncentralChisq(m, λ[i, 2]), 0.05) # quantile of alpha = 0.05

        ###     COMPUTE STATISTICS
        t1 = (μ1 - μ0) .+ Δ'
        t2 = (μ2 - μ0) .+ Δ'
        Coverage[1,:] += (sum(t1' .* (t1' * Σ0i), dims = 2) .< qm) / nsim
        Coverage[2,:] += (sum(t2' .* (t2' * Σi), dims = 2) .< qm) / nsim
        tmp1 = sum(t1' .* (t1' * cΣ0i), dims = 2)
        Coverage[3,:] += (tmp1 .< qc0) / nsim
        Coverage[4,:] += (tmp1 .< qc1) / nsim
        Coverage[5,:] += (sum(t2' .* (t2' * cΣi), dims = 2) .< qc2) / nsim

        ###     COMPUTE RELATIVE VOLUME
        # (via approximations as the determinant of large matrices is numerically unstable)
        Volume[i, 3] = prod(Σ0diag ./ cΣ0diag)^-0.5 * (qm / qc0)^(-m / 2)
        Volume[i, 4] = prod(Σ0diag ./ cΣ0diag)^-0.5 * (qm / qc1)^(-m / 2)
        Volume[i, 5] = prod(Σdiag ./ cΣdiag)^-0.5 * (qm / qc2)^(-m / 2)
    end

    ###     SAVE RESULTS
    filename = "data/unbalanced/" * string(now())
    mkpath(filename)
    Specs = DataFrame(nsim = nsim, m = m, ni = nj, mj = mj, n = N, lambda = λ0,
                        sv = τ0, se = σ0, sCv = τ1, sCe = σ1)
    CSV.write(filename * "/Specs.csv",  Specs, writeheader = true)
    CSV.write(filename * "/Coverage.csv",  DataFrame(Coverage), writeheader = false)
    CSV.write(filename * "/Lambda.csv",  DataFrame(λ), writeheader = false)
    CSV.write(filename * "/Volume.csv",  DataFrame(Volume), writeheader = false)

end


text:
# internal function for easy accessing a view of a column
__!(col1::DatasetColumn) =  col1.val
__!(col1::SubDatasetColumn) =  view(col1.val, col1.selected_index)
const SubOrDSCol = Union{SubDatasetColumn, DatasetColumn}

# we treat DatasetColumn as a one-column data set. and we need to manage every thing ourselves
# we don't encourage people to use ds[!, 1] syntax, manipulating a column of a data set should happen in modify/!
# isequal also use for == , since we don't want missing be annoying
Base.parent(col1::DatasetColumn) = col1.ds

Base.length(col1::SubOrDSCol) = length(__!(col1))
Base.size(col1::SubOrDSCol) = size(__!(col1))
Base.size(col1::SubOrDSCol, i::Integer) = size(__!(col1), i)
Base.isapprox(col1::SubOrDSCol, col2::SubOrDSCol) = isapprox(__!(col1), __!(col2))
Base.first(col1::SubOrDSCol) = first(__!(col1))
Base.last(col1::SubOrDSCol) = last(__!(col1))
Base.eltype(col1::SubOrDSCol) = eltype(__!(col1))
Base.ndims(col1::SubOrDSCol) = ndims(__!(col1))
Base.ndims(::Type{<:SubDatasetColumn}) = 1
Base.isassigned(col1::SubOrDSCol, i) = isassigned(__!(col1), i)
Base.identity(col1::SubOrDSCol) = identity(__!(col1))
Base.similar(col1::SubOrDSCol, args...) = similar(__!(col1), args...)
Base.copy(col1::SubOrDSCol) = copy(__!(col1))
Base.pairs(col1::SubOrDSCol) = pairs(IndexLinear(), __!(col1))
Base.iterate(col1::SubOrDSCol, kwargs...) = iterate(__!(col1), kwargs...)
PooledArrays.PooledArray(col1::SubOrDSCol; arg...) = PooledArray(__!(col1); arg...)
Base.convert(T::Type{<:AbstractVector}, col1::SubOrDSCol) = convert(T, __!(col1))
DataAPI.refarray(col::SubOrDSCol) = DataAPI.refarray(__!(col))
DataAPI.refpool(col::SubOrDSCol) = DataAPI.refpool(__!(col))

Base.isequal(col1::SubOrDSCol, y::Any) = isequal(__!(col1), y)
Base.isequal(y::Any, col1::SubOrDSCol) = isequal(y, __!(col1))
Base.isequal(col1::SubOrDSCol, col2::SubOrDSCol) = isequal(__!(col1), __!(col2))

Base.:(==)(col1::SubOrDSCol, y::Any) = (==)(__!(col1), y)
Base.:(==)(y::Any, col1::SubOrDSCol) = (==)(y, __!(col1))
Base.:(==)(col1::SubOrDSCol, col2::SubOrDSCol) = isequal(__!(col1), __!(col2))
function Base.fill!(col::SubOrDSCol, i)
    fill!(__!(col), i)
    removeformat!(col.ds, col.col)
    col.col ∈ index(parent(col.ds)).sortedcols && _reset_grouping_info!(parent(col.ds))
    _modified(_attributes(parent(col.ds)))
    col
end

function Base.pop!(col::SubOrDSCol)
    res = pop!(__!(col))
    col.col ∈ index(parent(col.ds)).sortedcols && _reset_grouping_info!(parent(col.ds))
    _modified(_attributes(parent(col.ds)))
    res
end
function Base.popfirst!(col::SubOrDSCol)
    res = popfirst!(__!(col))
    col.col ∈ index(parent(col.ds)).sortedcols && _reset_grouping_info!(parent(col.ds))
    _modified(_attributes(parent(col.ds)))
    res
end
function Base.pushfirst!(col::SubOrDSCol, x)
    res = pushfirst!(__!(col), x)
    col.col ∈ index(parent(col.ds)).sortedcols && _reset_grouping_info!(parent(col.ds))
    _modified(_attributes(parent(col.ds)))
    res
end

Base.:(*)(col1::SubOrDSCol, x::Any) = *(__!(col1), x)
Base.:(+)(col1::SubOrDSCol, x::Any) = +(__!(col1), x)
Base.:(/)(col1::SubOrDSCol, x::Any) = /(__!(col1), x)
Base.:(-)(col1::SubOrDSCol, x::Any) = -(__!(col1), x)
Base.:(*)(x::Any, col1::SubOrDSCol) = *(x, __!(col1))
Base.:(+)(x::Any, col1::SubOrDSCol) = +(x, __!(col1))
Base.:(/)(x::Any, col1::SubOrDSCol) = /(x, __!(col1))
Base.:(-)(x::Any, col1::SubOrDSCol) = -(x, __!(col1))
Base.:(*)(col2::SubOrDSCol, col1::SubOrDSCol) = *(__!(col2), __!(col1))
Base.:(+)(col2::SubOrDSCol, col1::SubOrDSCol) = +(__!(col2), __!(col1))
Base.:(/)(col2::SubOrDSCol, col1::SubOrDSCol) = /(__!(col2), __!(col1))
Base.:(-)(col2::SubOrDSCol, col1::SubOrDSCol) = -(__!(col2), __!(col1))
function Base.convert(::Type{T}, x::T) where T<:DatasetColumn
    x
end
function Base.convert(::Type{T}, x::T) where T<:SubDatasetColumn
    x
end

# threads is on for SubOrDSCol since it naturally shouldn't be used for unfavourable situations
Base.maximum(f, col::SubOrDSCol; threads = true) = maximum(f, __!(col), threads = threads)
Base.maximum(col::SubOrDSCol; threads = true) = maximum(identity, __!(col), threads = threads)
Base.minimum(f, col::SubOrDSCol; threads = true) = minimum(f, __!(col), threads = threads)
Base.minimum(col::SubOrDSCol; threads = true) = minimum(identity, __!(col), threads = threads)
Base.sum(f, col::SubOrDSCol; threads = true) = sum(f, __!(col), threads = threads)
Base.sum(col::SubOrDSCol; threads = true) = sum(identity, __!(col), threads = threads)
Statistics.mean(f, col::SubOrDSCol) = mean(f, __!(col))
Statistics.mean(col::SubOrDSCol) = mean(identity, __!(col))
Statistics.var(f, col::SubOrDSCol, dof = true) = var(f, __!(col), dof)
Statistics.var(col::SubOrDSCol, dof = true) = var(identity, __!(col), dof)
Statistics.std(f, col::SubOrDSCol, dof = true) = std(f, __!(col), dof)
Statistics.std(col::SubOrDSCol, dof = true) = std(identity, __!(col), dof)
Statistics.median(col::SubOrDSCol) = median(__!(col))
function Statistics.median!(col::SubOrDSCol)
    median!(__!(col))
    col.col ∈ index(parent(col.ds)).sortedcols && _reset_grouping_info!(parent(col.ds))
    _modified(_attributes(parent(col.ds)))
    col
end
Base.extrema(f, col::SubOrDSCol; threads = true) = extrema(f, __!(col), threads = threads)
Base.extrema(col::SubOrDSCol; threads = true) = extrema(identity, __!(col), threads = threads)
Base.argmax(col::SubOrDSCol; by = identity) = argmax(__!(col), by = by)
Base.argmin(col::SubOrDSCol; by = identity) = argmin(__!(col), by = by)
Base.findmax(f, col::SubOrDSCol) = findmax(f, __!(col))
Base.findmax(col::SubOrDSCol) = findmax(identity, __!(col))
Base.findmin(f, col::SubOrDSCol) = findmin(f, __!(col))
Base.findmin(col::SubOrDSCol) = findmin(identity, __!(col))
Base.cumsum(col::SubOrDSCol; missings = :ignore) = cumsum(__!(col), missings = missings)
Base.cumprod(col::SubOrDSCol; missings = :ignore) = cumprod(__!(col), missings = missings)
cummin(col::SubOrDSCol; missings = :ignore) = cummin(__!(col), missings = missings)
cummax(col::SubOrDSCol; missings = :ignore) = cummax(__!(col), missings = missings)

lag(col::SubOrDSCol; default = missing) = lag(__!(col), default = default)
lag(col::SubOrDSCol, k; default = missing) = lag(__!(col), k, default = default)
lead(col::SubOrDSCol; default = missing) = lead(__!(col), default = default)
lead(col::SubOrDSCol, k; default = missing) = lead(__!(col), k, default = default)

function lag!(col::SubOrDSCol; default = missing)
    lag!(__!(col), default = default)
    _modified(_attributes(parent(col.ds)))
    col.col ∈ index(parent(col.ds)).sortedcols && _reset_grouping_info!(parent(col.ds))
    col
end
function lag!(col::SubOrDSCol, k; default = missing)
    lag!(__!(col), k, default = default)
    _modified(_attributes(parent(col.ds)))
    col.col ∈ index(parent(col.ds)).sortedcols && _reset_grouping_info!(parent(col.ds))
    col
end
function lead!(col::SubOrDSCol; default = missing)
    lead!(__!(col), default = default)
    _modified(_attributes(parent(col.ds)))
    col.col ∈ index(parent(col.ds)).sortedcols && _reset_grouping_info!(parent(col.ds))
    col
end
function lead!(col::SubOrDSCol, k; default = missing)
    lead!(__!(col), k, default = default)
    _modified(_attributes(parent(col.ds)))
    col.col ∈ index(parent(col.ds)).sortedcols && _reset_grouping_info!(parent(col.ds))
    col
end

Base.Sort.defalg(col::SubOrDSCol) = Base.Sort.defalg(__!(col))
function Base.sort!(col::SubOrDSCol; alg::Base.Sort.Algorithm=Base.Sort.defalg(col), lt=isless, by=identity, rev::Bool=false, order::Base.Order.Ordering=Base.Order.Forward)
     sort!(__!(col), alg = alg, lt = lt, by = by, rev = rev, order = order)
     _modified(_attributes(parent(col.ds)))
     col.col ∈ index(parent(col.ds)).sortedcols && _reset_grouping_info!(parent(col.ds))
     col
end
function Base.sort(col::SubOrDSCol; alg::Base.Sort.Algorithm=Base.Sort.defalg(col), lt=isless, by=identity, rev::Bool=false, order::Base.Order.Ordering=Base.Order.Forward)
     sort(__!(col), alg = alg, lt = lt, by = by, rev = rev, order = order)
end

function Base.sortperm(col::SubOrDSCol; alg::Base.Sort.Algorithm=Base.Sort.DEFAULT_UNSTABLE, lt=isless, by=identity, rev::Bool=false, order::Base.Order.Ordering=Base.Order.Forward)
     sortperm(__!(col), alg = alg, lt = lt, by = by, rev = rev, order = order)
end


text:
# defines variables found in INDEXC.INC for use in Julia

const NDDIM = NFLPX + NENGX

const KUVINF =  1
const KUALFA =  2
const KUBETA =  3
const KUROTX =  4
const KUROTY =  5
const KUROTZ =  6
const KUBANK =  7
const KUELEV =  8
const KUHEAD =  9
const KUVIDT = 10
const KUALDT = 11
const KUBEDT = 12
const KUBADT = 13
const KUELDT = 14
const KUHEDT = 15
const KUWXDT = 16
const KUWYDT = 17
const KUWZDT = 18
const KUAXDT = 19
const KUAYDT = 20
const KUAZDT = 21
const KUVACX = 22
const KUVACY = 23
const KUVACZ = 24
const KURACX = 25
const KURACY = 26
const KURACZ = 27
const NUDIM  = 27

const KQVINF =  1
const KQALFA =  2
const KQBETA =  3
const KQROTX =  4
const KQROTY =  5
const KQROTZ =  6
const KQBANK =  7
const KQELEV =  8
const KQHEAD =  9
const KQVACX = 10
const KQVACY = 11
const KQVACZ = 12
const KQRACX = 13
const KQRACY = 14
const KQRACZ = 15
const KQVIDT = 16
const KQALDT = 17
const KQBEDT = 18
const KQBADT = 19
const KQELDT = 20
const KQHEDT = 21
const KQWXDT = 22
const KQWYDT = 23
const KQWZDT = 24
const KQAXDT = 25
const KQAYDT = 26
const KQAZDT = 27
const KQFORC = 28
const KQFORS = 29
const KQFORN = 30
const KQMOMC = 31
const KQMOMS = 32
const KQMOMN = 33
const KQPOSX = 34
const KQPOSY = 35
const KQPOSZ = 36
const KQVELX = 37
const KQVELY = 38
const KQVELZ = 39
const KQGAMA = 40
const NQDIM  = 40


text:
using FlujoFluidos
using Documenter

DocMeta.setdocmeta!(FlujoFluidos, :DocTestSetup, :(using FlujoFluidos); recursive=true)

makedocs(;
    modules=[FlujoFluidos],
    authors="Emilio Alvizo Velázquez",
    repo="https://github.com/EmilioAlvizo/FlujoFluidos.jl/blob/{commit}{path}#{line}",
    sitename="FlujoFluidos.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://EmilioAlvizo.github.io/FlujoFluidos.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/EmilioAlvizo/FlujoFluidos.jl",
)


text:
# test cross-validation
using RatPack
using Random

# first try my (largest current) dataset
file = "../data/nfl_2009.csv"
r0 = 0.0
irule = UpdateIterate( rule=UpdateElo(;r0 = r0, K = 32.0, θ=400.0/log(10.0) ), batch_size=1 )

# just do standard (sampled) ratings
r0 = 0.0
irule = UpdateIterate( rule=UpdateElo(;r0 = r0, K = 32.0, θ=400.0/log(10.0) ), batch_size=1 )
sample_ratings = update_ratings(irule, nfl_ratings, nfl_competitions)

# test prediction/scoring
srule = ScoreBrier(;normalise=true)
score( srule, nfl_competitions, irule, sample_ratings )

# try cross-validation
k = 10
Random.seed!(1)
S1 = cross_validate(irule, nfl_competitions;
                   k = k,
                   n_samples = 1000,
                   batch_size = 1,
                   srule = ScoreBrier(;normalise=true) )

Random.seed!(1)
S2 = cross_validate(irule, file;
                   k = k,
                   n_samples = 1000,
                   batch_size = 1,
                   srule = ScoreBrier(;normalise=true) )

println(" S1bar = $(sum(S1)/k), S2bar = $(sum(S2)/k)")


text:
using DataDeps

register(DataDep(
  "concrete-slump-test-flow",
  "https://archive.ics.uci.edu/ml/datasets/Concrete+Slump+Test",
  "https://archive.ics.uci.edu/ml/machine-learning-databases/concrete/slump/slump_test.data",
  "78cea9fbb28ad8fc17f389e216e724dc04642c6c363143e67f6bcb97e74525ee",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=10,
        feature_indices=2:8,
        id_indices=1,
        skipto=2,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
"""
    ClimaCore.DataLayouts

Notation:
- `i,j` are horizontal node indices within an element
- `k` is the vertical node index within an element
- `f` is the field index
- `v` is the vertical element index in a stack
- `h` is the element stack index

Data layout is specified by the order in which they appear, e.g. `IJKFVH`
indexes the underlying array as `[i,j,k,f,v,h]`

"""
module DataLayouts

import StaticArrays: SOneTo, MArray

# TODO:
#  - doc strings for each type
#  - printing
#  - should some of these be subtypes of AbstractArray?

import ..slab, ..column
export slab, column, IJFH, IJF, IFH, IF, VF, VIFH

include("struct.jl")

abstract type AbstractData{S} end

Base.size(data::AbstractData, i::Integer) = size(data)[i]

function Base.show(io::IO, data::AbstractData)
    indent_width = 2
    (rows, cols) = displaysize(io)
    println(io, summary(data))
    print(io, " "^indent_width)
    print(
        IOContext(
            io,
            :compact => true,
            :limit => true,
            :displaysize => (rows, cols - indent_width),
        ),
        vec(parent(data)),
    )
    return io
end

"""
    DataColumn{S}

Abstract type for data storage for a column. Objects `data` should define a
`data[k,v]`, returning a value of type `S`.
"""
abstract type DataColumn{S} <: AbstractData{S} end

"""
    DataSlab1D{S,Ni}

Abstract type for data storage for a slab of `Ni` values of type `S`.
Objects `data` should define a `data[i]`, returning a value of type `S`.
"""
abstract type DataSlab1D{S, Nij} <: AbstractData{S} end

"""
    DataSlab2D{S,Nij}

Abstract type for data storage for a slab of `Nij × Nij` values of type `S`.
Objects `data` should define a `data[i,j]`, returning a value of type `S`.
"""
abstract type DataSlab2D{S, Nij} <: AbstractData{S} end

"""
    Data1D{S,Ni}

Abstract type for data storage for a 1D field made up of `Ni` values of type `S`.

Objects `data` should define `slab(data, h)` to return a `DataSlab2D{S,Nij}` object.
"""
abstract type Data1D{S, Ni} <: AbstractData{S} end

"""
    Data2D{S,Nij}

Abstract type for data storage for a 2D field made up of `Nij × Nij` values of type `S`.

Objects `data` should define `slab(data, h)` to return a `DataSlab2D{S,Nij}` object.
"""
abstract type Data2D{S, Nij} <: AbstractData{S} end

"""
    Data1DX{S,Ni}

Abstract type for data storage for a 1D field with extruded columns.
The horizontal is made up of `Ni` values of type `S`.

Objects `data` should define `slab(data, v, h)` to return a
`DataSlab1D{S,Ni}` object, and a `column(data, i, h)` to return a `DataColumn`.
"""
abstract type Data1DX{S, Ni} <: AbstractData{S} end

"""
    Data2DX{S,Nij}

Abstract type for data storage for a 2D field with extruded columns.
The horizontal is made  is made up of `Nij × Nij` values of type `S`.

Objects `data` should define `slab(data, v, h)` to return a
`DataSlab2D{S,Nij}` object, and a `column(data, i, j, h)` to return a `DataColumn`.
"""
abstract type Data2DX{S, Nij} <: AbstractData{S} end

"""
    Data3D{S,Nij,Nk}

Abstract type for data storage for a 3D field made up of `Nij × Nij × Nk` values of type `S`.
"""
abstract type Data3D{S, Nij, Nk} <: AbstractData{S} end

Base.eltype(::AbstractData{S}) where {S} = S
@inline function Base.propertynames(::AbstractData{S}) where {S}
    filter(name -> sizeof(fieldtype(S, name)) > 0, fieldnames(S))
end
Base.parent(data::AbstractData) = getfield(data, :array)

Base.similar(data::AbstractData{S}) where {S} = similar(data, S)

function Base.copyto!(dest::D, src::D) where {D <: AbstractData}
    copyto!(parent(dest), parent(src))
    return dest
end

function nfields(data::AbstractData{S}) where {S}
    typesize(eltype(parent(data)), S)
end

# TODO: if this gets used inside kernels, move to a generated function?

@generated function _getproperty(
    data::AbstractData{S},
    ::Val{Name},
) where {S, Name}
    errorstring = "Invalid field name $(Name) for type $(S)"
    i = findfirst(isequal(Name), fieldnames(S))
    if i === nothing
        return :(error($errorstring))
    end
    static_idx = Val{i}()
    return :(Base.@_inline_meta; DataLayouts._property_view(data, $static_idx))
end

@inline function Base.getproperty(data::AbstractData{S}, name::Symbol) where {S}
    _getproperty(data, Val{name}())
end

struct IJKFVH{S, Nij, Nk, A} <: Data3D{S, Nij, Nk}
    array::A
end

function IJKFVH{S, Nij, Nk}(array::AbstractArray{T, 6}) where {S, Nij, Nk, T}
    @assert size(array, 1) == Nij
    @assert size(array, 2) == Nij
    @assert size(array, 3) == Nk
    IJKFVH{S, Nij, Nk, typeof(array)}(array)
end

@generated function _property_view(
    data::IJKFVH{S, Nij, Nk, A},
    idx::Val{Idx},
) where {S, Nij, Nk, A, Idx}
    SS = fieldtype(S, Idx)
    FT = eltype(A)
    offset = fieldtypeoffset(FT, S, Idx)
    nbytes = typesize(FT, SS)
    field_byterange = (offset + 1):(offset + nbytes)
    return :(IJKFVH{$SS, $Nij, $Nk}(
        view(parent(data), :, :, :, $field_byterange, :, :),
    ))
end

@inline function Base.getproperty(
    data::IJKFVH{S, Nij, Nk},
    i::Integer,
) where {S, Nij, Nk}
    array = parent(data)
    T = eltype(array)
    SS = fieldtype(S, i)
    offset = fieldtypeoffset(T, S, i)
    len = typesize(T, SS)
    IJKFVH{SS, Nij, Nk}(view(array, :, :, :, (offset + 1):(offset + len), :, :))
end

function Base.size(data::IJKFVH{S, Nij, Nk}) where {S, Nij, Nk}
    Nv = size(parent(data), 5)
    Nh = size(parent(data), 6)
    return (Nij, Nij, Nk, Nv, Nh)
end

struct IJFH{S, Nij, A} <: Data2D{S, Nij}
    array::A
end

function IJFH{S, Nij}(array::AbstractArray{T, 4}) where {S, Nij, T}
    @assert size(array, 1) == Nij
    @assert size(array, 2) == Nij
    IJFH{S, Nij, typeof(array)}(array)
end

rebuild(data::IJFH{S, Nij}, array) where {S, Nij} = IJFH{S, Nij}(array)
Base.copy(data::IJFH{S, Nij}) where {S, Nij} = IJFH{S, Nij}(copy(parent(data)))
function Base.size(data::IJFH{S, Nij}) where {S, Nij}
    Nv = 1
    Nh = size(parent(data), 4)
    (Nij, Nij, 1, Nv, Nh)
end
"""
    IJFH{S,Nij}(ArrayType, nelements)

Construct an IJFH structure given the backing `ArrayType`,
quadrature degrees of freedom `Nij`,
and the number of mesh elements `nelements`.
"""
function IJFH{S, Nij}(ArrayType, nelements) where {S, Nij}
    FT = eltype(ArrayType)
    IJFH{S, Nij}(ArrayType(undef, Nij, Nij, typesize(FT, S), nelements))
end

Base.length(data::IJFH) = size(parent(data), 4)

@generated function _property_view(
    data::IJFH{S, Nij, A},
    idx::Val{Idx},
) where {S, Nij, A, Idx}
    SS = fieldtype(S, Idx)
    FT = eltype(A)
    offset = fieldtypeoffset(FT, S, Idx)
    nbytes = typesize(FT, SS)
    field_byterange = (offset + 1):(offset + nbytes)
    return :(IJFH{$SS, $Nij}(view(parent(data), :, :, $field_byterange, :)))
end

@inline function Base.getproperty(data::IJFH{S, Nij}, i::Integer) where {S, Nij}
    array = parent(data)
    T = eltype(array)
    SS = fieldtype(S, i)
    offset = fieldtypeoffset(T, S, i)
    len = typesize(T, SS)
    IJFH{SS, Nij}(view(array, :, :, (offset + 1):(offset + len), :))
end

# 1D Data Layouts
Base.length(data::Data1D) = size(parent(data), 3)

struct IFH{S, Ni, A} <: Data1D{S, Ni}
    array::A
end
function Base.size(data::IFH{S, Ni}) where {S, Ni}
    Nv = 1
    Nh = size(parent(data), 3)
    (Ni, 1, 1, Nv, Nh)
end
function IFH{S, Ni}(array::AbstractArray{T, 3}) where {S, Ni, T}
    @assert size(array, 1) == Ni
    IFH{S, Ni, typeof(array)}(array)
end

rebuild(data::IFH{S, Ni}, array::AbstractArray{T, 3}) where {S, Ni, T} =
    IFH{S, Ni}(array)

Base.copy(data::IFH{S, Ni}) where {S, Ni} = IFH{S, Ni}(copy(parent(data)))

function IFH{S, Ni}(ArrayType, nelements) where {S, Ni}
    FT = eltype(ArrayType)
    IFH{S, Ni}(ArrayType(undef, Ni, typesize(FT, S), nelements))
end

@inline function slab(data::IFH{S, Ni}, h::Integer) where {S, Ni}
    @boundscheck (1 <= h <= length(data)) || throw(BoundsError(data, (h,)))
    IF{S, Ni}(view(parent(data), :, :, h))
end
@inline slab(data::IFH, v::Integer, h::Integer) = slab(data, h)

@generated function _property_view(
    data::IFH{S, Ni, A},
    i::Val{Idx},
) where {S, Ni, A, Idx}
    SS = fieldtype(S, Idx)
    FT = eltype(A)
    offset = fieldtypeoffset(FT, S, Idx)
    nbytes = typesize(FT, SS)
    field_byterange = (offset + 1):(offset + nbytes)
    return :(IFH{$SS, $Ni}(view(parent(data), :, $field_byterange, :)))
end

@inline function Base.getproperty(data::IFH{S, Ni}, f::Integer) where {S, Ni}
    array = parent(data)
    T = eltype(array)
    SS = fieldtype(S, f)
    offset = fieldtypeoffset(T, S, f)
    len = typesize(T, SS)
    IFH{SS, Ni}(view(array, :, (offset + 1):(offset + len), :))
end

"""
    IH1JH2{S, Nij}(data::AbstractMatrix{S})

Stores a 2D field in a matrix using a column-major format.
The primary use is for interpolation to a regular grid.
"""
struct IH1JH2{S, Nij, A} <: Data2D{S, Nij}
    array::A
end
function Base.size(data::IH1JH2{S, Nij}) where {S, Nij}
    Nv = 1
    Nh = div(length(parent(data)), Nij * Nij)
    (Nij, Nij, 1, Nv, Nh)
end

function IH1JH2{S, Nij}(array::AbstractMatrix{S}) where {S, Nij}
    @assert size(array, 1) % Nij == 0
    @assert size(array, 2) % Nij == 0
    IH1JH2{S, Nij, typeof(array)}(array)
end

Base.length(data::IH1JH2{S, Nij}) where {S, Nij} =
    div(length(parent(data)), Nij * Nij)

function Base.similar(
    data::IH1JH2{S, Nij, A},
    ::Type{Eltype},
) where {S, Nij, A, Eltype}
    array = similar(A, Eltype)
    return IH1JH2{Eltype, Nij}(array)
end

Base.copy(data::IH1JH2{S, Nij}) where {S, Nij} =
    IH1JH2{S, Nij}(copy(parent(data)))

@inline function slab(data::IH1JH2{S, Nij}, h::Integer) where {S, Nij}
    N1, N2 = size(parent(data))
    n1 = div(N1, Nij)
    n2 = div(N2, Nij)
    z2, z1 = fldmod(h - 1, n1)

    @boundscheck (1 <= h <= n1 * n2) || throw(BoundsError(data, (h,)))

    return view(parent(data), Nij * z1 .+ (1:Nij), Nij * z2 .+ (1:Nij))
end

#=
struct KFV{S,A} <: DataColumn{S}
  array::A
end
function KFV{S}(array::AbstractArray{T,3}) where {S,T}
  KFV{S,typeof(array)}(array)
end
=#

struct IJF{S, Nij, A} <: DataSlab2D{S, Nij}
    array::A
end

function IJF{S, Nij}(array::AbstractArray{T, 3}) where {S, Nij, T}
    @assert size(array, 1) == Nij
    @assert size(array, 2) == Nij
    IJF{S, Nij, typeof(array)}(array)
end

function Base.size(data::IJF{S, Nij}) where {S, Nij}
    return (Nij, Nij, 1, 1, 1)
end

@generated function _property_view(
    data::IJF{S, Nij, A},
    i::Val{Idx},
) where {S, Nij, A, Idx}
    SS = fieldtype(S, Idx)
    FT = eltype(A)
    offset = fieldtypeoffset(FT, S, Idx)
    nbytes = typesize(FT, SS)
    field_byterange = (offset + 1):(offset + nbytes)
    return :(IJF{$SS, $Nij}(view(parent(data), :, :, $field_byterange)))
end

@inline function Base.getproperty(data::IJF{S, Nij}, i::Integer) where {S, Nij}
    array = parent(data)
    T = eltype(array)
    SS = fieldtype(S, i)
    offset = fieldtypeoffset(T, S, i)
    len = typesize(T, SS)
    IJF{SS, Nij}(view(array, :, :, (offset + 1):(offset + len)))
end

# TODO: should this return a S or a 0-d box containing S?
#  - perhaps the latter, as then it is mutable?

# function column(ijfh::IJFH{S}, i::Integer, j::Integer, h) where {S}
#     get_struct(view(parent(ijfh), i, j, :, h), S)
# end

@inline function Base.getindex(
    ijf::IJF{S, Nij},
    i::Integer,
    j::Integer,
) where {S, Nij}
    @boundscheck (1 <= i <= Nij && 1 <= j <= Nij) ||
                 throw(BoundsError(ijf, (i, j)))
    @inbounds get_struct(view(parent(ijf), i, j, :), S)
end

@inline function Base.setindex!(
    ijf::IJF{S, Nij},
    val,
    i::Integer,
    j::Integer,
) where {S, Nij}
    @boundscheck (1 <= i <= Nij && 1 <= j <= Nij) ||
                 throw(BoundsError(ijf, (i, j)))
    set_struct!(view(parent(ijf), i, j, :), convert(S, val))
end

function Base.size(::DataSlab1D{<:Any, Ni}) where {Ni}
    return (Ni, 1, 1, 1, 1)
end

@propagate_inbounds function Base.getindex(slab::DataSlab1D, I::CartesianIndex)
    slab[I[1]]
end

@propagate_inbounds function Base.setindex!(
    slab::DataSlab1D,
    val,
    I::CartesianIndex,
)
    slab[I[1]] = val
end

Base.lastindex(::DataSlab1D{S, Ni}) where {S, Ni} = Ni

struct IF{S, Ni, A} <: DataSlab1D{S, Ni}
    array::A
end

function IF{S, Ni}(array::AbstractArray{T, 2}) where {S, Ni, T}
    @assert size(array, 1) == Ni
    IF{S, Ni, typeof(array)}(array)
end

@generated function _property_view(
    data::IF{S, Ni},
    idx::Val{Idx},
) where {S, Ni, Idx}
    SS = fieldtype(S, Idx)
    T = basetype(SS)
    offset = fieldtypeoffset(T, S, Idx)
    nbytes = typesize(T, SS)
    field_byterange = (offset + 1):(offset + nbytes)
    return :(IF{$SS, $Ni}(view(parent(data), :, $field_byterange)))
end

@inline function Base.getproperty(data::IF{S, Ni}, f::Integer) where {S, Ni}
    array = parent(data)
    T = eltype(array)
    SS = fieldtype(S, f)
    offset = fieldtypeoffset(T, S, f)
    len = typesize(T, SS)
    IF{SS, Ni}(view(array, :, (offset + 1):(offset + len)))
end

@inline function Base.getindex(data::IF{S, Ni}, i::Integer) where {S, Ni}
    @boundscheck (1 <= i <= Ni) || throw(BoundsError(data, (i,)))
    @inbounds get_struct(view(parent(data), i, :), S)
end

@inline function Base.setindex!(data::IF{S, Ni}, val, i::Integer) where {S, Ni}
    @boundscheck (1 <= i <= Ni) || throw(BoundsError(data, (i,)))
    set_struct!(view(parent(data), i, :), convert(S, val))
end

# TODO: should this return a S or a 0-d box containing S?
#  - perhaps the latter, as then it is mutable?
#function column(ijfh::IJFH{S}, i::Integer, j::Integer, h) where {S}
#    get_struct(view(parent(ijfh), i, j, :, h), S)
#end

@inline function slab(ijfh::IJFH{S, Nij}, h::Integer) where {S, Nij}
    @boundscheck (1 <= h <= length(ijfh)) || throw(BoundsError(ijfh, (h,)))
    IJF{S, Nij}(view(parent(ijfh), :, :, :, h))
end

@inline function slab(ijfh::IJFH{S, Nij}, v::Integer, h::Integer) where {S, Nij}
    @boundscheck (1 <= h <= length(ijfh)) || throw(BoundsError(ijfh, (h,)))
    IJF{S, Nij}(view(parent(ijfh), :, :, :, h))
end

@propagate_inbounds function Base.getindex(
    slab::DataSlab2D{S},
    I::CartesianIndex,
) where {S}
    slab[I[1], I[2]]
end

@propagate_inbounds function Base.setindex!(
    slab::DataSlab2D{S},
    val,
    I::CartesianIndex,
) where {S}
    slab[I[1], I[2]] = val
end

Base.size(::DataSlab2D{S, Nij}) where {S, Nij} = (Nij, Nij, 1, 1, 1)
Base.axes(::DataSlab2D{S, Nij}) where {S, Nij} = (SOneTo(Nij), SOneTo(Nij))

# Data column
Base.length(data::DataColumn) = size(parent(data), 1)
Base.size(data::DataColumn) = (1, 1, 1, length(data), 1)

struct VF{S, A} <: DataColumn{S}
    array::A
end

function VF{S}(array::AbstractArray{T, 2}) where {S, T}
    VF{S, typeof(array)}(array)
end

function VF{S}(array::AbstractVector{T}) where {S, T}
    @assert typesize(T, S) == 1
    VF{S}(reshape(array, (:, 1)))
end

function VF{S}(ArrayType, nelements) where {S}
    FT = eltype(ArrayType)
    VF{S}(ArrayType(undef, nelements, typesize(FT, S)))
end

function replace_basetype(data::VF{S}, ::Type{FT}) where {S, FT}
    SS = replace_basetype(S, FT)
    VF{SS}(similar(parent(data), FT))
end

Base.copy(data::VF{S}) where {S} = VF{S}(copy(parent(data)))
Base.lastindex(data::VF) = length(data)

@generated function _property_view(
    data::VF{S, A},
    idx::Val{Idx},
) where {S, A, Idx}
    SS = fieldtype(S, Idx)
    FT = eltype(A)
    offset = fieldtypeoffset(FT, S, Idx)
    nbytes = typesize(FT, SS)
    field_byterange = (offset + 1):(offset + nbytes)
    return :(VF{$SS}(view(parent(data), :, $field_byterange)))
end

@inline function Base.getproperty(data::VF{S}, i::Integer) where {S}
    array = parent(data)
    T = eltype(array)
    SS = fieldtype(S, i)
    offset = fieldtypeoffset(T, S, i)
    len = typesize(T, SS)
    VF{SS}(view(array, :, (offset + 1):(offset + len)))
end

@propagate_inbounds function Base.getindex(data::VF{S}, i::Integer) where {S}
    get_struct(view(parent(data), i, :), S)
end

@propagate_inbounds function Base.getindex(
    col::DataColumn,
    I::CartesianIndex{5},
)
    col[I[4]]
end

@propagate_inbounds function Base.setindex!(
    col::DataColumn,
    val,
    I::CartesianIndex{5},
)
    col[I[4]] = val
end

@inline function Base.setindex!(data::VF{S}, val, v::Integer) where {S}
    set_struct!(view(parent(data), v, :), convert(S, val))
end

column(data::VF, i, h) = data
column(data::VF, i, j, h) = column(data, i, h)

# combined 1D spectral element + extruded 1D FV column data layout

struct VIFH{S, Ni, A} <: Data1DX{S, Ni}
    array::A
end

function VIFH{S, Ni}(array::AbstractArray{T, 4}) where {S, Ni, T}
    @assert size(array, 2) == Ni
    VIFH{S, Ni, typeof(array)}(array)
end

function Base.size(data::VIFH{<:Any, Ni}) where {Ni}
    Nv = size(parent(data), 1)
    Nh = size(parent(data), 4)
    return (Ni, 1, 1, Nv, Nh)
end

function Base.length(data::VIFH)
    size(parent(data), 1) * size(parent(data), 4)
end
Base.copy(data::VIFH{S, Ni}) where {S, Ni} = VIFH{S, Ni}(copy(parent(data)))

@generated function _property_view(
    data::VIFH{S, Ni},
    idx::Val{Idx},
) where {S, Ni, Idx}
    SS = fieldtype(S, Idx)
    T = basetype(SS)
    offset = fieldtypeoffset(T, S, Idx)
    nbytes = typesize(T, SS)
    field_byterange = (offset + 1):(offset + nbytes)
    return :(VIFH{$SS, $Ni}(view(parent(data), :, :, $field_byterange, :)))
end

@inline function Base.getproperty(data::VIFH{S, Ni}, i::Integer) where {S, Ni}
    array = parent(data)
    T = eltype(array)
    SS = fieldtype(S, i)
    offset = fieldtypeoffset(T, S, i)
    len = typesize(T, SS)
    VIFH{SS, Ni}(view(array, :, :, (offset + 1):(offset + len), :))
end

function slab(data::VIFH{S, Ni}, v, h) where {S, Ni}
    IF{S, Ni}(view(parent(data), v, :, :, h))
end

function column(data::VIFH{S}, i, h) where {S}
    VF{S}(view(parent(data), :, i, :, h))
end

function column(data::VIFH{S}, i, j, h) where {S}
    @assert j == 1
    column(data, i, h)
end

@propagate_inbounds function Base.getindex(data::VIFH, I::CartesianIndex)
    data[I[1], I[4]]
end

@propagate_inbounds function Base.setindex!(data::VIFH, val, I::CartesianIndex)
    data[I[1], I[4]] = val
end

"""
    IV1JH2{S, Ni}(data::AbstractMatrix{S})

Stores values from an extruded 1D spectral field in a matrix using a column-major format.
The primary use is for interpolation to a regular grid.
"""
struct IV1JH2{S, Ni, A} <: Data1DX{S, Ni}
    array::A
end

function IV1JH2{S, Ni}(array::AbstractMatrix{S}) where {S, Ni}
    @assert size(array, 2) % Ni == 0
    IV1JH2{S, Ni, typeof(array)}(array)
end

function Base.size(data::IV1JH2{S, Ni}) where {S, Ni}
    Nv = size(parent(data), 1)
    Nh = div(size(parent(data), 2), Ni)
    (Ni, 1, 1, Nv, Nh)
end

Base.length(data::IV1JH2{S, Ni}) where {S, Ni} = div(length(parent(data)), Ni)

function Base.similar(
    data::IV1JH2{S, Ni, A},
    ::Type{Eltype},
) where {S, Ni, A, Eltype}
    array = similar(A, Eltype)
    return IV1JH2{Eltype, Ni}(array)
end

Base.copy(data::IV1JH2{S, Ni}) where {S, Ni} = IV1JH2{S, Ni}(copy(parent(data)))

@inline function slab(data::IV1JH2{S, Ni}, v::Integer, h::Integer) where {S, Ni}
    N1, N2 = size(parent(data))
    n1 = N1
    n2 = div(N2, Ni)
    _, z2 = fldmod(h - 1, n2)
    @boundscheck (1 <= v <= n1) && (1 <= h <= n2) ||
                 throw(BoundsError(data, (v, h)))
    return view(parent(data), v, Ni * z2 .+ (1:Ni))
end

# combined 2D spectral element + extruded 1D FV column data layout

struct VIJFH{S, Nij, A} <: Data2DX{S, Nij}
    array::A
end

function VIJFH{S, Nij}(array::AbstractArray{T, 5}) where {S, Nij, T}
    @assert size(array, 2) == size(array, 3) == Nij
    VIJFH{S, Nij, typeof(array)}(array)
end

function Base.size(data::VIJFH{<:Any, Nij}) where {Nij}
    Nv = size(parent(data), 1)
    Nh = size(parent(data), 5)
    return (Nij, Nij, 1, Nv, Nh)
end

function Base.length(data::VIJFH)
    size(parent(data), 1) * size(parent(data), 5)
end

function Base.copy(data::VIJFH{S, Nij}) where {S, Nij}
    VIJFH{S, Nij}(copy(parent(data)))
end

@generated function _property_view(
    data::VIJFH{S, Nij},
    idx::Val{Idx},
) where {S, Nij, Idx}
    SS = fieldtype(S, Idx)
    T = basetype(SS)
    offset = fieldtypeoffset(T, S, Idx)
    nbytes = typesize(T, SS)
    field_byterange = (offset + 1):(offset + nbytes)
    return :(VIJFH{$SS, $Nij}(view(parent(data), :, :, :, $field_byterange, :)))
end

@inline function Base.getproperty(
    data::VIJFH{S, Nij},
    i::Integer,
) where {S, Nij}
    array = parent(data)
    T = eltype(array)
    SS = fieldtype(S, i)
    offset = fieldtypeoffset(T, S, i)
    len = typesize(T, SS)
    VIJFH{SS, Nij}(view(array, :, :, :, (offset + 1):(offset + len), :))
end

function slab(data::VIJFH{S, Nij}, v, h) where {S, Nij}
    IJF{S, Nij}(view(parent(data), v, :, :, :, h))
end

function column(data::VIJFH{S}, i, j, h) where {S}
    VF{S}(view(parent(data), :, i, j, :, h))
end

@propagate_inbounds function Base.getindex(data::VIJFH, I::CartesianIndex{5})
    data[I[1], I[2], I[4]]
end

@propagate_inbounds function Base.setindex!(
    data::VIJFH,
    val,
    I::CartesianIndex{5},
)
    data[I[1], I[2], I[4]] = val
end

# broadcast machinery
include("broadcast.jl")

# GPU method specializations
include("cuda.jl")

end # module


text:
#!/usr/bin/env julia

using SimradRaw

function main(args)
    for arg in args
        for datagram in datagrams(arg)
            println(datagram)
        end
    end
end

main(ARGS)


text:
function getCost!(sys::System, pathEnsemble::Ensemble, i::Int)
    for k in eachindex(pathEnsemble)
        cost!(sys, pathEnsemble[k], i)
    end
end

function cost!(sys::System, p::Path, i::Int)
    for ds in eachindex(sys.data)
        for m in sys.model.obs
            @inbounds p.dm[ds] += ((sys.data[ds][i,m+1]-p.x[m]))^2
        end
    end
end

max_stars_repo_name: pzuliani/SPICE

max_issues_repo_name: pzuliani/SPICE

max_forks_repo_name: pzuliani/SPICE


text:
using Test
using MatrixProductStates, TensorNetworkTensors, LinearAlgebra, TensorOperations, KrylovKit
@testset "MatrixProductStates" begin
    include("mps.jl")
    include("dmrg.jl")
    include("tdvp.jl")
end


text:
struct SymbolicIntervalFV{F} <: AbstractSymbolicIntervalBounds{F}
    sym::SymbolicInterval{F}
    lbs::Vector{Vector{Float64}}
    ubs::Vector{Vector{Float64}}
    var_los::Matrix{Float64} # symbolic lower bounds of fresh vars
    var_his::Matrix{Float64} # symbolic upper bounds of fresh vars
    max_vars::Int64
end

function init_symbolic_interval_fv(network::NetworkNegPosIdx, input_set; max_vars=10)
    sym = init_deep_poly_symbolic_interval(input_set)
    lbs = [fill(-Inf, n_nodes(l)) for l in network.layers]
    ubs = [fill( Inf, n_nodes(l)) for l in network.layers]
    # symbolic bounds on variables have number of inputs plus constant
    var_his = zeros(max_vars, dim(input_set) + 1)
    var_los = zeros(max_vars, dim(input_set) + 1)
    return SymbolicIntervalFV(sym, lbs, ubs, var_los, var_his, max_vars)
end

function init_symbolic_interval_fv(s::SymbolicIntervalFV, input_set::AbstractHyperrectangle; max_vars=10)
    sym = init_deep_poly_symbolic_interval(input_set)
    lbs = [copy(lb) for lb in s.lbs]
    ubs = [copy(ub) for ub in s.ubs]
    var_his = zeros(max_vars, dim(input_set) + 1)
    var_los = zeros(max_vars, dim(input_set) + 1)
    return SymbolicIntervalFV(sym, lbs, ubs, var_los, var_his, max_vars)
end

function substitute_variables(sym_lo, sym_hi, var_los, var_his, n_in, n_vars)
    # for sym_lo
    # should we change the mask to [:, n_in + 1: n_in + n_vars] ???
    var_terms = sym_lo[:, n_in + 1: end - 1]

    var_terms⁺ = max.(var_terms, 0)
    var_terms⁻ = min.(var_terms, 0)

    subs_lb = var_terms⁺ * var_los[1:n_vars, :] .+ var_terms⁻ * var_his[1:n_vars, :] .+ sym_lo[:, (1:n_in) ∪ [end]]

    # for sym_hi
    var_terms .= sym_hi[:, n_in + 1: end - 1]

    var_terms⁺ .= max.(var_terms, 0)
    var_terms⁻ .= min.(var_terms, 0)

    subs_ub = var_terms⁺ * var_his[1:n_vars, :] .+ var_terms⁻ * var_los[1:n_vars, :] .+ sym_hi[:, (1:n_in) ∪ [end]]

    return subs_lb, subs_ub
end

# TODO: only calculate upper bound, don't use bounds which calculates lo and up
function LazySets.ρ(d::AbstractArray{T,1} where T, sym::A where A<:SymbolicIntervalFV)
    n_sym = size(sym.sym.Low, 2) - 1
    n_in = dim(domain(sym))
    current_n_vars = n_sym - n_in

    # substitute before interval_map as symbolic intervals are exact for affine
    # transformations and after substitution, dimension of symlo, symhi is smaller
    subs_sym_lo, subs_sym_hi = substitute_variables(sym.sym.Low, sym.sym.Up, sym.var_los, sym.var_his, n_in, current_n_vars)
    tmp_sym = SymbolicInterval(subs_sym_lo, subs_sym_hi, domain(sym))
    d_Low, d_Up = interval_map(reshape(d, 1, length(d)), tmp_sym.Low, tmp_sym.Up)
    sym_prime = SymbolicInterval(d_Low, d_Up, domain(sym))
    lo, up = bounds(sym_prime)
    # the result should only be one-dimensional
    return up[1]
end

######## Splitting

function split_symbolic_interval_bounds(s::SymbolicIntervalFV{<:Hyperrectangle}, index)
    lbs, ubs = low(domain(s)), high(domain(s))
    split_val = 0.5 * (lbs[index] + ubs[index])

    high1 = copy(ubs)
    high1[index] = split_val
    low2 = copy(lbs)
    low2[index] = split_val

    domain1 = Hyperrectangle(low=lbs, high=high1)
    domain2 = Hyperrectangle(low=low2, high=ubs)

    s1 = init_symbolic_interval_fv(s, domain1, max_vars=s.max_vars)
    s2 = init_symbolic_interval_fv(s, domain2, max_vars=s.max_vars)

    return [s1, s2]
end


text:
using Pkg
Pkg.activate("/home/Morgan/FrankWolfe.jl")
# Pkg.activate("/home/Morgan/FrankWolfe.jl/fw-rde")

# parse command line arguments if given
if length(ARGS) > 0
    subdir = ARGS[1]
# otherwise prompt user to specify
else
    print("Please enter sub directory to run RDE in: ")
    subdir = readline()
end
# input validation
while !isdir(joinpath(@__DIR__, subdir))
    print("Invalid directory $subdir. Please enter sub directory to run RDE in: ")
    global subdir = readline()
end

using PyCall
pushfirst!(PyVector(pyimport("sys")["path"]), joinpath(@__DIR__, subdir))

import FrankWolfe
import FrankWolfe: LpNormLMO
include("custom_oralces.jl")
include(joinpath(@__DIR__, subdir, "config.jl"))  # load indices, rates, max_iter
cd(subdir)

# Get the Python side of RDE
rde = pyimport("rde")

for idx in indices

    # Load data sample and distortion functional
    x, fname = rde.get_data_sample(idx)
#     rde.store_test(x, "untargeted_ktest")

    f, df_s, node, pred = rde.get_distortion(x)

    # Wrap objective and gradiet functions
    function func(s, p)
        if !(s isa Vector{eltype(x)})
            s = convert(Vector{eltype(x)}, s)
        end
        if !(p isa Vector{eltype(x)})
            p = convert(Vector{eltype(x)}, p)
        end
        return f(s)
    end

    function grad_s!(storage, s, p)
        if !(s isa Vector{eltype(x)})
            s = convert(Vector{eltype(x)}, s)
        end
        if !(p isa Vector{eltype(x)})
            p = convert(Vector{eltype(x)}, p)
        end
        g = df_s(s)
        return @. storage = g
    end

    function grad_p!(storage, s, p)
        if !(s isa Vector{eltype(x)})
            s = convert(Vector{eltype(x)}, s)
        end
        if !(p isa Vector{eltype(x)})
            p = convert(Vector{eltype(x)}, p)
        end
        g = df_s(p)
        return @. storage = g
    end

    all_s = zeros(eltype(x), (length(rates), length(x)))

    # rates = [50, 100, 150, 200, 250, 300, 350, 400, 1000, 2000]

    for rate in rates
        print("\nrate: ")
        print(rate)
        print("\n")
        # Run FrankWolfe
        println("Running sample $idx with rate $rate")

        s0 = similar(x[:])
        s0 .= 0.0
        lmo_s = NonNegKSparseLMO(rate, 1.0)

        p0 = similar(x[:])
        p0 .= 0.0
        lmo_p = LpNormLMO{Float64,Inf}(d)

        @time s, v_s, p, v_p, primal, dual_gap_s, dual_gap_p  = FrankWolfe.frank_wolfe_2var(
        #@time s, v, primal, dual_gap = FrankWolfe.frank_wolfe(
        #@time s, v, primal, dual_gap = FrankWolfe.away_frank_wolfe(
        #@time s, v, primal, dual_gap = FrankWolfe.blended_conditional_gradient(
        #@time s, v, primal, dual_gap = FrankWolfe.lazified_conditional_gradient(
            (s, p) -> func(s, p),
            (storage, s, p) -> grad_s!(storage, s, p),
            (storage, s, p) -> grad_p!(storage, s, p),
            lmo_s,
            lmo_p,
            s0,
            p0,
            ;fw_arguments...
        )
        # reset adaptive step size if necessary
        if fw_arguments.line_search isa FrankWolfe.MonotonousNonConvexStepSize
            fw_arguments.line_search.factor = 0
        end

        # Store single rate result
        all_s[indexin(rate, rates)[1], :] = s
        # rde.store_single_result(s, idx, fname, rate)
        rde.store_single_result(s, idx, "untargeted_ktest", rate)

        # rde.store_test(s, "untargeted_test")

    end

    break

    # Store multiple rate results
    rde.store_collected_results(all_s, idx, node, pred, fname, rates)

end


text:
@diffusion_process FitzHughNagumoConjug{T} begin
    :dimensions
    process --> 2
    wiener --> 1

    :parameters
    (ϵ, s, γ, β, σ) --> T

    :additional
    constdiff --> true
end

@conjugate_gaussian FitzHughNagumoConjug begin
    :intercept --> (-x[2],)
    :ϵ --> (x[1]-x[1]^3+(1-3*x[1]^2)*x[2],)
    :s --> (one(x[1]),)
    :γ --> (-x[1],)
    :β --> (-one(x[1]),)
    :hypo_a_inv --> 1.0/P.σ^2
    :nonhypo --> 2:2
end

DiffusionDefinition.b(t, x, P::FitzHughNagumoConjug) = @SVector [
    x[2],
    (
        (P.ϵ - P.γ)*x[1]
        - P.ϵ*(x[1]^3 + (3.0*x[1]^2 - 1.0)*x[2])
        + P.s
        - P.β
        - x[2]
    )
]

σ(t, x, P::FitzHughNagumoConjug) = @SMatrix [0.0; P.σ]


text:
function bottleneckAnalysisMaxThroughput(this::PredictionCenter)

	meanAbsError = {};
	meanRelError = {};
	errorHeader = {};
	extra = {};
	title = ""
	legends={}
	Xdata={}
	Ydata={}
	Xlabel=""
	Ylabel=""

	testTPS = this.testSampleTPS;
	testTransactionCount = this.testSampleTransactionCount;

	range = (1:15000)'';
	maxFlushRate = 6000;

	#%cfFlushRateApprox_conf = struct('io_conf', this.testConfig.io_conf, 'workloadName', this.workloadName);
	cfFlushRateApprox_conf = typeOfFlushRate(this.ioConf,"TPCC");

	myFlushRateThroughput = findClosestValue("cfFlushRateApprox", (1:6000)''*this.testMixture, maxFlushRate, cfFlushRateApprox_conf);
	modelP = barzanLinSolve(this.trainConfig.averageCpuUsage, this.trainConfig.transactionCount);

	#%[testMaxThroughputIdx testMaxThroughput] = findMaxThroughput(testTPS);
	trainMaxThroughputIdx, trainMaxThroughput = findMaxThroughput(this.trainConfig.TPSUngrouped);

	actualThr = [];

	idx=[]
	if isdefined(:trainMaxThroughputIdx) && !isempty(trainMaxThroughputIdx)
		idx=(1:trainMaxThroughputIdx)';
	else
		idx=(1:size(this.trainConfig.averageCpuUsage,1))';
	end

	#%CPU-based throughput with classification
	cpuC = barzanLinInvoke(modelP, range*this.testMixture);

	tempLen=length(cpuC)
	cpuCLThroughput=[]
	while tempLen>=1
		if cpuC[tempLen]>88 && cpuC[tempLen]<90
			cpuCLThroughput = tempLen
			break
		end
		tempLen -= 1
	end

	tempLen=length(cpuC)
	cpuCUThroughput=[]
	while tempLen>=1
		if cpuC[tempLen]>98 && cpuC[tempLen]<100
			cpuCUThroughput = tempLen
			break
		end
		tempLen -= 1
	end

	#%CPU-based without classification
	cpuTModel = barzanLinSolve(this.trainConfig.averageCpuUsage, this.trainConfig.TPS);
	cpuT = barzanLinInvoke(cpuTModel, range);

	tempLen=length(cpuT)
	cpuTLThroughput=[]
	while tempLen>=1
		if cpuT[tempLen]>88 && cpuT[tempLen]<90
			cpuTLThroughput = tempLen
			break
		end
		tempLen -= 1
	end

	tempLen=length(cpuT)
	cpuTUThroughput=[]
	while tempLen>=1
		if cpuT[tempLen]>98 && cpuT[tempLen]<100
			cpuTUThroughput = tempLen
			break
		end
		tempLen -= 1
	end
	myModelP = barzanLinSolve(this.trainConfig.averageCpuUsage[vec(idx),:], this.trainConfig.transactionCount[vec(idx),:]);
	myCpuC = barzanLinInvoke(myModelP, range*this.testMixture);

	tempLen=length(myCpuC)
	myCpuCLThroughput=[]
	while tempLen>=1
		if myCpuC[tempLen]>88 && myCpuC[tempLen]<90
			myCpuCLThroughput = tempLen
			break
		end
		tempLen -= 1
	end

	tempLen=length(myCpuC)
	myCpuCUThroughput=[]
	while tempLen>=1
		if myCpuC[tempLen]>98 && myCpuC[tempLen]<100
			myCpuCUThroughput = tempLen
			break
		end
		tempLen -= 1
	end

	tempLen=length(myCpuC)
	myCpuCLThroughput=[]
	while tempLen>=1
		if myCpuC[tempLen]>44 && myCpuC[tempLen]<45
			myCpuCLThroughput = tempLen
			break
		end
		tempLen -= 1
	end

	tempLen=length(myCpuC)
	myCpuCUThroughput=[]
	while tempLen>=1
		if myCpuC[tempLen]>59 && myCpuC[tempLen]<50
			myCpuCUThroughput = tempLen
			break
		end
		tempLen -= 1
	end

	#Our IO-based throughput
	#% cfFlushRateApprox_conf = struct('io_conf', this.testConfig.io_conf, 'workloadName', this.workloadName);
	cfFlushRateApprox_conf = typeOfFlushRate(this.ioConf, "TPCC");
	myFlushRateThroughput = findClosestValue("cfFlushRateApprox", (1:6000)''*this.testMixture, maxFlushRate, cfFlushRateApprox_conf);

	#%Lock-based throughput
	#%getConcurrencyLebel_conf = struct('lock_conf', this.testConfig.lock_conf, 'workloadName', this.workloadName);
	getConcurrencyLebel_conf = typeOfConcurrencyLevel(this.lockConf, "TPCC");
	concurrencyThroughput = findClosestValue("getConcurrencyLevel", (1:10000)''*this.testMixture, 160, getConcurrencyLebel_conf);

	#%Linear IO-based throughput
	modelFlushRate = barzanLinSolve(this.trainConfig.pagesFlushed, this.trainConfig.transactionCount);
	linFlushRate = barzanLinInvoke(modelFlushRate, range*this.testMixture);
	for i=length(linFlushRate):-1:1
		if linFlushRate[i]<maxFlushRate
			linFlushRateThroughput = i
			break
		end
	end
	if isempty(linFlushRateThroughput)
		linFlushRateThroughput=0
	end

	myMaxThroughput1 = min(myCpuCLThroughput, myFlushRateThroughput, concurrencyThroughput);
	if myMaxThroughput1 == myCpuCLThroughput
		PredReasonIdx1 = 1
	elseif myMaxThroughput1 == myFlushRateThroughput
		PredReasonIdx1 = 2
	elseif myMaxThroughput1 == concurrencyThroughput
		PredReasonIdx1 = 3
	end
	myMaxThroughput2 = min(myCpuCUThroughput, myFlushRateThroughput, concurrencyThroughput);
	if myMaxThroughput2 == myCpuCUThroughput
		PredReasonIdx2 = 1
	elseif myMaxThroughput2 == myFlushRateThroughput
		PredReasonIdx2 = 2
	elseif myMaxThroughput2 == concurrencyThroughput
		PredReasonIdx2 = 3
	end

	Xdata = {[1:size(testTPS, 1)]''};
	Ydata = {testTPS};
	legends = {"Signal Generated From User Input"};

	num_row = size(Xdata[1], 1);
	num_col = size(Xdata[1], 2);

	if !isempty(actualThr)
		push!(Xdata,[1:size(testTPS, 1)]'')
		push!(Ydata,repmat([actualThr], num_row, num_col));
		push!(legends,"Actual Max Throughput");
	end
	if !isempty(cpuCLThroughput)
		push!(Xdata,[1:size(testTPS, 1)]'')
		push!(Ydata,repmat([cpuCLThroughput], num_row, num_col));
		push!(legends,"Max Throughput based on adjusted LR for CPU+classification");
	end
	if !isempty(cpuCUThroughput)
		push!(Xdata,[1:size(testTPS, 1)]'')
		push!(Ydata,repmat([cpuCUThroughput], num_row, num_col));
		push!(legends,"Max Throughput on LR for CPU+classification");
	end
	if !isempty(cpuTLThroughput)
		push!(Xdata,[1:size(testTPS, 1)]'')
		push!(Ydata,repmat([cpuTLThroughput], num_row, num_col));
		push!(legends,"Max Throughput on adjusted LR for CPU");
	end
	if !isempty(cpuTUThroughput)
		push!(Xdata,[1:size(testTPS, 1)]'')
		push!(Ydata,repmat([cpuTUThroughput], num_row, num_col));
		push!(legends,"Max Throughput based on LR for CPU");
	end
	if !isempty(myFlushRateThroughput)
		push!(Xdata,[1:size(testTPS, 1)]'')
		push!(Ydata,repmat([myFlushRateThroughput], num_row, num_col));
		push!(legends,"Max Throughput based on our flush rate model");
	end
	if !isempty(linFlushRateThroughput)
		push!(Xdata,[1:size(testTPS, 1)]'')
		push!(Ydata,repmat([linFlushRateThroughput], num_row, num_col));
		push!(legends,"Max Throughput based on LR for flush rate");
	end
	if !isempty(concurrencyThroughput)
		push!(Xdata,[1:size(testTPS, 1)]'')
		push!(Ydata,repmat([concurrencyThroughput], num_row, num_col));
		push!(legends,"Max Throughput based on our contention model");
	end

	finalMaxThroughput = minimum([cpuCLThroughput cpuCUThroughput cpuTLThroughput cpuTUThroughput myFlushRateThroughput linFlushRateThroughput concurrencyThroughput]);

	extra = {[finalMaxThroughput]};

	title = "Bottleneck Analysis: Max Throughput";
	Ylabel = "TPS";
	Xlabel = "Time";

	return title,legends,Xdata,Ydata,Xlabel,Ylabel,meanAbsError,meanRelError,errorHeader,extra
end # end function

max_stars_repo_name: barzan/dbseer

max_issues_repo_name: barzan/dbseer

max_forks_repo_name: barzan/dbseer


text:
import Bibliography
using Test

for file in ["test.bib"]
    test_import = Bibliography.import_bibtex("../examples/$file")
    selection = ["CitekeyArticle", "CitekeyBook"]
    test_select = Bibliography.select(test_import, selection)
    Bibliography.export_web(test_import)
    println("Test import\n $(Bibliography.export_bibtex("result.bib", test_import))\n")
    println("Test select\n $(Bibliography.export_bibtex("result.bib", test_select))\n")
    rm("result.bib")
end