text:
# TODO: uncomment vexity checks once SDP on vars/constraints changes vexity of problem
@testset "SDP Atoms: " for solver in solvers
    if true # can_solve_sdp(solver)
        @testset "sdp variables" begin
            y = Variable((2,2), :Semidefinite)
            p = Problem{BigFloat}(:minimize, y[1,1])
            # @fact vexity(p) --> ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 0 atol=TOL

            y = Variable((3,3), :Semidefinite)
            p = Problem{BigFloat}(:minimize, y[1,1], y[2,2]==1)
            # @fact vexity(p) --> ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 0 atol=TOL

            # Solution is obtained as y[2,2] -> infinity
            # This test fails on Mosek. See
            # https://github.com/JuliaOpt/Mosek.jl/issues/29
            # y = Variable((2, 2), :Semidefinite)
            # p = Problem{BigFloat}(:minimize, y[1, 1], y[1, 2] == 1)
            # # @fact vexity(p) --> ConvexVexity()
            # solve!(p, solver)
            # @fact p.optval --> roughly(0, TOL)

            y = Semidefinite(3)
            p = Problem{BigFloat}(:minimize, sum(diag(y)), y[1, 1] == 1)
            # @fact vexity(p) --> ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 1 atol=TOL

            y = Variable((3, 3), :Semidefinite)
            p = Problem{BigFloat}(:minimize, tr(y), y[2,1]<=4, y[2,2]>=3)
            # @fact vexity(p) --> ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 3 atol=TOL

            x = Variable(Positive())
            y = Semidefinite(3)
            p = Problem{BigFloat}(:minimize, y[1, 2], y[2, 1] == 1)
            # @fact vexity(p) --> ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 1 atol=TOL
        end

        @testset "sdp constraints" begin
            # This test fails on Mosek
            x = Variable(Positive())
            y = Variable((3, 3))
            p = Problem{BigFloat}(:minimize, x + y[1, 1], isposdef(y), x >= 1, y[2, 1] == 1)
            # @fact vexity(p) --> ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 1 atol=TOL
        end

        @testset "nuclear norm atom" begin
            y = Semidefinite(3)
            p = Problem{BigFloat}(:minimize, nuclearnorm(y), y[2,1]<=4, y[2,2]>=3, y[3,3]<=2)
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 3 atol=TOL
            @test evaluate(nuclearnorm(y)) ≈ 3 atol=TOL
        end

        @testset "operator norm atom" begin
            y = Variable((3,3))
            p = Problem{BigFloat}(:minimize, opnorm(y), y[2,1]<=4, y[2,2]>=3, sum(y)>=12)
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 4 atol=TOL
            @test evaluate(opnorm(y)) ≈ 4 atol=TOL
        end

        @testset "sigma max atom" begin
            y = Variable((3,3))
            p = Problem{BigFloat}(:minimize, sigmamax(y), y[2,1]<=4, y[2,2]>=3, sum(y)>=12)
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 4 atol=TOL
            @test evaluate(sigmamax(y)) ≈ 4 atol=TOL
        end

        @testset "lambda max atom" begin
            y = Semidefinite(3)
            p = Problem{BigFloat}(:minimize, lambdamax(y), y[1,1]>=4)
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 4 atol=TOL
            @test evaluate(lambdamax(y)) ≈ 4 atol=TOL
        end

        @testset "lambda min atom" begin
            y = Semidefinite(3)
            p = Problem{BigFloat}(:maximize, lambdamin(y), tr(y)<=6)
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 2 atol=TOL
            @test evaluate(lambdamin(y)) ≈ 2 atol=TOL
        end

        @testset "matrix frac atom" begin
            x = [1, 2, 3]
            P = Variable(3, 3)
            p = Problem{BigFloat}(:minimize, matrixfrac(x, P), P <= 2*eye(3), P >= 0.5 * eye(3))
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 7 atol=TOL
            @test (evaluate(matrixfrac(x, P)))[1] ≈ 7 atol=TOL
        end

        @testset "matrix frac atom both arguments variable" begin
            x = Variable(3)
            P = Variable(3, 3)
            p = Problem{BigFloat}(:minimize, matrixfrac(x, P), lambdamax(P) <= 2, x[1] >= 1)
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 0.5 atol=TOL
            @test (evaluate(matrixfrac(x, P)))[1] ≈ 0.5 atol=TOL
        end

        @testset "sum largest eigs" begin
            x = Semidefinite(3)
            p = Problem{BigFloat}(:minimize, sumlargesteigs(x, 2), x >= 1)
            solve!(p, solver)
            @test p.optval ≈ 3 atol=TOL
            @test evaluate(x) ≈ ones(3, 3) atol=TOL

            x = Semidefinite(3)
            p = Problem{BigFloat}(:minimize, sumlargesteigs(x, 2), [x[i,:] >= i for i=1:3]...)
            solve!(p, solver)
            @test p.optval ≈ 8.4853 atol=TOL

            x1 = Semidefinite(3)
            p1 = Problem{BigFloat}(:minimize, lambdamax(x1), x1[1,1]>=4)
            solve!(p1, solver)

            x2 = Semidefinite(3)
            p2 = Problem{BigFloat}(:minimize, sumlargesteigs(x2, 1), x2[1,1]>=4)
            solve!(p2, solver)

            @test p1.optval ≈ p2.optval atol=TOL

            x1 = Semidefinite(3)
            p1 = Problem{BigFloat}(:minimize, lambdamax(x1), [x1[i,:] >= i for i=1:3]...)
            solve!(p1, solver)

            x2 = Semidefinite(3)
            p2 = Problem{BigFloat}(:minimize, sumlargesteigs(x2, 1), [x2[i,:] >= i for i=1:3]...)
            solve!(p2, solver)

            @test p1.optval ≈ p2.optval atol=TOL

            println(p1.optval)
        end

        @testset "kron atom" begin
            id = eye(4)
            X = Semidefinite(4)
            W = kron(id, X)
            p = Problem{BigFloat}(:maximize, tr(W), tr(X) ≤ 1)
            @test vexity(p) == AffineVexity()
            solve!(p, solver)
            @test p.optval ≈ 4 atol=TOL
        end

        @testset "Partial trace" begin
            A = Semidefinite(2)
            B = [1 0; 0 0]
            ρ = kron(B, A)
            constraints = [partialtrace(ρ, 1, [2; 2]) == BigFloat[0.09942819 0.29923607; 0.29923607 0.90057181], ρ in :SDP]
            p = Problem{BigFloat}(:minimize, Constant(0), constraints)
            solve!(p, solver)
            @test evaluate(ρ) ≈ BigFloat[0.09942819 0.29923607 0 0; 0.299237 0.900572 0 0; 0 0 0 0; 0 0 0 0] atol=TOL
            @test evaluate(partialtrace(ρ, 1, [2; 2])) ≈ [0.09942819 0.29923607; 0.29923607 0.90057181] atol=TOL

            function rand_normalized(n)
                A = 5*randn(n, n) + im*5*randn(n, n)
                A / tr(A)
            end

            As = [ rand_normalized(3) for _ = 1:5]
            Bs = [ rand_normalized(2) for _ = 1:5]
            p = rand(5)

            AB = sum(i -> p[i]*kron(As[i],Bs[i]), 1:5)
            @test partialtrace(AB, 2, [3, 2]) ≈ sum( p .* As )
            @test partialtrace(AB, 1, [3, 2]) ≈ sum( p .* Bs )

            A, B, C = rand(5,5), rand(4,4), rand(3,3)
            ABC = kron(kron(A, B), C)
            @test kron(A,B)*tr(C) ≈ partialtrace(ABC, 3, [5, 4, 3])

            # Test 281
            A = rand(6,6)
            expr = partialtrace(Constant(A), 1, [2, 3])
            @test size(expr) == size(evaluate(expr))

            @test_throws ArgumentError partialtrace(rand(6, 6), 3, [2, 3])
            @test_throws ArgumentError partialtrace(rand(6, 6), 1, [2, 4])
            @test_throws ArgumentError partialtrace(rand(3, 4), 1, [2, 3])
        end

        @testset "Optimization with complex variables" begin
            @testset "Real Variables with complex equality constraints" begin
                n = 10 # variable dimension (parameter)
                m = 5 # number of constraints (parameter)
                xo = rand(n)
                A = randn(m,n) + im*randn(m,n)
                b = A * xo
                x = Variable(n)
                p1 = Problem{BigFloat}(:minimize, sum(x), A*x == b, x>=0)
                solve!(p1, solver)
                x1 = x.value

                p2 = Problem{BigFloat}(:minimize, sum(x), real(A)*x == real(b), imag(A)*x==imag(b), x>=0)
                solve!(p2, solver)
                x2 = x.value
                @test x1 == x2
            end

            @testset "Complex Variable with complex equality constraints" begin
                n = 10 # variable dimension (parameter)
                m = 5 # number of constraints (parameter)
                xo = rand(n)+im*rand(n)
                A = randn(m,n) + im*randn(m,n)
                b = A * xo
                x = ComplexVariable(n)
                p1 = Problem{BigFloat}(:minimize, real(sum(x)), A*x == b, real(x)>=0, imag(x)>=0)
                solve!(p1, solver)
                x1 = x.value

                xr = Variable(n)
                xi = Variable(n)
                p2 = Problem{BigFloat}(:minimize, sum(xr), real(A)*xr-imag(A)*xi == real(b), imag(A)*xr+real(A)*xi == imag(b), xr>=0, xi>=0)
                solve!(p2, solver)
                #x2 = xr.value + im*xi.value
                real_diff = real(x1) - xr.value

                @test real_diff ≈ zeros(10, 1) atol=TOL
                imag_diff = imag(x1) - xi.value
                @test imag_diff ≈ zeros(10, 1) atol=TOL
                #@fact x1==x2 --> true
            end

            @testset "norm2 atom" begin
                a = big"2.0"+big"4.0"*im
                x = ComplexVariable()
                objective = norm2(a-x)
                c1 = real(x)>=0
                p = Problem{BigFloat}(:minimize, objective,c1)
                solve!(p, solver);
                @test p.optval ≈ 0 atol=TOL
                @test evaluate(objective) ≈ 0 atol=TOL
                real_diff = real(x.value) - real(a)
                imag_diff = imag(x.value) - imag(a)
                @test real_diff ≈ 0 atol=TOL
                @test imag_diff ≈ 0 atol=TOL
            end

            @testset "sumsquares atom" begin
                a = [2+4im;4+6im]
                x = ComplexVariable(2)
                objective = sumsquares(a-x)
                c1 = real(x)>=0
                p = Problem{BigFloat}(:minimize, objective,c1)
                solve!(p, solver)
                @test p.optval ≈ 0 atol=TOL
                @test evaluate(objective) ≈ zeros(1, 1) atol=TOL
                real_diff = real.(x.value) - real.(a)
                imag_diff = imag.(x.value) - imag.(a)
                @test real_diff ≈ zeros(2, 1) atol=TOL
                @test imag_diff ≈ zeros(2, 1) atol=TOL
            end

            @testset "abs atom" begin
                a = [5-4im]
                x = ComplexVariable()
                objective = abs(a-x)
                c1 = real(x)>=0
                p = Problem{BigFloat}(:minimize, objective,c1)
                solve!(p, solver)
                @test p.optval ≈ 0 atol=TOL
                @test evaluate(objective) ≈ zeros(1) atol=TOL
                real_diff = real(x.value) .- real(a)
                imag_diff = imag(x.value) .- imag(a)
                @test real_diff ≈ zeros(1) atol=TOL
                @test imag_diff ≈ zeros(1) atol=TOL
            end

            @testset "Complex Semidefinite constraint" begin
                n = 3
                A = rand(n,n) + im*rand(n,n)
                A = A + A' # now A is hermitian
                x = ComplexVariable(n,n)
                objective = sumsquares(A - x)
                c1 = x in :SDP
                p = Problem{BigFloat}(:minimize, objective, c1)
                solve!(p, solver)
                # test that X is approximately equal to posA:
                l,v = eigen(A)
                posA = v*Diagonal(max.(l,0))*v'

                real_diff = real.(x.value) - real.(posA)
                imag_diff = imag.(x.value) - imag.(posA)
                @test real_diff ≈ zeros(n, n) atol=TOL
                @test imag_diff ≈ zeros(n, n) atol=TOL
            end
        end
    end
end


text:
abstract type AlignmentResults end

struct GMAlignmentResult{D,S,T,N,F<:AbstractAffineMap,X<:AbstractGMM{D,S},Y<:AbstractGMM{D,T}} <: AlignmentResults
    x::X
    y::Y
    upperbound::T
    lowerbound::T
    tform::F
    tform_params::NTuple{N,T}
    obj_calls::Int
    num_splits::Int
    num_blocks::Int
    stagnant_splits::Int
end

# Keyword arguments:\n
#     nsplits     - an integer representing the number of splits that should be made along each dimension during branching\n
#     initblock   - a `Block` that defines the searchspace, which defaults to the smallest space gauranteed to contain the global minimum\n
#     rot         - a three-tuple containing a rotation position, which is passed to the `blockfun`\n
#     trl         - a three-tuple containing a translation position, which is passed to the `blockfun`\n
#     blockfun    - the function used for generating `Block`s that define search subspaces (i.e. fullBlock, rotBlock, trlBlock)\n
#     objfun      - the objective function used for local alignment\n
#     atol        - absolute tolerance. Search terminates when the upper bound is within `atol` of the lower bound\n
#     rtol        - relative tolerance. Search terminates when the upper bound is within `rtol*lb` of the lower bound `lb`\n
#     maxblocks   - the maximum number of `Block`s that can be held in the priority queue before search termination\n
#     maxsplits   - the maximum number of `Block` splits that are allowed before search termination\n
#     maxevals    - the maximum number of objective function evaluations allowed before search termination\n
#     maxstagnant - the maximum number of `Block` splits allowed without improvement before search termination\n
#     threads     - when true, utilizes multithreading for performing `Block` splitting. More useful when `nsplits` is large.\n
"""
    result = gogma_align(gmmx, gmmy; nsplits=2, initblock=nothing,
                         rot=nothing, trl=nothing, blockfun=fullBlock, objfun=alignment_objective,
                         rtol=0.01, maxblocks=5e8, maxevals=Inf, maxstagnant=Inf, threads=false)

Finds the globally optimal rigid transform for alignment between two isotropic Gaussian mixtures, `gmmx`
and `gmmy`, using the [GOGMA algorithm](https://arxiv.org/abs/1603.00150).

Returns a `GMAlignmentResult` that contains the maximized overlap of the two GMMs (the upperbound on the objective function),
a lower bound on the alignment objective function, an `AffineMap` which aligns `gmmx` with `gmmy`, and information about the
number of evaluations during the alignment procedure.
"""
function gogma_align(gmmx::AbstractGMM, gmmy::AbstractGMM;
                     nsplits=2, initblock=nothing, rot=nothing, trl=nothing,
                     blockfun=fullBlock, objfun=alignment_objective, tformfun=AffineMap,
                     atol=0.1, rtol=0, maxblocks=5e8, maxsplits=Inf, maxevals=Inf, maxstagnant=Inf, threads=false)
    if isodd(nsplits)
        throw(ArgumentError("`nsplits` must be even"))
    end
    if dims(gmmx) != dims(gmmy)
        throw(ArgumentError("Dimensionality of the GMMs must be equal"))
    end

    # prepare pairwise widths and weights
    pσ, pϕ = pairwise_consts(gmmx, gmmy)

    # initialization
    if isnothing(initblock)
        initblock = blockfun(gmmx, gmmy, nothing, pσ, pϕ, rot, trl)
    end
    ndims = dims(initblock)
    ub, bestloc = initblock.upperbound, initblock.center    # best-so-far objective value and transformation
    t = promote_type(numbertype(gmmx), numbertype(gmmy))
    lb = typemin(t)
    pq = PriorityQueue{Block{ndims, t}, Tuple{t,t}}()
    enqueue!(pq, initblock, (initblock.lowerbound, initblock.upperbound))

    # split cubes until convergence
    ndivisions = 0
    sinceimprove = 0
    evalsperdiv = length(gmmx)*length(gmmy)*nsplits^ndims
    while !isempty(pq)
        if (length(pq) > maxblocks) || (ndivisions*evalsperdiv > maxevals) || (sinceimprove > maxstagnant) || (ndivisions > maxsplits)
            break
        end
        ndivisions += 1
        sinceimprove += 1

        # take the block with the lowest lower bound
        bl, (lb, blub) = dequeue_pair!(pq)

        # if the best solution so far is close enough to the best possible solution, end
        if abs((ub - lb)/lb) < rtol || abs(ub-lb) < atol
            return GMAlignmentResult(gmmx, gmmy, ub, lb, tformfun(bestloc...), bestloc, ndivisions*evalsperdiv, ndivisions, length(pq), sinceimprove)
        end

        # split up the block into `nsplits` smaller blocks across each dimension
        subrngs = subranges(bl.ranges, nsplits)
        sblks = fill(Block{ndims,t}(), nsplits^ndims)
        if threads
            Threads.@threads for i=1:length(subrngs)
                sblks[i] = blockfun(gmmx, gmmy, subrngs[i], pσ, pϕ, rot, trl)
            end
        else
            for i=1:length(subrngs)
                sblks[i] = blockfun(gmmx, gmmy, subrngs[i], pσ, pϕ, rot, trl)
            end
        end

        # reset the upper bound if appropriate
        minub, ubidx = findmin([sblk.upperbound for sblk in sblks])
        if minub < ub
            ub, bestloc = local_align(gmmx, gmmy, sblks[ubidx], pσ, pϕ; objfun=objfun, rot=rot, trl=trl)
            sinceimprove = 0
        end

        # only add sub-blocks to the queue if they present possibility for improvement
        for sblk in sblks
            if sblk.lowerbound < ub
                enqueue!(pq, sblk, (sblk.lowerbound, sblk.upperbound))
            end
        end
    end
    if isempty(pq)
        return GMAlignmentResult(gmmx, gmmy, ub, lb, tformfun(bestloc...), bestloc, ndivisions*evalsperdiv, ndivisions, length(pq), sinceimprove)
    else
        return GMAlignmentResult(gmmx, gmmy, ub, dequeue_pair!(pq)[2][1], tformfun(bestloc...), bestloc, ndivisions*evalsperdiv, ndivisions, length(pq), sinceimprove)
    end
end

"""
    result = rot_gogma_align(gmmx, gmmy; kwargs...)

Finds the globally optimal rotation for alignment between two isotropic Gaussian mixtures, `gmmx`
and `gmmy`, using the [GOGMA algorithm](https://arxiv.org/abs/1603.00150), for a given rotation, `rot`.
That is, only rigid translation is allowed.

For details about keyword arguments, see `gogma_align()`.
"""
function rot_gogma_align(gmmx::AbstractGMM, gmmy::AbstractGMM; kwargs...)
    return gogma_align(gmmx, gmmy; blockfun=rotBlock, objfun=rot_alignment_objective, tformfun=LinearMap, kwargs...)
end

"""
    result = trl_gogma_align(gmmx, gmmy; kwargs...)

Finds the globally optimal translation for alignment between two isotropic Gaussian mixtures, `gmmx`
and `gmmy`, using the [GOGMA algorithm](https://arxiv.org/abs/1603.00150), for a given translation, `trl`.
That is, only rigid rotation is allowed.

For details about keyword arguments, see `gogma_align()`.
"""
function trl_gogma_align(gmmx::AbstractGMM, gmmy::AbstractGMM;  kwargs...)
    return gogma_align(gmmx, gmmy; blockfun=trlBlock, objfun=trl_alignment_objective, tformfun=Translation, kwargs...)
end

ext: jl


text:
import DandelionWebSockets: backoff_min, backoff_max, reset

type FakeBackoff <: AbstractBackoff
    values::Vector{Float64}
    min::Float64
    max::Float64
    i::Int

    FakeBackoff(values::Vector{Float64}, min::Float64, max::Float64) =
        new(values, min, max, 1)
end

function (b::FakeBackoff)()
    v = b.values[b.i]
    b.i += 1
    v
end

backoff_min(b::FakeBackoff) = b.min
backoff_max(b::FakeBackoff) = b.max
reset(b::FakeBackoff) = b.i = 1

facts("Reconnect") do
    context("Backoff") do
        # Requirement
        # @7_2_3-1 Client backoff
        # @7_2_3-2 First client reconnect
        # @7_2_3-4 Subsequent reconnection attempts

        min = 30.0
        max = 200.0
        backoff_range = max - min
        backoff = Backoff(min, max)

        values = [backoff() for i in range(1, 12)]

        # The first three attempts won't reach half the backoff range.
        for i in range(1, 3)
            @fact values[i] < 0.5 * backoff_range --> true
        end

        # After twelve attempts the backoff is over 90% of the backoff range.
        @fact values[12] > 0.9 * backoff_range --> true
    end

    context("Backoff reset") do
        backoff = Backoff(0.0, 100.0)
        values = [backoff() for i in range(1, 12)]
        reset(backoff)
        @fact backoff() --> 0.0
    end

    context("Backoff limits") do
        min = 30.0
        max = 200.0
        backoff = Backoff(min, max)
        values = [backoff() for i in range(1, 100)]
        for v in values
            @fact v >= min --> true
            @fact v <= max --> true
        end

        @fact backoff_min(backoff) --> min
        @fact backoff_max(backoff) --> max
    end

    context("RandomizedBackoff") do
        rmin = 5.0
        rmax = 13.0
        r = 5.0
        backoff = FakeBackoff([5.0, 7.0, 9.0, 12.0, 13.0], rmin, rmax)
        expected_rng_values = [1.0, -3.0, 3.0, 1.0, 5.0]
        rng = FakeRNG{Float64}([x/(2*r) + 0.5 for x in expected_rng_values])
        expected = [6.0, 5.0, 12.0, 13.0, 13.0]

        rand_backoff = RandomizedBackoff(backoff, rng, r)
        values = [rand_backoff() for i in range(1, 5)]
        for i in range(1, 5)
            @fact values[i] --> roughly(expected[i])
        end
    end

    context("Retry") do
        expected_delays = [1.0, 2.0, 3.0, 4.0, 5.0]
        backoff = FakeBackoff(expected_delays, 0.0, 5.0)
        actual_delays = Vector{Float64}()

        retry_fun = () -> nothing
        timer_fun = (actual_retry_fun, delay) -> begin
            push!(actual_delays, delay)
            nothing
        end

        r = Retry(backoff, retry_fun; timer_fun=timer_fun)

        retry(r)
        retry(r)
        retry(r)
        retry(r)
        retry(r)

        @fact actual_delays --> expected_delays
    end

    context("Retry reset") do
        backoff_values = [1.0, 2.0, 3.0, 4.0, 5.0]
        backoff = FakeBackoff(backoff_values, 0.0, 5.0)
        expected_delays = [1.0, 2.0, 1.0, 2.0]
        actual_delays = Vector{Float64}()

        retry_fun = () -> nothing
        timer_fun = (actual_retry_fun, delay) -> begin
            push!(actual_delays, delay)
            nothing
        end

        r = Retry(backoff, retry_fun; timer_fun=timer_fun)

        retry(r)
        retry(r)
        reset(r)
        retry(r)
        retry(r)

        @fact actual_delays --> expected_delays
    end

    context("Set reconnect function") do
        backoff_values = [1.0, 2.0, 3.0, 4.0, 5.0]
        backoff = FakeBackoff(backoff_values, 0.0, 5.0)
        retry_fun = () -> nothing
        r = Retry(backoff, retry_fun)

        new_retry_fun = () -> 1+1
        set_function(r, new_retry_fun)
    end

    context("Default timer") do
        backoff_value = 0.2
        backoff = FakeBackoff([backoff_value], 0.0, 5.0)
        called = 0
        retry_fun = () -> called += 1
        r = Retry(backoff, retry_fun)

        retry(r)
        sleep(backoff_value * 2.0)
        @fact called --> 1
    end
end


text:
mutable struct LeastSquares{T} <: AbstractLoss
    A::Matrix{T}
    b::Vector{T}
    temp::Vector{T}

    function (::Type{LeastSquares})(A::Matrix, b::Vector)
        T = promote_type(eltype(Q), eltype(b))
        new{T}(A, b)
    end

    function (::Type{LeastSquares})(N::Integer, d::Integer)
        @assert d ≥ 1 "d must be at least 1"
        @assert N ≥ d "N must be at least $(d)"

        A = randn(N, d)
        A *= 2/norm(A)
        b = Float64[rand() ≤ 0.5 ? -1. : 1. for idx in 1:N]
        T = promote_type(eltype(A), eltype(b))
        new{T}(A, b, copy(b))
    end
end

nfeatures(loss::LeastSquares) = size(loss.A,2)
nsamples(loss::LeastSquares) = size(loss.A,1)

function loss!(loss::LeastSquares{T}, x::Vector{T}, g::Vector{T}) where T
    mul!(loss.temp,loss.A,x)
    result = convert(T, 0.5) * norm(loss.temp-loss.b,2)^2
    mul!(g, adjoint(loss.A), loss.temp)
    return result
end


text:
module binaryc

include("./file.jl")

@doc raw"""
Take the command line arguments and return 0 if there is none
else, check if the first argument endswith a binaryc file extension

if so, create a BinaryCFile object and compiler the compiled_code
else, throw an exception
"""
function parse_command_arguments(argument::Array{String, 1})::Int64
	if(length(argument) == 0)
		return 0
	end
	required_command = argument[1]
	if(endswith(required_command, ".binc"))
		file = BinaryCFile(required_command)
		is_valid_file(file)
	else
		error("$required_command not a binary c file")
	end

	return 1
end

parse_command_arguments(ARGS)

end

max_stars_repo_name: Technetium-Inc/BinaryC

max_issues_repo_name: Technetium-Inc/BinaryC

max_forks_repo_name: Technetium-Inc/BinaryC


text:
using ProtoBuf

run(ProtoBuf.protoc(`-I=/Users/peter/onnx/onnx --julia_out=/Users/peter/.julia/dev/OnnxIO/src /Users/peter/onnx/onnx/onnx.proto`))


text:
using InteractiveUtils

struct DoubleNumWiltonSauterQStrat{R,S}
    outer_rule_far::R
    inner_rule_far::R
    outer_rule_near::R
    inner_rule_near::R
    sauter_schwab_common_tetr::S
    sauter_schwab_common_face::S
    sauter_schwab_common_edge::S
    sauter_schwab_common_vert::S
end

struct DoubleNumQStrat{R}
    outer_rule::R
    inner_rule::R
end

# defaultquadstrat(op, tfs, bfs) = DoubleNumWiltonSauterQStrat(2,3,6,7,6,5,4,3)
defaultquadstrat(op, tfs, bfs) = defaultquadstrat(op, refspace(tfs), refspace(bfs))
# defaultquadstrat(op, ::RefSpace, ::RefSpace) = error("No default quadstrat set.")

struct SingleNumQStrat
    quad_rule::Int
end

function quadinfo(op, tfs, bfs; quadstrat=defaultquadstrat(op, tfs, bfs))

    tels, tad = assemblydata(tfs)
    bels, bad = assemblydata(bfs)

    tref = refspace(tfs)
    bref = refspace(bfs)

    i, τ = 1, first(tels)
    j, σ = 1, first(bels)

    @show quadstrat
    println(@which BEAST.quaddata(op,tref,bref,tels,bels,quadstrat))

    qd = quaddata(op,tref,bref,tels,bels,quadstrat)
    println(@which quadrule(op,tref,bref,i,τ,j,σ,qd,quadstrat))

    nothing
end


text:
# reciprocal.jl: data structures needed to store reciprocal space data
# TODO: I'm sure there's a Julia package that will generate paths in k-space for band structures
# perhaps we should try to integrate that functionality?

"""
    AbstractKPoints

Supertype for data structures that store k-point information associated with a structure.
"""
abstract type AbstractKPoints
end

"""
    KPointGrid{N} <: AbstractKPoints

Stores information about a k-point grid.

The storage of a k-point grid matches the convention used in abinit's `kptrlatt` variable. With
this convention, the k-point grid can be interpreted as defining a supercell in real space. Each
vector comprising the grid is given in terms of the real space primitive basis vectors. The shift
of the k-point grid off the reciprocal space origin (Γ) is also given.
"""
struct KPointGrid{N} <: AbstractKPoints
    # k-point grid
    grid::SVector{SVector{N,Int}}
    # shift of the k-point grid off Γ
    shift::SVector{N, Float64}

    function KPointGrid{N}(
        grid::AbstractVector{<:AbstractVector{<:Real}},
        shift::AbstractVector{<:Real}
    )
        # Verify that the grid vectors are linearly independent
        # TODO: this might need to be done independently for exterior constructors
        @assert _is_linearly_independent(grid...) "k-point grid vectors are not linearly independent"
        # Bring the shift within the first Brillouin zone
        shift = shift - floor.(shift)
        return new(grid, shift)
    end
end

"""
    KPointList{N}

Stores information about individual k-points used in a calculation. This may be useful for
calculations of band structures.
"""
struct KPointList{N} <: AbstractKPoints
    list::Vector{SVector{N,Float64}}
end

"""
    nkpt(k::KPointList{N}) -> Int

Returns the number of k-points in a `KPointList`.
"""
nkpt(k::KPointList{N}) = size(k.list)


text:
import SparseArrays: SparseVector, sparsevec
import LinearAlgebra: cholesky, Matrix

struct BasisFunction{PT,RT}
    center::PT
    coefficients::SparseVector{RT, Int}
end

function BasisFunction(center, coefficients)
    return BasisFunction{typeof(center), eltype(coefficients)}(center, coefficients)
end

function real_type(::Type{BasisFunction{PT, RT}})  where {PT, RT}
    return RT
end

function center(bf::BasisFunction)
    return bf.center
end

function coefficients(bf::BasisFunction)
    return bf.coefficients
end

function create_basis_vector(id_arrays, weights, N)
    @assert size(id_arrays) == size(weights)
    out_weights = vcat([weights[k] * one.(id_arrays[k]) for k = 1 : length(weights)]...)
    return sparsevec(vcat(id_arrays...), out_weights, N)
end

# Input are the coarsest-scale domains
function compute_basis_functions(domains::AbstractVector{Domain{PT}}, N = maximum(id.(gather_descendants(domains)))) where PT<:AbstractVector{<:Real}
    centers = Vector{Vector{PT}}(undef,1)
    centers[1] = Vector{PT}(undef, 0)
    basis_vectors = Vector{Vector{SparseVector{eltype(PT), Int}}}(undef,1)
    basis_vectors[1] = Vector{SparseVector{eltype(PT), Int}}(undef, 0)
    # recursive functions that computes the basis functions orthogonal in ℓ^2 inner product. Here, k is the level of domain
    function recurse(dm, k)
        # If the domain is elementrary, we return an array containint its index.
        n_ch = length(children(dm))
        if iselementary(dm)
            return [id(dm)]
        else
            # Extend length of output arrays if first domain on new level
            if k + 1 > length(basis_vectors)
                # make sure we didn't skip a level
                @assert length(basis_vectors) == length(centers) == k
                push!(centers, Vector{PT}(undef, 0))
                push!(basis_vectors, Vector{SparseVector{eltype(PT), Int}}(undef, 0))
            end

            # Compute coefficients of linear combination
            M = zeros(eltype(PT), n_ch + 1, n_ch + 1)
            M[1, 1] = sum(weight.(children(dm)))
            # setting the inner products.
            for i = 1 : n_ch;
                M[1, i + 1] = M[i + 1, 1] = M[i + 1, i + 1] = weight(children(dm)[i])
            end
            # Adding 1 to the last entry of M to ensure positivity#
            # This should not affect the result, since the affected columns are dropped
            M[end, end] += 1
            # Compute QR factorization that will provide coefficients
            R = inv(Matrix(cholesky(M).L'))

            # apply the function to the children to obtain the children arrays
            id_arrays = map(ch -> recurse(ch, k + 1), children(dm))
            id_arrays = vcat([vcat(id_arrays...)], id_arrays)

            # If on first level, return normalized average as basis function

            if k == 1
                push!(basis_vectors[k], create_basis_vector(id_arrays, R[:, 1], N))
                push!(centers[k], center(dm))
            end

            # add the remaining basis functions
            for ch = 2 : n_ch
                push!(basis_vectors[k + 1], create_basis_vector(id_arrays, R[:, ch], N))
                push!(centers[k + 1], center(dm))
            end
            return id_arrays[1]
        end
    end
    for dm in domains
        recurse(dm, 1)
    end

    # combine basis_vectors and centers to basis functions
    # cannot broadcast directly, because we are dealing with arrays of arrays of basis functions due to the different scales.
    out = Vector{Vector{BasisFunction{eltype(centers[1]), eltype(basis_vectors[1][1])}}}(undef, length(basis_vectors))
    for k = 1 : length(basis_vectors)
        out[k] = BasisFunction.(centers[k], basis_vectors[k])
    end

    return out
end


text:
include("scripts.jl")

function scr_ml(ARGS)

    # PARSING ARGUMENTS
    id_scr = parse(Int32, ARGS[1])
    # number of biotracers
    n_bio = parse(Int64, ARGS[2])
    # number of indiv used for building distribution
    n_distr = parse(Int64, ARGS[3]) # NB: not used here
    # number of sample used to test
    n_sample = parse(Int64, ARGS[4])
    # number of replicates
    n_rep = parse(Int64, ARGS[5])
    # maximum column combintaions
    mxcb = parse(Int64, ARGS[6])
    # path to data (pca-transformed or not)
    id_rep = parse(Int64, ARGS[7])
    # path to data
    inp="../data/data_f_cs.csv"
    # FIXED PARAMETERS
    # path to res folder
    # number of regions considered
    n_reg = 3
    n_tot_indiv = 30
    n_tot_bio = 17
    # data augmentation parameters
    aug = Int64(1e3)
    sig = .01

    df = get_data(false, inp)

    out = "./"

    println("Data loaded ✓")

    if id_scr == 1
        # NB n_distr is ignored
        scr_ndistr(df, n_bio, n_sample, mxcb, n_rep, n_tot_bio, n_tot_indiv,
            n_reg, aug, sig, id_rep)
    elseif id_scr == 2
        scr_nsample(df, n_bio, n_distr, n_sample, mxcb, n_rep, n_tot_bio,
        n_tot_indiv, n_reg, aug, sig, out, pca)
    elseif id_scr == 3
        scr_noise(df, n_bio, n_distr, n_sample, mxcb, n_rep, n_tot_bio,
        n_tot_indiv, n_reg, aug, sig, id_rep)
    else
        # NB n_bio used as n_bio_max
        scr_nbio_v2(df, n_bio, n_distr, n_sample, mxcb, n_rep, n_tot_bio,
        n_tot_indiv, n_reg, aug, sig, out, pca)
    end

    return 0
end

println("All scripts loaded ✓")
scr_ml(ARGS)


text:
module StochasticVI

end


text:
using Documenter
using unfold
using DocStringExtensions
makedocs(sitename="Unfold.jl",
        prettyurls = get(ENV, "CI", nothing) == "true",
        pages = [
            "index.md",
            "LM Tutorial" =>"lm_tutorial.md",
            "LMM Tutorial" =>"lmm_tutorial.md",
        ])


text:
using ProductizeBayes

bernoulli_model = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
";

data = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

sm = VersionedModel("bernoulli", bernoulli_model)
rc = stan_sample(sm.stan_model; data, num_chains=2, seed=-1)


text:
using LevelSetEvaluation

"""
Description:
    This example creates an automotive level set estimator, and uses it to
    evaluate the impact of two input variables (noise in the longitudinal and
    lateral acceleration) on an output variable (number of hard brakes).

    It then evaluates the impact of sensor noise on a uav collision avoidance
    simulator. In this case, the sensor noise is over the relative angle and
    heading of the other uav, and the target value is the probability of a
    near mid-air collision.
"""

function grid_evaluation!(lse::LevelSetEvaluator, params_1::Vector{Float64},
        params_2::  Vector{Float64})
    num_runs = length(params_1)
    targets = zeros(num_runs, num_runs)
    for (i, x1) in enumerate(params_1)
        for (j, x2) in enumerate(params_2)
            tic()
            targets[i, j] = evaluate_parameters!(lse, [x1,x2], 1, 1)
            print(@sprintf("input: (%03f, %03f)\toutput: %05f\t", x1, x2, targets[i,j]))
            toc()
        end
    end
    return targets
end

# auto
auto_lse = build_evaluator("auto")
num_runs = 2
param_1_values = collect(linspace(0., 2., num_runs))
param_2_values = collect(linspace(0, .5, num_runs))
targets = grid_evaluation!(auto_lse, param_1_values, param_2_values)
println(reshape(targets, (num_runs, num_runs)))

# uav
uav_lse = build_evaluator("uav")
num_runs = 5
param_1_values = collect(linspace(0, pi / 2, num_runs))
param_2_values = collect(linspace(0, pi / 2, num_runs))
targets = grid_evaluation!(uav_lse, param_1_values, param_2_values)
println(reshape(targets, (num_runs, num_runs)))


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -12.75118444054646 <= q <= 9.039664285972906)

                     add_NL_constraint(m, :(softplus(0.715796138346545 + -0.1516695982843861*softplus(-0.317423658088571 + -0.8843064794993589*softplus(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3])) + 0.6219431368805965*softplus(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3])) + -0.24757066957523488*softplus(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3]))) + -0.0820632819191216*softplus(-0.8960501996212993 + -0.7682728024607406*softplus(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3])) + -0.3378399703859083*softplus(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3])) + 0.4153867033969405*softplus(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3]))) + -0.0571486344530765*softplus(0.054021626349418206 + -0.05819989275612114*softplus(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3])) + 0.8875139231924525*softplus(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3])) + 0.11280184637312685*softplus(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3])))) + softplus(-0.8464220616292248 + -0.9331432460947604*softplus(-0.317423658088571 + -0.8843064794993589*softplus(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3])) + 0.6219431368805965*softplus(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3])) + -0.24757066957523488*softplus(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3]))) + 0.425513343659472*softplus(-0.8960501996212993 + -0.7682728024607406*softplus(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3])) + -0.3378399703859083*softplus(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3])) + 0.4153867033969405*softplus(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3]))) + 0.1421259709107252*softplus(0.054021626349418206 + -0.05819989275612114*softplus(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3])) + 0.8875139231924525*softplus(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3])) + 0.11280184637312685*softplus(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3])))) + softplus(-0.5205460270098605 + 0.13812073353255405*softplus(-0.317423658088571 + -0.8843064794993589*softplus(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3])) + 0.6219431368805965*softplus(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3])) + -0.24757066957523488*softplus(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3]))) + 0.7749706249828683*softplus(-0.8960501996212993 + -0.7682728024607406*softplus(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3])) + -0.3378399703859083*softplus(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3])) + 0.4153867033969405*softplus(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3]))) + 0.5007112790542019*softplus(0.054021626349418206 + -0.05819989275612114*softplus(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3])) + 0.8875139231924525*softplus(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3])) + 0.11280184637312685*softplus(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# NONLEPUS CONTROL
#integrator = DiscretizeRelax(prob, DynamicBoundspODEsDiscrete.LohnerContractor{4}(), h = 0.01, skip_step2 = false, relax = use_relax)
#integrator = DiscretizeRelax(prob, DynamicBoundspODEsDiscrete.LohnerContractor{4}(), h = 0.025, skip_step2 = false, relax = use_relax)

# LEPUS CONTROL
#integrator = DiscretizeRelax(prob, DynamicBoundspODEsDiscrete.LohnerContractor{6}(), h = 0.02,
#                             repeat_limit = 1, skip_step2 = false, step_limit = 5, relax = use_relax)

# TODO: need to adjust coefficient from hoe test...

#using Revise
using IntervalArithmetic, TaylorSeries
setrounding(Interval, :none)
import Base: literal_pow, ^
import IntervalArithmetic.pow
function ^(x::Interval{Float64}, n::Integer)  # fast integer power
    if n < 0
        return 1/IntervalArithmetic.pow(x, -n)
    end
    isempty(x) && return x
    if iseven(n) && 0 ∈ x
        return IntervalArithmetic.hull(zero(x),
                    hull(Base.power_by_squaring(Interval(mig(x)), n),
                        Base.power_by_squaring(Interval(mag(x)), n))
            )
    else
      return IntervalArithmetic.hull( Base.power_by_squaring(Interval(x.lo), n),
                    Base.power_by_squaring(Interval(x.hi), n) )
    end
end

using DynamicBoundsBase, Plots, DifferentialEquations#, Cthulhu
using DynamicBoundspODEsDiscrete, BenchmarkTools

println(" ")
println(" ------------------------------------------------------------ ")
println(" ------------- PACKAGE EXAMPLE       ------------------------ ")
println(" ------------------------------------------------------------ ")

use_relax = false
lohners_type = 2
prob_num = 1
ticks = 50.0
steps = 50.0
tend = 0.02*steps/ticks # lo 7.6100

if prob_num == 1
    x0(p) = [9.0]
    function f!(dx, x, p, t)
        dx[1] = p[1] - x[1]^2 #x[1]*x[1]
        nothing
    end
    tspan = (0.0, tend)
    pL = [-1.0]
    pU = [1.0]

elseif prob_num == 2
    x0(p) = [1.2; 1.1]
    function f!(dx, x, p, t)
        dx[1] = p[1]*x[1]*(one(typeof(p[1])) - x[2])
        dx[2] = p[1]*x[2]*(x[1] - one(typeof(p[1])))
        nothing
    end
    tspan = (0.0, tend)
    pL = [2.95]
    pU = [3.05]
end

prob = DynamicBoundsBase.ODERelaxProb(f!, tspan, x0, pL, pU)
tol = 1E-5

if lohners_type == 1
    integrator = DiscretizeRelax(prob, DynamicBoundspODEsDiscrete.LohnerContractor{5}(), h = 1/ticks,
                                 repeat_limit = 1, skip_step2 = false, step_limit = steps, relax = use_relax, tol= tol)
elseif lohners_type == 2
    integrator = DiscretizeRelax(prob, DynamicBoundspODEsDiscrete.HermiteObreschkoff(3, 3), h = 1/ticks,
                             repeat_limit = 1, skip_step2 = false, step_limit = steps, relax = use_relax, tol= tol)
elseif lohners_type == 3
    function iJx!(dx, x, p, t)
        dx[1] = -2.0*x[1]
        nothing
    end
    function iJp!(dx, x, p, t)
        dx[1] = one(p[1])
        nothing
    end
    integrator = DiscretizeRelax(prob, DynamicBoundspODEsDiscrete.AdamsMoulton(2), h = 1/ticks,
                                 repeat_limit = 1, step_limit = steps, skip_step2 = false,
                                 relax = false, Jx! = iJx!, Jp! = iJp!, tol= tol)
end

ratio = rand(1)
pstar = pL.*ratio .+ pU.*(1.0 .- ratio)
setall!(integrator, ParameterValue(), [0.0])
DynamicBoundsBase.relax!(integrator)
integrate!(integrator)

t_vec = integrator.time
if !use_relax
    lo_vec = getfield.(getindex.(integrator.storage[:],1), :lo)
    hi_vec = getfield.(getindex.(integrator.storage[:],1), :hi)
else
    lo_vec = getfield.(getfield.(getindex.(integrator.storage[:],1), :Intv), :lo)
    hi_vec = getfield.(getfield.(getindex.(integrator.storage[:],1), :Intv), :hi)
end

plt = plot(t_vec , lo_vec, label="Interval Bounds 0.0", marker = (:hexagon, 2, 0.6, :green), linealpha = 0.0, legend=:bottomleft)
plot!(plt, t_vec , hi_vec, label="", linealpha = 0.0, marker = (:hexagon, 2, 0.6, :green))

prob = ODEProblem(f!, [9.0], tspan, [-1.0])
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)
plot!(plt, sol.t , sol[1,:], label="", linecolor = :red, linestyle = :solid, lw=1.5)

prob = ODEProblem(f!, [9.0], tspan,[1.0])
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)
plot!(plt, sol.t , sol[1,:], label="", linecolor = :red, linestyle = :solid, lw=1.5)

ylabel!("x[1] (M)")
xlabel!("Time (seconds)")
display(plt)


text:
module Pljuto
using HypertextLiteral

function cljs_bootstrap()
    println("Loading CLJS dependencies...")
    @htl("""
    <script src="https://cdn.jsdelivr.net/gh/borkdude/scittle@0.0.2/js/scittle.js" type="application/javascript"></script>
    <script crossorigin src="https://unpkg.com/react@17/umd/react.production.min.js"></script>
    <script crossorigin src="https://unpkg.com/react-dom@17/umd/react-dom.production.min.js"></script>
    <script src="https://cdn.jsdelivr.net/gh/borkdude/scittle@0.0.2/js/scittle.reagent.js" type="application/javascript"></script>
    <script> scittle.core.disable_auto_eval(true);</script>""")
end

macro cljs_str(cljs::String)
    @htl("""<script>return scittle.core.eval_string($(cljs));</script>""")
end

export @cljs,@cljs_str,cljs_bootstrap
end


text:
const DEFAULT_VOIGT_ORDER = ([1], [1 3; 4 2], [1 6 5; 9 2 4; 8 7 3])
"""
    tovoigt(A::Union{SecondOrderTensor, FourthOrderTensor}; kwargs...)

Converts a tensor to "Voigt"-format.

Keyword arguments:
 - `offdiagscale`: determines the scaling factor for the offdiagonal elements.
   `frommandel` can also be used for "Mandel"-format which sets `offdiagscale = √2`.
   This argument is only applicable for `SymmetricTensor`s.
 - `order`: matrix of the linear indices determining the Voigt order. The default
   index order is `[11, 22, 33, 23, 13, 12, 32, 31, 21]`, corresponding to
   `order = [1 6 5; 9 2 4; 8 7 3]`.

See also [`tovoigt!`](@ref) and [`fromvoigt`](@ref).

```jldoctest
julia> tovoigt(Tensor{2,3}(1:9))
9-element Array{Int64,1}:
 1
 5
 9
 8
 7
 4
 6
 3
 2

julia> tovoigt(SymmetricTensor{2,3}(1:6); offdiagscale = 2)
6-element Array{Int64,1}:
  1
  4
  6
 10
  6
  4

julia> tovoigt(Tensor{4,2}(1:16))
4×4 Array{Int64,2}:
 1  13   9  5
 4  16  12  8
 3  15  11  7
 2  14  10  6
```
"""
function tovoigt end
@inline function tovoigt(A::Tensor{2, dim, T, M}; order=DEFAULT_VOIGT_ORDER[dim]) where {dim, T, M}
    @inbounds tovoigt!(Vector{T}(undef, M), A; order=order)
end
@inline function tovoigt(A::Tensor{4, dim, T, M}; order=DEFAULT_VOIGT_ORDER[dim]) where {dim, T, M}
    @inbounds tovoigt!(Matrix{T}(undef, Int(√M), Int(√M)), A; order=order)
end
@inline function tovoigt(A::SymmetricTensor{2, dim, T, M}; offdiagscale::T=one(T), order=DEFAULT_VOIGT_ORDER[dim]) where {dim, T, M}
    @inbounds tovoigt!(Vector{T}(undef, M), A; offdiagscale=offdiagscale, order=order)
end
@inline function tovoigt(A::SymmetricTensor{4, dim, T, M}; offdiagscale::T=one(T), order=DEFAULT_VOIGT_ORDER[dim]) where {dim, T, M}
    @inbounds tovoigt!(Matrix{T}(undef, Int(√M), Int(√M)), A; offdiagscale=offdiagscale, order=order)
end

"""
    tovoigt!(v::Array, A::Union{SecondOrderTensor, FourthOrderTensor}; kwargs...)

Converts a tensor to "Voigt"-format using the following index order:
`[11, 22, 33, 23, 13, 12, 32, 31, 21]`.

Keyword arguments:
 - `offset`: offset index for where in the array `A` the tensor should be stored.
   For 4th order tensors the keyword arguments are `offset_i` and `offset_j`,
   respectively. Defaults to `0`.
 - `offdiagscale`: determines the scaling factor for the offdiagonal elements.
   `frommandel` can also be used for "Mandel"-format which sets `offdiagscale = √2`.
   This argument is only applicable for `SymmetricTensor`s.
 - `order`: matrix of the linear indices determining the Voigt order. The default
   index order is `[11, 22, 33, 23, 13, 12, 32, 31, 21]`.

See also [`tovoigt`](@ref) and [`fromvoigt`](@ref).

```jldoctest
julia> T = rand(Tensor{2,2})
2×2 Tensor{2,2,Float64,4}:
 0.794026  0.200586
 0.854147  0.298614

julia> x = zeros(4);

julia> tovoigt!(x, T)
4-element Array{Float64,1}:
 0.7940257103317943
 0.2986142783434118
 0.20058603493384108
 0.8541465903790502

julia> x = zeros(5);

julia> tovoigt!(x, T; offset=1)
5-element Array{Float64,1}:
 0.0
 0.7940257103317943
 0.2986142783434118
 0.20058603493384108
 0.8541465903790502
```
"""
function tovoigt! end
Base.@propagate_inbounds function tovoigt!(v::AbstractVector, A::Tensor{2, dim}; offset::Int=0, order=DEFAULT_VOIGT_ORDER[dim]) where {dim}
    for j in 1:dim, i in 1:dim
        v[offset + order[i, j]] = A[i, j]
    end
    return v
end
Base.@propagate_inbounds function tovoigt!(v::AbstractMatrix, A::Tensor{4, dim}; offset_i::Int=0, offset_j::Int=0, order=DEFAULT_VOIGT_ORDER[dim]) where {dim}
    for l in 1:dim, k in 1:dim, j in 1:dim, i in 1:dim
        v[offset_i + order[i, j], offset_j + order[k, l]] = A[i, j, k, l]
    end
    return v
end
Base.@propagate_inbounds function tovoigt!(v::AbstractVector{T}, A::SymmetricTensor{2, dim}; offdiagscale::T=one(T), offset::Int=0, order=DEFAULT_VOIGT_ORDER[dim]) where {T, dim}
    for j in 1:dim, i in 1:j
        v[offset + order[i, j]] = i == j ? A[i, j] : A[i, j] * offdiagscale
    end
    return v
end
Base.@propagate_inbounds function tovoigt!(v::AbstractMatrix{T}, A::SymmetricTensor{4, dim}; offdiagscale::T=one(T), offset_i::Int=0, offset_j::Int=0, order=DEFAULT_VOIGT_ORDER[dim]) where {T, dim}
    for l in 1:dim, k in 1:l, j in 1:dim, i in 1:j
        v[offset_i + order[i, j], offset_j + order[k, l]] =
            (i == j && k == l) ? A[i, j, k, l] :
            (i == j || k == l) ? A[i, j, k, l] * offdiagscale :
                                 A[i, j, k, l] * (offdiagscale * offdiagscale)
    end
    return v
end

@inline function tomandel(A::SymmetricTensor{o, dim, T}; order=DEFAULT_VOIGT_ORDER[dim]) where {o, dim, T}
    return @inbounds tovoigt(A; offdiagscale=T(√2), order=order)
end
Base.@propagate_inbounds function tomandel!(v::AbstractVector{T}, A::SymmetricTensor{2}; offset::Int=0, order=DEFAULT_VOIGT_ORDER[dim]) where {T}
    tovoigt!(v, A; offdiagscale=T(√2), offset=offset, order=order)
end
Base.@propagate_inbounds function tomandel!(v::AbstractMatrix{T}, A::SymmetricTensor{4}; offset_i::Int=0, offset_j::Int=0, order=DEFAULT_VOIGT_ORDER[dim]) where {T}
    tovoigt!(v, A; offdiagscale=T(√2), offset_i=offset_i, offset_j=offset_j, order=order)
end

"""
    fromvoigt(S::Type{<:AbstractTensor}, A::Array{T}; kwargs...)

Converts an array `A` stored in Voigt format to a Tensor of type `S`.

Keyword arguments:
 - `offset`: offset index for where in the array `A` the tensor starts. For 4th order
   tensors the keyword arguments are `offset_i` and `offset_j`, respectively.
   Defaults to `0`.
 - `offdiagscale`: determines the inverse scaling factor on the offdiagonal elements.
   `frommandel` can also be used for "Mandel"-format which sets `offdiagscale = √2`.
   This argument is only applicable for `SymmetricTensor`s.

See also [`tovoigt`](@ref).

```jldoctest
julia> fromvoigt(Tensor{2,3}, 1.0:1.0:9.0)
3×3 Tensor{2,3,Float64,9}:
 1.0  6.0  5.0
 9.0  2.0  4.0
 8.0  7.0  3.0
```
"""
Base.@propagate_inbounds function fromvoigt(TT::Type{<: Tensor{2, dim}}, v::AbstractVector; offset::Int=0, order=DEFAULT_VOIGT_ORDER[dim]) where {dim}
    return TT(function (i, j); return v[offset + order[i, j]]; end)
end
Base.@propagate_inbounds function fromvoigt(TT::Type{<: Tensor{4, dim}}, v::AbstractMatrix; offset_i::Int=0, offset_j::Int=0, order=DEFAULT_VOIGT_ORDER[dim]) where {dim}
    return TT(function (i, j, k, l); return v[offset_i + order[i, j], offset_j + order[k, l]]; end)
end
Base.@propagate_inbounds function fromvoigt(TT::Type{<: SymmetricTensor{2, dim}}, v::AbstractVector{T}; offdiagscale::T = T(1), offset::Int=0, order=DEFAULT_VOIGT_ORDER[dim]) where {dim, T}
    return TT(function (i, j)
            i > j && ((i, j) = (j, i))
            i == j ? (return v[offset + order[i, j]]) :
                     (return v[offset + order[i, j]] / offdiagscale)
        end)
end
Base.@propagate_inbounds function fromvoigt(TT::Type{<: SymmetricTensor{4, dim}}, v::AbstractMatrix{T}; offdiagscale::T = T(1), offset_i::Int=0, offset_j::Int=0, order=DEFAULT_VOIGT_ORDER[dim]) where {dim, T}
    return TT(function (i, j, k, l)
            i > j && ((i, j) = (j, i))
            k > l && ((k, l) = (l, k))
            i == j && k == l ? (return v[offset_i + order[i, j], offset_j + order[k, l]]) :
            i == j || k == l ? (return v[offset_i + order[i, j], offset_j + order[k, l]] / offdiagscale) :
                               (return v[offset_i + order[i, j], offset_j + order[k, l]] / (offdiagscale * offdiagscale))
        end)
end

Base.@propagate_inbounds function frommandel(TT::Type{<: SymmetricTensor{2,dim}}, v::AbstractVector{T}; offset::Int=0, order=DEFAULT_VOIGT_ORDER[dim]) where {T,dim}
    fromvoigt(TT, v, offdiagscale=T(√2), offset=offset, order=order)
end
Base.@propagate_inbounds function frommandel(TT::Type{<: SymmetricTensor{4,dim}}, v::AbstractMatrix{T}; offset_i::Int=0, offset_j::Int=0, order=DEFAULT_VOIGT_ORDER[dim]) where {T,dim}
    fromvoigt(TT, v, offdiagscale=T(√2), offset_i=offset_i, offset_j=offset_j, order=order)
end


text:
# Users of this module should first read the document
# https://github.com/alan-turing-institute/MLJ.jl/blob/master/doc/adding_new_models.md

module MLJBase

export MLJType, Model, Supervised, Unsupervised, Deterministic, Probabilistic
export Rows, Cols, Schema, retrieve, getrows
export fit, update, clean!, info, coerce
export predict, predict_mean, predict_mode
export transform, inverse_transform, se, evaluate, best
export target_kind, target_quantity, inputs_can_be, is_pure_julia
export package_url, package_name, package_uuid

export HANDLE_GIVEN_ID, @show, @constant  # from show.jl
export UnivariateNominal, average         # from distributions.jl

# methods from other packages to be rexported:
export pdf, mean, mode

import Base.==
using Query
import TableTraits
import DataFrames                # TODO: get rid of this dependency
import Distributions
import Distributions: pdf, mode
using CategoricalArrays

# from Standard Library:
using Statistics

## CONSTANTS

# the directory containing this file:
const srcdir = dirname(@__FILE__)
# horizontal space for field names in `MLJType` object display:
const COLUMN_WIDTH = 24
# how deep to display fields of `MLJType` objects:
const DEFAULT_SHOW_DEPTH = 1

## ABSTRACT TYPES

# overarching MLJ type:
abstract type MLJType end

# for storing hyperparameters:
abstract type Model <: MLJType end

abstract type Supervised{R} <: Model end # parameterized by fit-result type `R`
abstract type Unsupervised <: Model  end

# supervised models that `predict` probability distributions are of:
abstract type Probabilistic{R} <: Supervised{R} end

# supervsied models that `predict` point-values are of:
abstract type Deterministic{R} <: Supervised{R} end

# for displaying objects of `MLJType`:
include("show.jl")

# probability distributions and methods not provided by
# Distributions.jl package:
include("distributions.jl")

# convenience methods for manipulating categorical and tabular data
include("data.jl")

## THE MODEL INTERFACE

# every model interface must implement a `fit` method of the form
# `fit(model, verbosity, X, y) -> fitresult, cache, report` or
# `fit(model, verbosity, X, ys...) -> fitresult, cache, report` (multivariate case)
# or, one the simplified versions
# `fit(model, X, y) -> fitresult`
# `fit(model, X, ys...) -> fitresult`
fit(model::Model, verbosity::Int, args...) = fit(model, args...), nothing, nothing

# each model interface may optionally overload the following refitting
# method:
update(model::Model, verbosity, fitresult, cache, args...) =
    fit(model, verbosity, args...)

# methods dispatched on a model and fit-result are called *operations*.
# supervised models must implement this operation:
function predict end

# probabilistic supervised models may also overload one or more of
# these operations delevering point predictions:
predict_mode(model::Supervised, fitresult, Xnew) =
    mode.(predict(model, fitresult, Xnew))
predict_mean(model::Supervised, fitresult, Xnew) =
    mean.(predict(model, fitresult, Xnew))
predict_median(model::Supervised, fitresult, Xnew) =
    median.(predict(model, fitresult, Xnew))

# unsupervised methods must implement this operation:
function transform end

# unsupervised methods may implement this operation:
function inverse_transform end

# operations implemented by some meta-models:
function se end
function evaluate end
function best end

# a model wishing invalid hyperparameters to be corrected with a
# warning should overload this method (return value is the warning
# message):
clean!(model::Model) = ""

# supervised models may need to overload the following method to
# ensure the input data supplied by user (something implementing the
# Queryverse iterable table API) is coerced into the form required by
# its `fit` method and operations:
coerce(model::Model, Xtable) = Xtable

# if the return type `TABLE` of `coerce` is not a Queryverse
# iterable table or a `Matrix` (with rows
# corresponding to patterns and columns corresponding to features)
# then users will not be able to use MLJ's performant `EnsembleModel`
# on `model` unless one overloads the following method for type
# `TABLE`:
getrows(model::Model, X, r) = retrieve(X, Rows, r)
# here `r` is any integer, unitrange or colon `:`, and the right-hand
# side defined in `data.jl`.

# fallback trait declarations:
target_kind(::Type{<:Supervised}) = :unknown
target_quantity(::Type{<:Supervised}) = :univariate
inputs_can_be(::Type{<:Model}) = Symbol[]
is_pure_julia(::Type{<:Model}) = :unknown
package_name(::Type{<:Model}) = "unknown"
package_uuid(::Type{<:Model}) = "unknown"
package_url(::Type{<:Model}) = "unknown"

_response(::Type{<:Supervised}) = :unknown
_response(::Type{<:Deterministic}) = :deterministic
_response(::Type{<:Probabilistic}) = :probabilistic

target_is(modeltype::Type{<:Supervised}) =
    [_response(modeltype), target_kind(modeltype), target_quantity(modeltype)]

if VERSION < v"1.0.0"
    import Base.info
end

function info(modeltype::Type{<:Model})

    message = "$modeltype has a bad trait declaration."

    if modeltype <: Supervised
        target_kind(modeltype) in [:numeric, :binary, :multiclass, :unknown] ||
            error(message*"target_kind must return :numeric, :binary, :multiclass (or :unknown).")
        target_quantity(modeltype) in [:univariate, :multivariate] ||
            error(message*"target_quantity must return :univariate or :multivariate")
    end

    issubset(Set(inputs_can_be(modeltype)), Set([:numeric, :nominal, :missing])) ||
        error(message*"inputs_can_be must return a vector with entries from [:numeric, :nominal, :missing]")
    is_pure_julia(modeltype) in [:yes, :no, :unknown] ||
        error(message*"is_pure_julia must return :yes, :no (or :unknown).")

    # modelnamesplit = split(string(modeltype.name), '.')
    # if length(modelnamesplit) > 1
    #     modelnamesplit = modelnamesplit[2:end]
    # end
    # modelname = string(reduce((a, b)->"$a.$b", modelnamesplit))

    d = Dict{Symbol,Union{Symbol,Vector{Symbol},String}}()
    d[:model_name] = string(modeltype)
    d[:inputs_can_be] = inputs_can_be(modeltype)
    d[:is_pure_julia] = is_pure_julia(modeltype)
    d[:package_name] = package_name(modeltype)
    d[:package_uuid] = package_uuid(modeltype)
    d[:package_url] = package_url(modeltype)
    if modeltype <: Supervised
        d[:target_is] = target_is(modeltype)
    end

    return d
end
info(model::Supervised) = info(typeof(model))

# models are `==` if they have the same type and their field values are `==`:
function ==(m1::M, m2::M) where M<:Model
    ret = true
    for fld in fieldnames(M)
        ret = ret && getfield(m1, fld) == getfield(m2, fld)
    end
    return ret
end

end # module


text:
##############################################################################
#
# Text output
#
##############################################################################

function escapedprint(io::IO, x::Any, escapes::AbstractString)
    ourshowcompact(io, x)
end

function escapedprint(io::IO, x::AbstractString, escapes::AbstractString)
    escape_string(io, x, escapes)
end

function printtable(io::IO,
                    df::AbstractDataFrame;
                    header::Bool = true,
                    separator::Char = ',',
                    quotemark::Char = '"',
                    missingstring::AbstractString = "missing")
    n, p = size(df)
    etypes = eltypes(df)
    if header
        cnames = _names(df)
        for j in 1:p
            print(io, quotemark)
            print(io, cnames[j])
            print(io, quotemark)
            if j < p
                print(io, separator)
            else
                print(io, '\n')
            end
        end
    end
    quotestr = string(quotemark)
    for i in 1:n
        for j in 1:p
            if !ismissing(df[j][i])
                if ! (etypes[j] <: Real)
                    print(io, quotemark)
                    escapedprint(io, df[i, j], quotestr)
                    print(io, quotemark)
                else
                    print(io, df[i, j])
                end
            else
                print(io, missingstring)
            end
            if j < p
                print(io, separator)
            else
                print(io, '\n')
            end
        end
    end
    return
end

function printtable(df::AbstractDataFrame;
                    header::Bool = true,
                    separator::Char = ',',
                    quotemark::Char = '"',
                    missingstring::AbstractString = "missing")
    printtable(stdout,
               df,
               header = header,
               separator = separator,
               quotemark = quotemark,
               missingstring = missingstring)
    return
end
##############################################################################
#
# HTML output
#
##############################################################################

function html_escape(cell::AbstractString)
    cell = replace(cell, "&"=>"&")
    cell = replace(cell, "<"=>"<")
    cell = replace(cell, ">"=>">")
    return cell
end

function Base.show(io::IO, ::MIME"text/html", df::AbstractDataFrame)
    cnames = _names(df)
    write(io, "
")
    write(io, "
")
    write(io, "
")
    write(io, "	")
    for column_name in cnames
        write(io, "$(html_escape(String(column_name)))	")
    end
    write(io, "
")
    write(io, "
")
    write(io, "	")
    for j in 1:ncol(df)
        s = html_escape(compacttype(eltype(df[j])))
        write(io, "$s	")
    end
    write(io, "
")
    write(io, "
")
    write(io, "
")
    haslimit = get(io, :limit, true)
    n = size(df, 1)
    if haslimit
        tty_rows, tty_cols = displaysize(io)
        mxrow = min(n,tty_rows)
    else
        mxrow = n
    end
    for row in 1:mxrow
        write(io, "
")
        write(io, "$row	")
        for column_name in cnames
            if isassigned(df[column_name], row)
                cell = sprint(ourshowcompact, df[row, column_name])
            else
                cell = sprint(ourshowcompact, Base.undef_ref_str)
            end
            write(io, "$(html_escape(cell))	")
        end
        write(io, "
")
    end
    if n > mxrow
        write(io, "
")
        write(io, "⋮	")
        for column_name in cnames
            write(io, "⋮	")
        end
        write(io, "
")
    end
    write(io, "
")
    write(io, "
")
end

##############################################################################
#
# LaTeX output
#
##############################################################################

function latex_char_escape(char::Char)
    if char == '\\'
        return "\\textbackslash{}"
    elseif char == '~'
        return "\\textasciitilde{}"
    else
        return string('\\', char)
    end
end

function latex_escape(cell::AbstractString)
    replace(cell, ['\\','~','#','$','%','&','_','^','{','}']=>latex_char_escape)
end

function Base.show(io::IO, ::MIME"text/latex", df::AbstractDataFrame)
    nrows = size(df, 1)
    ncols = size(df, 2)
    cnames = _names(df)
    alignment = repeat("c", ncols)
    write(io, "\\begin{tabular}{r|")
    write(io, alignment)
    write(io, "}\n")
    write(io, "\t& ")
    header = join(map(c -> latex_escape(string(c)), cnames), " & ")
    write(io, header)
    write(io, "\\\\\n")
    write(io, "\t\\hline\n")
    for row in 1:nrows
        write(io, "\t")
        write(io, @sprintf("%d", row))
        for col in 1:ncols
            write(io, " & ")
            cell = isassigned(df[col], row) ? df[row,col] : Base.undef_ref_str
            if !ismissing(cell)
                if showable(MIME("text/latex"), cell)
                    show(io, MIME("text/latex"), cell)
                else
                    print(io, latex_escape(sprint(ourshowcompact, cell)))
                end
            end
        end
        write(io, " \\\\\n")
    end
    write(io, "\\end{tabular}\n")
end

##############################################################################
#
# MIME
#
##############################################################################

function Base.show(io::IO, ::MIME"text/csv", df::AbstractDataFrame)
    printtable(io, df, true, ',')
end

function Base.show(io::IO, ::MIME"text/tab-separated-values", df::AbstractDataFrame)
    printtable(io, df, true, '\t')
end

##############################################################################
#
# DataStreams-based IO
#
##############################################################################

using DataStreams, WeakRefStrings

struct DataFrameStream{T}
    columns::T
    header::Vector{String}
end
DataFrameStream(df::DataFrame) = DataFrameStream(Tuple(columns(df)), string.(names(df)))

# DataFrame Data.Source implementation
Data.schema(df::DataFrame) =
    Data.Schema(Type[eltype(A) for A in columns(df)], string.(names(df)), size(df, 1))

Data.isdone(source::DataFrame, row, col, rows, cols) = row > rows || col > cols
function Data.isdone(source::DataFrame, row, col)
    cols = length(source)
    return Data.isdone(source, row, col, cols == 0 ? 0 : length(source.columns[1]), cols)
end

Data.streamtype(::Type{DataFrame}, ::Type{Data.Column}) = true
Data.streamtype(::Type{DataFrame}, ::Type{Data.Field}) = true

Data.streamfrom(source::DataFrame, ::Type{Data.Column}, ::Type{T}, row, col) where {T} =
    source[col]
Data.streamfrom(source::DataFrame, ::Type{Data.Field}, ::Type{T}, row, col) where {T} =
    source[col][row]

# DataFrame Data.Sink implementation
Data.streamtypes(::Type{DataFrame}) = [Data.Column, Data.Field]
Data.weakrefstrings(::Type{DataFrame}) = true

allocate(::Type{T}, rows, ref) where {T} = Vector{T}(undef, rows)
allocate(::Type{CategoricalString{R}}, rows, ref) where {R} = CategoricalArray{String, 1, R}(undef, rows)
allocate(::Type{Union{CategoricalString{R}, Missing}}, rows, ref) where {R} =
    CategoricalArray{Union{String, Missing}, 1, R}(undef, rows)
allocate(::Type{CategoricalValue{T, R}}, rows, ref) where {T, R} =
    CategoricalArray{T, 1, R}(undef, rows)
allocate(::Type{Union{Missing, CategoricalValue{T, R}}}, rows, ref) where {T, R} =
    CategoricalArray{Union{Missing, T}, 1, R}(undef, rows)
allocate(::Type{WeakRefString{T}}, rows, ref) where {T} =
    WeakRefStringArray(ref, WeakRefString{T}, rows)
allocate(::Type{Union{Missing, WeakRefString{T}}}, rows, ref) where {T} =
    WeakRefStringArray(ref, Union{Missing, WeakRefString{T}}, rows)
allocate(::Type{Missing}, rows, ref) = missings(rows)

# Construct or modify a DataFrame to be ready to stream data from a source with `sch`
function DataFrame(sch::Data.Schema{R}, ::Type{S}=Data.Field,
                   append::Bool=false, args...;
                   reference::Vector{UInt8}=UInt8[]) where {R, S <: Data.StreamType}
    types = Data.types(sch)
    if !isempty(args) && args[1] isa DataFrame && types == Data.types(Data.schema(args[1]))
        # passing in an existing DataFrame Sink w/ same types as source
        sink = args[1]
        sinkrows = size(Data.schema(sink), 1)
        # are we appending and either column-streaming or there are an unknown # of rows
        if append && (S == Data.Column || !R)
            sch.rows = sinkrows
            # dont' need to do anything because:
              # for Data.Column, we just append columns anyway (see Data.streamto! below)
              # for Data.Field, unknown # of source rows, so we'll just push! in streamto!
        else
            # need to adjust the existing sink
            # similar to above, for Data.Column or unknown # of rows for Data.Field,
                # we'll append!/push! in streamto!, so we empty! the columns
            # if appending, we want to grow our columns to be able to include every row
                # in source (sinkrows + sch.rows)
            # if not appending, we're just "re-using" a sink, so we just resize it
                # to the # of rows in the source
            newsize = ifelse(S == Data.Column || !R, 0,
                        ifelse(append, sinkrows + sch.rows, sch.rows))
            foreach(col->resize!(col, newsize), columns(sink))
            sch.rows = newsize
        end
        # take care of a possible reference from source by addint to WeakRefStringArrays
        if !isempty(reference)
            foreach(col-> col isa WeakRefStringArray && push!(col.data, reference),
                sink.columns)
        end
        sink = DataFrameStream(sink)
    else
        # allocating a fresh DataFrame Sink; append is irrelevant
        # for Data.Column or unknown # of rows in Data.Field, we only ever append!,
            # so just allocate empty columns
        rows = ifelse(S == Data.Column, 0, ifelse(!R, 0, sch.rows))
        names = Data.header(sch)
        sink = DataFrameStream(
                Tuple(allocate(types[i], rows, reference) for i = 1:length(types)), names)
        sch.rows = rows
    end
    return sink
end

DataFrame(sink, sch::Data.Schema, ::Type{S}, append::Bool;
          reference::Vector{UInt8}=UInt8[]) where {S} =
    DataFrame(sch, S, append, sink; reference=reference)

@inline Data.streamto!(sink::DataFrameStream, ::Type{Data.Field}, val,
                      row, col::Int) =
    (A = sink.columns[col]; row > length(A) ? push!(A, val) : setindex!(A, val, row))
@inline Data.streamto!(sink::DataFrameStream, ::Type{Data.Field}, val,
                       row, col::Int, ::Type{Val{false}}) =
    push!(sink.columns[col], val)
@inline Data.streamto!(sink::DataFrameStream, ::Type{Data.Field}, val,
                       row, col::Int, ::Type{Val{true}}) =
    sink.columns[col][row] = val
@inline function Data.streamto!(sink::DataFrameStream, ::Type{Data.Column}, column,
                       row, col::Int, knownrows)
    append!(sink.columns[col], column)
end

Data.close!(df::DataFrameStream) =
    DataFrame(collect(AbstractVector, df.columns), Symbol.(df.header), makeunique=true)


text:
using QuantumPropagators
using Test
using SafeTestsets

# Note: comment outer @testset to stop after first @safetestset failure
@time @testset verbose = true "QuantumPropagators" begin

    print("\n* Spec.Rad. (test_specrad.jl):")
    @time @safetestset "Spec.Rad." begin
        include("test_specrad.jl")
    end

    print("\n* Cheby (test_cheby.jl):")
    @time @safetestset "Cheby" begin
        include("test_cheby.jl")
    end

    print("\n* Newton (test_newton.jl):")
    @time @safetestset "Newton" begin
        include("test_newton.jl")
    end

    print("\n* Exp (test_expprop.jl):")
    @time @safetestset "Exp" begin
        include("test_expprop.jl")
    end

    print("\n* Propagate (test_propagate.jl):")
    @time @safetestset "Propagate" begin
        include("test_propagate.jl")
    end

    print("\n")

end;


text:
module TestWithprogressMacro

using Logging
using ProgressLogging
using ProgressLogging: ProgressLevel, ROOTID
using Test
using Test: collect_test_logs

@testset "specify name in `@logprogress`" begin
    logs, = collect_test_logs(min_level = ProgressLevel) do
        @withprogress @logprogress "hello" 0.1
    end
    @test length(logs) == 3
    @test logs[1].kwargs[:progress] === NaN
    @test logs[2].kwargs[:progress] === 0.1
    @test logs[3].kwargs[:progress] === "done"
    @test logs[1].message == ""
    @test logs[2].message == "hello"
    @test logs[3].message == ""
    @test [l.message.progress.fraction for l in logs] == [nothing, 0.1, nothing]
    @test [l.message.progress.done for l in logs] == [false, false, true]
    @test length(unique([l.message.progress.id for l in logs])) == 1
    @test length(unique([l.id for l in logs])) == 1
    @test unique([l.message.progress.parentid for l in logs]) == [ROOTID]
end

@testset "specify name in `@withprogress`" begin
    logs, = collect_test_logs(min_level = ProgressLevel) do
        @withprogress name = "hello" @logprogress 0.1
    end
    @test length(logs) == 3
    @test logs[1].kwargs[:progress] === NaN
    @test logs[2].kwargs[:progress] === 0.1
    @test logs[3].kwargs[:progress] === "done"
    @test logs[1].message == "hello"
    @test logs[2].message == "hello"
    @test logs[3].message == "hello"
    @test [l.message.progress.fraction for l in logs] == [nothing, 0.1, nothing]
    @test [l.message.progress.done for l in logs] == [false, false, true]
    @test length(unique([l.message.progress.id for l in logs])) == 1
    @test length(unique([l.id for l in logs])) == 1
    @test unique([l.message.progress.parentid for l in logs]) == [ROOTID]
end

@testset "keyword argument when no name" begin
    logs, = collect_test_logs(min_level = ProgressLevel) do
        @withprogress @logprogress 0.1 message = "hello"
    end
    @test length(logs) == 3
    @test logs[1].kwargs[:progress] === NaN
    @test logs[2].kwargs[:progress] === 0.1
    @test logs[3].kwargs[:progress] === "done"
    @test logs[2].kwargs[:message] === "hello"
    @test [l.message.progress.fraction for l in logs] == [nothing, 0.1, nothing]
    @test [l.message.progress.done for l in logs] == [false, false, true]
    @test length(unique([l.message.progress.id for l in logs])) == 1
    @test length(unique([l.id for l in logs])) == 1
    @test unique([l.message.progress.parentid for l in logs]) == [ROOTID]
end

@testset "change name" begin
    logs, = collect_test_logs(min_level = ProgressLevel) do
        @withprogress name = "name" @logprogress "hello" 0.1
    end
    @test length(logs) == 3
    @test logs[1].kwargs[:progress] === NaN
    @test logs[2].kwargs[:progress] === 0.1
    @test logs[3].kwargs[:progress] === "done"
    @test logs[1].message == "name"
    @test logs[2].message == "hello"
    @test logs[3].message == "name"
    @test [l.message.progress.fraction for l in logs] == [nothing, 0.1, nothing]
    @test [l.message.progress.done for l in logs] == [false, false, true]
    @test length(unique([l.message.progress.id for l in logs])) == 1
    @test length(unique([l.id for l in logs])) == 1
    @test unique([l.message.progress.parentid for l in logs]) == [ROOTID]
end

@testset "nested" begin
    logs, = collect_test_logs(min_level = ProgressLevel) do
        @withprogress begin
            @logprogress "hello" 0.1
            @withprogress begin
                @logprogress "world" 0.2
            end
        end
    end

    @test length(logs) == 6

    ids = unique([l.id for l in logs])
    @test length(ids) == 2

    @test Tuple(
        (l.id, l.message.progress.parentid, string(l.message), l.kwargs[:progress])
        for l in logs
    ) === (
        (ids[1], ROOTID, "", NaN),
        (ids[1], ROOTID, "hello", 0.1),
        (ids[2], ids[1], "", NaN),
        (ids[2], ids[1], "world", 0.2),
        (ids[2], ids[1], "", "done"),
        (ids[1], ROOTID, "", "done"),
    )
end

@testset "`@logprogress 1.0` is not `done`" begin
    logs, = collect_test_logs(min_level = ProgressLevel) do
        @withprogress @logprogress 1.0
    end
    @test [l.message.progress.done for l in logs] == [false, false, true]
end

@testset "invalid input" begin
    local err
    @test try
        @eval @withprogress invalid_argument = "" nothing
    catch err
        err
    end isa Exception  # unfortunately `LoadError`, not an `ArgumentError`
    msg = sprint(showerror, err)
    @test occursin("Unsupported optional arguments:", msg)
    @test occursin("invalid_argument", msg)
end

module IsolatedNamespace
import ProgressLogging
using Test: @test, @testset, collect_test_logs
@testset "IsolatedNamespace" begin
    logs, = collect_test_logs(min_level = ProgressLogging.ProgressLevel) do
        ProgressLogging.@withprogress ProgressLogging.@logprogress "hello" 0.1
    end
    @test length(logs) == 3
    @test logs[1].kwargs[:progress] === NaN
    @test logs[2].kwargs[:progress] === 0.1
    @test logs[3].kwargs[:progress] === "done"
    @test logs[1].message == ""
    @test logs[2].message == "hello"
    @test logs[3].message == ""
    @test [l.message.progress.fraction for l in logs] == [nothing, 0.1, nothing]
    @test [l.message.progress.done for l in logs] == [false, false, true]
    @test length(unique([l.message.progress.id for l in logs])) == 1
    @test length(unique([l.id for l in logs])) == 1
    @test unique([l.message.progress.parentid for l in logs]) == [ProgressLogging.ROOTID]
end
end  # module IsolatedNamespace

end  # module


text:
export quaddata

"""
    quaddata(operator, test_refspace, trial_refspace, test_elements, trial_elements)

Returns an object cashing data required for the computation of boundary element
interactions. It is up to the client programmer to decide what (if any) data is
cached. For double numberical quadrature, storing the integration points for
example can significantly speed up matrix assembly.

- `operator` is an integration kernel.
- `test_refspace` and `trial_refspace` are reference space objects. `quadata`
is typically overloaded on the type of these local spaces of shape functions.
(See the implementation in `maxwell.jl` for an example).
- `test_elements` and `trial_elements` are iterable collections of the geometric
elements on which the finite element space are defined. These are provided to
allow computation of the actual integrations points - as opposed to only their
coordinates.
"""
function quaddata end


text:
using HorizonSideRobots
include("../basics.jl")

md"""
            task_4 c перегородками
"""

 """Решение задачи 17"""
 function solve_16!(r)
    visited = Set()
    (M, N), (I, J) = get_field_size_and_position!(r)
    cells = Set()

    for i in 1:M
        for j in N:-1:i
            push!(cells, (M-i+1, j))
        end
    end

    mark_some_cells_with_dfs!(r, visited, (I, J), false, cells)

    show(r)
end


text:
function Base.getindex(repo::Repo{T},id::Int) where T
    return repo.cargo[id]
end

function Base.getindex(cls::RepoClass{T},pos::Int) where T
    return cls.repo.cargo[repo_id(cls.repo,cls.c,pos)]
end

function Base.getindex(repo::Repo{T},cls::Int,pos::Int) where T
    p=repo.starts[cls]+pos-1
    id=repo.ids[p]
    return repo.cargo[id]
end

function Base.setindex!(repo::Repo{T},val::T,id::Int) where T
    repo.cargo[id]=val
end

function Base.setindex!(repo::Repo{T},val::T,c::Int,pos::Int) where T
    p=repo.starts[c]+pos-1
    id=repo.ids[p]
    repo.cargo[id]=val
end

function Base.iterate(repo::Repo{T}, n=Int(1)) where T
    if n>repo.N return nothing; end
    return (repo.cargo[n],n+1)
end

function Base.iterate(cls::RepoClass{T}, n=(0,Array{Int,1}(undef,0))) where T
    i=n[1]
    a=n[2]
    if (i==0)
        count=cls.repo.counts[cls.c]
        a=Array{Int,1}(undef,count)
        for j in 1:count
            a[j]=repo_id(cls,j)
        end
        i=1
    end
    if i>length(a) return nothing; end
    return (cls.repo.cargo[a[i]],(i+1,a))
end

function Base.eltype(::Repo{T}) where T
    return T
end

function Base.eltype(::RepoClass{T}) where T
    return T
end

function Base.size(repo::Repo{T}) where T
    return repo.N
end

function Base.size(rc::RepoClass{T}) where T
    repo=rc.repo
    return repo.counts[rc.c]
end

function Base.firstindex(repo::Repo{T}) where T
    return 1
end

function Base.firstindex(cls::RepoClass{T}) where T
    return 1
end

function Base.lastindex(repo::Repo{T}) where T
    return repo.N
end

function Base.lastindex(cls::RepoClass{T}) where T
    return cls.repo.counts[cls.c]
end

function Base.IndexStyle(repo::Repo{T}) where T
    return IndexLinear()
end

function Base.IndexStyle(cls::RepoClass{T}) where T
    return IndexLinear()
end

function Base.length(repo::Repo{T}) where T
    return repo.N
end

function Base.length(cls::RepoClass{T}) where T
    return cls.repo.counts[cls.c]
end


text:
using Documenter, MetacommunityDynamics
using DynamicGrids
using DynamicGrids: CellRule, SetCellRule, applyrule,Rule,Ruleset;

makedocs(
    sitename="MetacommunityDynamics.jl",
    authors="Michael Catchen",
    modules=[MetacommunityDynamics],
    pages=[
        "Index" => "index.md",
        "Library" => [
            "Internal" => "lib/internal.md",
            "Public API" => "lib/public.md",
        ],
        ]
)

deploydocs(
    deps=Deps.pip("pygments", "python-markdown-math"),
    repo="github.com/gottacatchenall/MetacommunityDynamics.jl.git",
    devbranch="main",
    push_preview=true
)