text:
using SymPy
using Test
using LinearAlgebra

@testset "Matrix" begin
    ## matrices
    x, y = @vars x y
    A = [x 1; 1 x]
    B = [x 1; 0 2x]
    v = [x, 2]

    ## These fail for older installations of SymPy
    @test simplify(det(A)) == x^2 - 1

    ## we use inverse for A[:inv]()

    # aliased to use inverse
    @test simplify.(inv(A) * A) ==  [1 0; 0 1]
    @test simplify.(A * inv(A)) == [1 0; 0 1]
    ##XXX    @test simplify.(A[:inv]() - inv(A)) == zeros(2, 2)

    @test SymPy.adjoint(B) == [adjoint(x) 0; 1 adjoint(2x)]
    @test SymPy.adjoint(B) == B'
    @test A.dual() == sympy.zeros(2, 2)

    A1 = Sym[25 15 -5; 15 18 0; -5 0 11]
    r = A1.cholesky()
    @test r*r.transpose() == A1

    #    s = LUsolve(A, v)
    s = A.LUsolve(B)
    @test simplify.(A * s) == B

    # norm
    @test norm(A) == sqrt(2 * abs(x)^2 + 2)
    # test norm for different subtypes of AbstractArray
    ## XXX @test norm(A) == norm(Symmetric(A)) LinearAlgebra.Symmetric no long works
    @test norm(A) == norm(view(A, :, :))

    # abs
    @test  all(convert.(Bool, abs.(A) .≧ 0))
    @test  abs.(A) == abs.(view(A, :, :))

    # is_lower, is_square, is_symmetric much slower than julia only counterparts. May deprecate, but for now they are here
    @test A.is_lower == istril(A)
    @test A.is_square == true
    @test A.is_symmetric() != issymmetric(A)

    @vars x real=true
    A = [x 1; 1 x]
    @test A.is_symmetric() == issymmetric(A)

    @test Set(eigvals(A)) == Set([x-1, x+1])

    # issue with transpose being non-typestable
    @test eltype(transpose(A)) == Sym
    @test eltype(Symmetric(A)) == Sym

    #numerical tests
    M = Sym[1 0 0; 0 1 0; 0 0 x]
    evecs = eigvecs(M)
    @test evecs[:,1] == [1, 0, 0]

    A = Sym[1 2 3; 3 6 2; 2 0 1]
    q, r = A.QRdecomposition()
    @test q * r == A
    @test abs(det(q)) == 1

    # for v0.4, the vector type is not correctly inferred
    #L = Vector{Sym}[Sym[2,3,5], Sym[3,6,2], Sym[8,3,6]]
    # L = collect(sympy.Matrix.(([[2,3,5]], [[3,6,2]], [[8,3,6]]]))
    L = collect(sympy.Matrix.([[2,3,5], [3,6,2], [8,3,6]]))
    out = sympy.GramSchmidt(L, true)  # qualify, as L not SymbolicObject type
    for i = 1:3, j = i:3
        @test out[i].dot(out[j]) == (i == j ? 1 : 0)
    end

    A = Sym[4 3; 6 3]
    L, U, _ = A.LUdecomposition()
    @test L == Sym[1 0; 3//2 1]

    A = Sym[1 0; 0 1] * 2
    B = Sym[1 2; 3 4]
    @test A.diagonal_solve(B) == B/2

    M = Sym[1 2 0; 0 3 0; 2 -4 2]
    P, D = M.diagonalize()
    @test D == [1 0 0; 0 2 0; 0 0 3]
    @test P == [-1 0 -1; 0 0 -1; 2 1  2]
    @test D == inv(P) * M * P

    # test SymPy's expm against Julia's expm
    A = [1 2 0; 0 3 0; 2 -4 2]
    M = Sym.(A)
    ## no exp(M)!
    U = M.exp() - exp(A)
    @test maximum(abs.(N.(U))) <= 1e-12

    M = [x y; 1 0]
    @test integrate.(M, x) == [x^2/2 x*y; x 0]
    @test integrate.(M, Ref((x, 0, 2))) == [2 2y; 2 0]

    M = Sym[1 3 0; -2 -6 0; 3 9 6]
    @test M.nullspace()[1] ==  reshape(Sym[-3, 1, 0], 3, 1)

    M = Sym[1 2 0; 0 3 0; 2 -4 2]
    # M.cofactor  uses 0-based indexing!
    i, j = 1, 2
    @test M.cofactor(i, j) == (-1)^(i+j) * det(M[setdiff(1:3, i+1), setdiff(1:3, j+1)])
    @test M.adjugate() / M.det() == M.inv()

    M = Sym[ 6  5 -2 -3;
            -3 -1  3  3;
             2  1 -2 -3;
            -1  1  5  5]

    P, J = M.jordan_form()
    @test J == [2 1 0 0;
                0 2 0 0;
                0 0 2 1;
                0 0 0 2]
    @test J == inv(P) * M * P

    ρ, ϕ = symbols("rho, phi")
    X = [ρ*cos(ϕ), ρ*sin(ϕ), ρ^2]
    Y = [ρ, ϕ]
    @test X.jacobian(Y) ==   [cos(ϕ) -ρ*sin(ϕ);
                              sin(ϕ)  ρ*cos(ϕ);
                                  2ρ       0]
    X = [ρ*cos(ϕ), ρ*sin(ϕ)]
    @test convert(Matrix{Sym}, X.jacobian(Y)) == [cos(ϕ) -ρ*sin(ϕ);
                                                  sin(ϕ)  ρ*cos(ϕ)]
    @test X.jacobian(Y) == view(X, :, :).jacobian(view(Y, :, :))

    ## Issue   #359
    if VERSION >= v"1.2"
        @vars a
        @test eltype([a 1; 1 a] + I) == Sym
        @test eltype([a 1; 1 a] + 2I) == Sym
    end

    ## Issue #397 adjoint losing type
    A = ones(Sym, 1, 1)
    @test eltype(A * A') == Sym

    ## Issue 413 wth matrix exponential; SymMatrix multiplication
    @vars a
    A,A1 = [1 0; 0 a], [1 1; a 2]
    @test exp(A) == [exp(Sym(1)) 0; 0 exp(a)]
    B,B1 = convert(SymMatrix,A), convert(SymMatrix,A1)
    @test B*B1 == convert(SymMatrix, A*A1)
end


text:
using BSplines
using Test
using LinearAlgebra
using SparseArrays

function vecdist(a::AbstractVector, b::AbstractVector,
                 ϵ = eps(eltype(a)))
    δ = √(sum(abs2, a-b))
    δ, δ/√(sum(abs2, a .+ ϵ))
end

@testset "Knot sets" begin
    @testset "Linear" begin
        t = LinearKnotSet(7, 0, 1, 10)
        @test order(t) == 7
        @test numintervals(t) == 10
        @test length(t) == 23
        @test first(t) == 0
        @test last(t) == 1
        for i = 1:7
            @test t[i] == 0
            @test t[end-i+1] == 1
        end
        tt = range(0, stop=1, length=11)
        for i = 7:17
            @test t[i] == tt[i-6]
        end

        ttt = collect(t)
        @test length(t) == length(ttt)
        for i ∈ eachindex(t)
            @test t[i] == ttt[i]
        end

        @test eltype(t) <: Real
    end

    @testset "Exponential" begin
        t = ExpKnotSet(7, -8, 1, 31)
        @test order(t) == 7
        @test numintervals(t) == 31
        @test first(t) == 0
        @test last(t) == 10

        t′ = ExpKnotSet(7, -8, 1, 31, include0=false)
        @test order(t′) == 7
        @test numintervals(t′) == 31
        @test first(t′) == 1e-8
        @test last(t′) == 10

        t′′ = ExpKnotSet(7, -8, 1, 31, base=2, include0=false)
        @test order(t′′) == 7
        @test numintervals(t′′) == 31
        @test first(t′′) == 1/256
        @test last(t′′) == 2
    end
end

@testset "Quadrature" begin
    t = LinearKnotSet(1, 0, 1, 2)
    x,w = BSplines.lgwt(t)
    @test all(w .== 1/4)
    @test x == [-1,1,-1,1]/(4*√3) + [1,1,3,3]/4
end

@testset "Basis" begin
    k = 3
    t = LinearKnotSet(k, 0, 1, 2)
    basis = BSplines.Basis(t)
    @testset "Eval on subintervals" begin
        x₁ = range(-1,stop=-0.5,length=10)
        Bᵢ₁ = basis(x₁)
        @test norm(Bᵢ₁) == 0

        function testbasis(x)
            Bᵢ = basis(x)
            B̃ = spzeros(Float64, length(x), 4)
            B̃[:,1] = (x .>= 0) .* (x .< 0.5) .* ((2x .- 1).^2)
            B̃[:,2] = (x .>= 0) .* (x .< 0.5) .* (2/3*(1 .- (3x .- 1).^2)) +
                (x .>= 0.5) .* (x .< 1)  .* (2*(x .- 1).^2)
            B̃[:,3] = (x .>= 0) .* (x .< 0.5) .* (2*x.^2) +
                (x .>= 0.5) .* (x .< 1)  .* (2/3*(1 .- (3x .- 2).^2))
            B̃[:,4] = (x .>= 0.5) .* (x .< 1) .* ((2x .- 1).^2)
            for j = 1:4
                δ,δr = vecdist(Bᵢ[:,j],B̃[:,j])
                @test δ < 1e-15
                @test δr < 1e-15
            end
        end
        testbasis(range(0,stop=1,length=50))
        testbasis(range(-0.5,stop=0.6,length=40))
        testbasis(range(0.5,stop=1.6,length=40))
    end

    @testset "Scalar operators" begin
        @testset "Overlap matrix" begin
            B = basis(I)
            Br = [ 1/10  7/120  1/120      0
                  7/120    1/6   1/10  1/120
                  1/120   1/10    1/6  7/120
                      0  1/120  7/120   1/10]
            @test norm(B-Br) < 1e-15
        end
        @testset "Position operator" begin
            B = basis(x->x)
            Br = [ 1/120   1/96   1/480       0
                    1/96  7/120    1/20   1/160
                   1/480   1/20  13/120  23/480
                       0  1/160  23/480  11/120]
            @test norm(B-Br) < 1e-15
        end
    end

    @testset "Derivative operators" begin
        # References generated using AnalyticBSplines.jl
        ∂n = [derop(basis, n) for n ∈ 1:k-1]
        ∂ref = [[ -1//2    5//12   1//12   0
                  -5//12   0//1    1//3   1//12
                  -1//12  -1//3    0//1   5//12
                    0     -1//12  -5//12  1//2],
                [ 4//3  -2//1   2//3   0
                  2//1  -8//3   0//1  2//3
                  2//3   0//1  -8//3  2//1
                   0     2//3  -2//1  4//3]]
        for n ∈ 1:k-1
            @test norm(∂n[n]-∂ref[n]) < 1e-15
        end

        k′ = 7
        t′ = LinearKnotSet(k′, 0, 1, 3)
        basis′ = BSplines.Basis(t′)

        ∂′n = [derop(basis′, n) for n ∈ 1:k′-1]
        ∂′ref = [
 [    -1//2          10625//29568     283021//2395008    23827//1197504     479//199584        5//27216          1//149688         0                0
  -10625//29568          0//1          35867//177408    118891//1064448   26167//709632    11017//1419264     7141//8515584       1//946176         0
 -283021//2395008   -35867//177408         0//1           5129//39424     11645//101376   245743//4257792    20455//1216512    7141//8515584       1//149688
  -23827//1197504  -118891//1064448    -5129//39424          0//1         11035//118272   218689//2128896   245743//4257792   11017//1419264       5//27216
    -479//199584    -26167//709632    -11645//101376    -11035//118272        0//1         11035//118272     11645//101376    26167//709632      479//199584
      -5//27216     -11017//1419264  -245743//4257792  -218689//2128896  -11035//118272        0//1           5129//39424    118891//1064448   23827//1197504
      -1//149688     -7141//8515584   -20455//1216512  -245743//4257792  -11645//101376    -5129//39424          0//1         35867//177408   283021//2395008
        0               -1//946176     -7141//8515584   -11017//1419264  -26167//709632  -118891//1064448   -35867//177408        0//1         10625//29568
        0                 0               -1//149688        -5//27216      -479//199584   -23827//1197504  -283021//2395008  -10625//29568         1//2],
[    90//11      -27843//2464      169651//66528      2347//4752        377//5544         7//1188           1//4158            0                0
  16509//2464     -4413//616      -104639//133056   191011//266112    72047//177408  118771//1064448    31501//2128896        3//78848          0
 169651//66528  -104639//133056    -46861//22176    -16405//29568     58403//177408   45935//118272    121843//709632     31501//2128896       1//4158
   2347//4752    191011//266112    -16405//29568     -3083//3168     -36415//88704    39203//177408     45935//118272    118771//1064448       7//1188
    377//5544     72047//177408     58403//177408   -36415//88704     -5807//7392    -36415//88704      58403//177408     72047//177408      377//5544
      7//1188    118771//1064448    45935//118272    39203//177408   -36415//88704    -3083//3168      -16405//29568     191011//266112     2347//4752
      1//4158     31501//2128896   121843//709632    45935//118272    58403//177408  -16405//29568     -46861//22176    -104639//133056   169651//66528
       0              3//78848      31501//2128896  118771//1064448   72047//177408  191011//266112   -104639//133056     -4413//616       16509//2464
       0               0                1//4158          7//1188        377//5544      2347//4752      169651//66528     -27843//2464         90//11],
[  -108//1     41193//224    -178259//2016    10723//1008     41//24        43//252        1//126          0             0
 -23049//224     162//1       -78221//1344   -13759//2688   4773//1792    2023//1536    5039//21504       9//7168        0
 -93901//2016  78221//1344         0//1       -9741//896    -905//256     4079//3584   10065//7168     5039//21504      1//126
 -10723//1008  13759//2688      9741//896         0//1     -4621//896    -5055//1792    4079//3584     2023//1536      43//252
    -41//24    -4773//1792       905//256      4621//896       0//1      -4621//896     -905//256      4773//1792      41//24
    -43//252   -2023//1536     -4079//3584     5055//1792   4621//896        0//1      -9741//896    -13759//2688   10723//1008
     -1//126   -5039//21504   -10065//7168    -4079//3584    905//256     9741//896        0//1      -78221//1344   93901//2016
       0          -9//7168     -5039//21504   -2023//1536  -4773//1792   13759//2688   78221//1344     -162//1      23049//224
       0            0             -1//126       -43//252     -41//24    -10723//1008  178259//2016   -41193//224      108//1],
[  1080//1    -224775//112    413855//336     -58255//168     1045//28       185//42         5//21            0              0
 138105//112   -62505//28     851525//672    -354385//1344  -16085//896    56255//5376   34625//10752      135//3584         0
 232415//336  -781435//672     57055//112       7365//448   -45905//896   -25215//1792   25535//3584     34625//10752       5//21
  32465//168  -354385//1344     7365//448       8255//112      595//64     -3815//128   -25215//1792     56255//5376      185//42
   1045//28    -16085//896    -45905//896        595//64      5055//112      595//64    -45905//896     -16085//896      1045//28
    185//42     56255//5376   -25215//1792     -3815//128      595//64      8255//112     7365//448    -354385//1344    32465//168
      5//21     34625//10752   25535//3584    -25215//1792  -45905//896     7365//448    57055//112    -781435//672    232415//336
       0          135//3584    34625//10752    56255//5376  -16085//896  -354385//1344  851525//672     -62505//28     138105//112
       0             0             5//21         185//42      1045//28    -58255//168   413855//336    -224775//112      1080//1],
[   -7290//1    1578285//112   -1077255//112    203895//56     -26055//28       1395//14         45//7             0              0
 -1170045//112    40095//2     -3001995//224   2125575//448   -866295//896     -3375//256    127305//3584      3645//3584         0
  -873225//112  3274155//224     -18225//2     1109295//448     26325//896   -421605//1792   -10125//512     127305//3584       45//7
  -158535//56   2228985//448   -1109295//448         0//1      202095//448     26325//896   -421605//1792     -3375//256      1395//14
   -19305//28    866295//896     -26325//896   -202095//448         0//1      202095//448     26325//896    -866295//896     19305//28
    -1395//14      3375//256     421605//1792   -26325//896   -202095//448         0//1     1109295//448   -2228985//448    158535//56
      -45//7    -127305//3584     10125//512    421605//1792   -26325//896  -1109295//448     18225//2     -3274155//224    873225//112
         0        -3645//3584   -127305//3584     3375//256    866295//896  -2125575//448   3001995//224     -40095//2     1170045//112
         0             0            -45//7       -1395//14      26055//28    -203895//56    1077255//112   -1578285//112      7290//1],
[174960//7    -98415//2      70605//2     -107595//7      38070//7       -9180//7       1080//7            0            0
  98415//2   -677970//7    1944135//28    -241245//8    1185435//112    -79245//32    105435//448     10935//448        0
 115155//2  -3158865//28   1117395//14   -1852875//56    157545//16    -176175//224  -419985//448    105435//448    1080//7
 215595//7   -474525//8    2229525//56    -186705//14     26325//56     325215//112  -176175//224    -79245//32    13500//7
  72090//7  -2080485//112   157545//16      26325//56    -57105//14      26325//56    157545//16   -2080485//112   72090//7
  13500//7    -79245//32   -176175//224    325215//112    26325//56    -186705//14   2229525//56    -474525//8    215595//7
   1080//7    105435//448  -419985//448   -176175//224   157545//16   -1852875//56   1117395//14   -3158865//28   115155//2
       0       10935//448   105435//448    -79245//32   1185435//112   -241245//8    1944135//28    -677970//7     98415//2
       0            0         1080//7       -9180//7      38070//7     -107595//7      70605//2      -98415//2    174960//7]]

        for n ∈ 1:k′-1
            @test norm(∂′n[n]-∂′ref[n])/abs(1e-16 + norm(∂′ref[n])) < 1e-14
        end
    end

    @testset "Dirichlet0 conditions" begin
        k = 3
        t = LinearKnotSet(k, 0, 2, 4)
        basis = BSplines.Basis(t, 0, 0)

        ∂n = [derop(basis, n) for n ∈ 1:k-1]
        ∂ref = [[ 0    0       0       0      0     0
                  0   0//1    3//8    1//24   0     0
                  0  -3//8    0//1    5//12  1//24  0
                  0  -1//24  -5//12   0//1   3//8   0
                  0    0     -1//24  -3//8   0//1   0
                  0    0       0       0      0     0],
                [ 0    0      0      0      0    0
                  0  -8//3   1//3   1//3    0    0
                  0   1//3  -2//1   2//3   1//3  0
                  0   1//3   2//3  -2//1   1//3  0
                  0    0     1//3   1//3  -8//3  0
                  0    0      0      0      0    0]]
        for n ∈ 1:k-1
            @test norm(∂n[n]-∂ref[n]) < 1e-15
        end
    end
end

@testset "Splines" begin
    k = 7
    t = LinearKnotSet(k, 0, 1, 10)
    basis = BSplines.Basis(t)
    x = range(first(t), stop=last(t), length=101)[1:end-1] # Due to issue #2
    B = basis(x)
    for n = 1:k-1
        c = rand(n)
        f = x -> sum(c[i]*x^(i-1) for i in eachindex(c))
        S = Spline(f, basis)
        @test norm(S(B) - f.(x))/abs(1e-16 + norm(f.(x))) < 1e-14
    end
end

ext: jl


text:
using DifferentialEquations, Plots
A = 1
B = 1
C = 4
x0 = 1
function f(du, u, p, t)
    du[1] = -0.5 * u[1]^2 + 0.5 * C^2 * u[2] / B
    du[2] = -A + u[1] * u[2]
end

function bc1!(residual, u, p, t)
    residual[1] = u[1][1] - x0
    residual[2] = u[end][2]
end
tspan = (0,1)
bvp1 = BVProblem(f, bc1!, [0, 0], tspan)
sol1 = solve(bvp1, GeneralMIRK4(), dt = 0.05)
plot(sol1)


text:
"""
    createexperiment(mlf::MLFlow; name=missing, artifact_location=missing, tags=missing)

Creates an MLFlow experiment.

# Arguments
- `mlf`: [`MLFlow`](@ref) configuration.
- `name`: experiment name. If not specified, MLFlow sets it.
- `artifact_location`: directory where artifacts of this experiment will be stored. If not specified, MLFlow uses its default configuration.
- `tags`: a Dictionary of key-values which tag the experiment.

# Returns
An object of type [`MLFlowExperiment`](@ref).

"""
function createexperiment(mlf::MLFlow; name=missing, artifact_location=missing, tags=missing)
    endpoint = "experiments/create"
    if ismissing(name)
        name = string(UUIDs.uuid4())
    end
    result = mlfpost(mlf, endpoint; name=name, artifact_location=artifact_location, tags=tags)
    experiment_id = parse(Int, result["experiment_id"])
    getexperiment(mlf, experiment_id)
end

"""
    getexperiment(mlf::MLFlow, experiment_id::Integer)

Retrieves an MLFlow experiment by id.

# Arguments
- `mlf`: [`MLFlow`](@ref) configuration.
- `experiment_id`: Experiment identifier.

# Returns
An instance of type [`MLFlowExperiment`](@ref)

"""
function getexperiment(mlf::MLFlow, experiment_id::Integer)
    try
        result = _getexperimentbyid(mlf, experiment_id)
        return MLFlowExperiment(result)
    catch e
        if isa(e, HTTP.ExceptionRequest.StatusError) && e.status == 404
            return missing
        end
        throw(e)
    end
end
"""
    getexperiment(mlf::MLFlow, experiment_name::String)

Retrieves an MLFlow experiment by name.

# Arguments
- `mlf`: [`MLFlow`](@ref) configuration.
- `experiment_name`: Experiment name.

# Returns
An instance of type [`MLFlowExperiment`](@ref)

"""
function getexperiment(mlf::MLFlow, experiment_name::String)
    try
        result = _getexperimentbyname(mlf, experiment_name)
        return MLFlowExperiment(result)
    catch e
        if isa(e, HTTP.ExceptionRequest.StatusError) && e.status == 404
            return missing
        end
        throw(e)
    end
end
function _getexperimentbyid(mlf::MLFlow, experiment_id::Integer)
    endpoint = "experiments/get"
    arguments = (:experiment_id => experiment_id, )
    mlfget(mlf, endpoint; arguments...)["experiment"]
end
function _getexperimentbyname(mlf::MLFlow, experiment_name::String)
    endpoint = "experiments/get-by-name"
    arguments = (:experiment_name => experiment_name, )
    mlfget(mlf, endpoint; arguments...)["experiment"]
end

"""
    getorcreateexperiment(mlf::MLFlow, experiment_name::String)

Gets an experiment if one alrady exists, or creates a new one.

# Arguments
- `mlf`: [`MLFlow`](@ref) configuration.
- `experiment_name`: Experiment name.

# Returns
An instance of type [`MLFlowExperiment`](@ref)

"""
function getorcreateexperiment(mlf::MLFlow, experiment_name::String)
    exp = getexperiment(mlf, experiment_name)
    if ismissing(exp)
        exp = createexperiment(mlf, name=experiment_name)
    end
    exp
end

"""
    deleteexperiment(mlf::MLFlow, experiment_id::Integer)

Deletes an MLFlow experiment.

# Arguments
- `mlf`: [`MLFlow`](@ref) configuration.
- `experiment_id`: experiment identifier.

# Returns

`true` if successful. Otherwise, raises exception.
"""
function deleteexperiment(mlf::MLFlow, experiment_id::Integer)
    endpoint = "experiments/delete"
    try
        result = mlfpost(mlf, endpoint; experiment_id=experiment_id)
    catch e
        if isa(e, HTTP.ExceptionRequest.StatusError) && e.status == 404
            # experiment already deleted
            return true
        end
        throw(e)
    end
    true
end
"""
    deleteexperiment(mlf::MLFlow, experiment::MLFlowExperiment)

Deletes an MLFlow experiment.

# Arguments
- `mlf`: [`MLFlow`](@ref) configuration.
- `experiment`: an object of type [`MLFlowExperiment`](@ref)

Dispatches to `deleteexperiment(mlf::MLFlow, experiment_id::Integer)`.

"""
deleteexperiment(mlf::MLFlow, experiment::MLFlowExperiment) =
    deleteexperiment(mlf, experiment.experiment_id)

"""
    listexperiments(mlf::MLFlow)

Returns a list of MLFlow experiments.

TODO: not yet entirely implemented
"""
function listexperiments(mlf::MLFlow)
    endpoint = "experiments/list"
    mlfget(mlf, endpoint)
end

max_stars_repo_name: JuliaAI/MLFlowClient

max_issues_repo_name: JuliaAI/MLFlowClient

max_forks_repo_name: JuliaAI/MLFlowClient


text:
# create a new 'type'
struct Person
    name::String
    age::Int
    gender::Char
end

p1 = Person("Richard", 21, 'M')
println("I have a person : $(p1.name) / age : $(p1.age) / gender: $(p1.gender)")

# types are 'immutable' by default. To ensure that you can update attributes of a type use the 'mutable' keyword
mutable struct Car
    manufacturer::String
    model::String;
    registered::Bool
end
c1 = Car("Ford", "Escort", true)
println(c1)
c1.registered = false
println(c1)

# to check the 'type' of a variable use the 'isa' method
println(isa(c1, Car))

max_stars_repo_name: rpillar/JuliaExamples

max_issues_repo_name: rpillar/JuliaExamples

max_forks_repo_name: rpillar/JuliaExamples


text:
#Pkg.checkout("EcoEvo3D","model1")
using EcoEvo3D

function main(ngenana_,ngenretana_)
	seed = 177;#seed for random numbers (to control the outputs)
	nreal = 2;#Number of realizations
	Gmax = 2000;#Maximum number of Generations
	J = 50000;#Total Number of individuals in the system
        nGenAna = ngenana_
	anaG = J*nGenAna;#Threshold to consider speciation (anagenesis and cladogenesis)
        nGenRetG = ngenretana_; #parameter to define gene flow retard in anagenetic speciation (10% of retard, minimum value of retard is 1)
	retG = nGenRetG*J;#gene flow retard in anagenetic speciation

	srand(seed);
	v =  (0.000005 + rand()* 0.00005);#Cladogenesis speciation
	mr = (0.00001 + rand()* 0.0001);#Regional migration -
	ml = (0.008 + rand()* 0.008);#Local Migration

	distmatfile = "upstream_cost0.001_size8_m1.txt";#Name of the file that contains the location of the points with clownfishes.
	verticesdata = "VerticesData.txt";#Name of the file that contains the size and the height of the points with clownfishes.
#	verticesdata = "vert.txt";#Name of the file that contains the size and the height of the points with clownfishes.
	model = 1;

	EcoEvo3D.dynamic(seed,nreal,Gmax,J,v,mr,ml,anaG,retG,distmatfile,verticesdata,model);
end

if(length(ARGS)==0)
    println();
    println("Please give the correct parameters:");
    println("To run: julia test.jl <ngenana> <ngenretana>");
    println();
    println("where: ");
    println("- ngenana: number of generations for anagenesis speciation");
    println("- ngenretana: number of generations in the retard of an anagenesis speciation");
    println();
else
    main(float(ARGS[1]),float(ARGS[2]))
end


text:
# https://www.codewars.com/kata/556deca17c58da83c00002db

function tribonacci((a, b, c), n::Int)::Array{Float64}
  if n <= 3
    [a, b, c][1:n]
  else
    d = a + b + c
    append!([a], tribonacci((b, c, d), n-1))
  end
end

max_stars_repo_name: hiljusti/codewars-solutions

max_issues_repo_name: hiljusti/codewars-solutions

max_forks_repo_name: hiljusti/codewars-solutions


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""
    ObjsPrimaryOwner(;
        email=nothing,
        id=nothing,
    )

    - email::String
    - id::String
"""
mutable struct ObjsPrimaryOwner <: SwaggerModel
    email::Any # spec type: Union{ Nothing, String } # spec name: email
    id::Any # spec type: Union{ Nothing, String } # spec name: id

    function ObjsPrimaryOwner(;email=nothing, id=nothing)
        o = new()
        validate_property(ObjsPrimaryOwner, Symbol("email"), email)
        setfield!(o, Symbol("email"), email)
        validate_property(ObjsPrimaryOwner, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        o
    end
end # type ObjsPrimaryOwner

const _property_map_ObjsPrimaryOwner = Dict{Symbol,Symbol}(Symbol("email")=>Symbol("email"), Symbol("id")=>Symbol("id"))
const _property_types_ObjsPrimaryOwner = Dict{Symbol,String}(Symbol("email")=>"String", Symbol("id")=>"String")
Base.propertynames(::Type{ ObjsPrimaryOwner }) = collect(keys(_property_map_ObjsPrimaryOwner))
Swagger.property_type(::Type{ ObjsPrimaryOwner }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ObjsPrimaryOwner[name]))}
Swagger.field_name(::Type{ ObjsPrimaryOwner }, property_name::Symbol) =  _property_map_ObjsPrimaryOwner[property_name]

function check_required(o::ObjsPrimaryOwner)
    (getproperty(o, Symbol("email")) === nothing) && (return false)
    (getproperty(o, Symbol("id")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ObjsPrimaryOwner }, name::Symbol, val)
end


text:
module Bridges

using Compat

using MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities

const MOIB = MOI.Bridges # used in macro

const SVF = MOI.SingleVariable
const VVF = MOI.VectorOfVariables
const SAF{T} = MOI.ScalarAffineFunction{T}
const VAF{T} = MOI.VectorAffineFunction{T}
const SQF{T} = MOI.ScalarQuadraticFunction{T}
const VQF{T} = MOI.VectorQuadraticFunction{T}

const VI = MOI.VariableIndex
const CI = MOI.ConstraintIndex

include("bridge.jl")
include("bridgeoptimizer.jl")
include("singlebridgeoptimizer.jl")
include("lazybridgeoptimizer.jl")

# This is used by JuMP and removes the need to update JuMP everytime a bridge is added
MOIU.@model(AllBridgedConstraints,
            (),
            (MOI.Interval,),
            (MOI.SecondOrderCone, MOI.RotatedSecondOrderCone, MOI.GeometricMeanCone,
             MOI.PositiveSemidefiniteConeSquare,
             MOI.LogDetConeTriangle, MOI.RootDetConeTriangle),
            (),
            (MOI.SingleVariable,),
            (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
            (MOI.VectorOfVariables,),
            (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction))
"""
    fullbridgeoptimizer(model::MOI.ModelLike, ::Type{T}) where T

Returns a `LazyBridgeOptimizer` bridging `model` for every bridge defined in this package and for the coefficient type `T`.
"""
function fullbridgeoptimizer(model::MOI.ModelLike, ::Type{T}) where T
    bridgedmodel = MOIB.LazyBridgeOptimizer(model, AllBridgedConstraints{T}())
    add_bridge(bridgedmodel, MOIB.SplitIntervalBridge{T})
    add_bridge(bridgedmodel, MOIB.QuadtoSOCBridge{T})
    add_bridge(bridgedmodel, MOIB.GeoMeanBridge{T})
    add_bridge(bridgedmodel, MOIB.SquarePSDBridge{T})
    add_bridge(bridgedmodel, MOIB.LogDetBridge{T})
    add_bridge(bridgedmodel, MOIB.RootDetBridge{T})
    add_bridge(bridgedmodel, MOIB.RSOCBridge{T})
    add_bridge(bridgedmodel, MOIB.RSOCtoPSDBridge{T})
    add_bridge(bridgedmodel, MOIB.SOCtoPSDBridge{T})
    bridgedmodel
end

include("intervalbridge.jl")
@bridge SplitInterval SplitIntervalBridge () (MOI.Interval,) () () (MOI.SingleVariable,) (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction) () ()
include("rsocbridge.jl")
@bridge RSOC RSOCBridge () () (MOI.RotatedSecondOrderCone,) () () () (MOI.VectorOfVariables,) (MOI.VectorAffineFunction,)
include("quadtosocbridge.jl")
@bridge QuadtoSOC QuadtoSOCBridge () (MOI.LessThan, MOI.GreaterThan) () () () (MOI.ScalarQuadraticFunction,) () ()
include("geomeanbridge.jl")
@bridge GeoMean GeoMeanBridge () () (MOI.GeometricMeanCone,) () () () (MOI.VectorOfVariables,) (MOI.VectorAffineFunction,)
include("squarepsdbridge.jl")
@bridge SquarePSD SquarePSDBridge () () (MOI.PositiveSemidefiniteConeSquare,) () () () (MOI.VectorOfVariables,) (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
include("detbridge.jl")
@bridge LogDet LogDetBridge () () (MOI.LogDetConeTriangle,) () () () (MOI.VectorOfVariables,) (MOI.VectorAffineFunction,)
@bridge RootDet RootDetBridge () () (MOI.RootDetConeTriangle,) () () () (MOI.VectorOfVariables,) (MOI.VectorAffineFunction,)
include("soctopsdbridge.jl")
@bridge SOCtoPSD SOCtoPSDBridge () () (MOI.SecondOrderCone,) () () () (MOI.VectorOfVariables,) (MOI.VectorAffineFunction,)
@bridge RSOCtoPSD RSOCtoPSDBridge () () (MOI.RotatedSecondOrderCone,) () () () (MOI.VectorOfVariables,) (MOI.VectorAffineFunction,)

end # module


text:
using Dates
using ArgParse

function parse_commandline()
    s = ArgParseSettings(
        fromfile_prefix_chars = "@",
    )

    @add_arg_table! s begin
        "--input", "-i"
            help = "one or more settings files"
            nargs = '+' # At least one input is required
            arg_type = String
            required = true
        "--output", "-o"
            help = "output folder suffix; timestamp will be prepended"
            arg_type = String
            required = true
        "--julia"
            help = "path to julia binary"
            nargs = '*' # Zero or more inputs
            arg_type = String
            default = String["julia"]
        "--decaes-version", "-v"
            help = "tag (e.g. 0.3 or v0.3), branch (e.g. master), or git commit (e.g. 338795e)"
            nargs = '*' # Zero or more inputs
            arg_type = String
            default = String["master"]
        "--threads", "-t"
            help = "number of threads"
            nargs = '*' # Zero or more inputs
            arg_type = Any
            default = Any["auto"]
        "--optimize", "-O"
            help = "optimization level"
            arg_type = Int
            nargs = '*' # Zero or more inputs
            default = Int[2]
        "--warmup"
            help = "number of warmup runs"
            arg_type = Int
            default = 1
        "--min-runs"
            help = "minimum number of runs"
            arg_type = Int
            default = 3
        # "--show-output"
        #     help = "print stdout and stderr of the benchmark instead of suppressing"
        #     action = :store_true
    end

    return parse_args(s)
end

function main()
    args = parse_commandline()
    dateformat = "yyyy-mm-dd-T-HH-MM-SS"
    timestamp = Dates.format(Dates.now(), dateformat)
    outfolder = mkpath(joinpath(dirname(args["output"]), timestamp * "_" * basename(args["output"])))
    outpath(xs...) = joinpath(outfolder, xs...)

    # Benchmarking command
    cmd = `
    hyperfine
        "JULIA_NUM_THREADS={threads} {julia} --project=$(joinpath(@__DIR__, ".bench.tmp"))/{julia}/{version} --startup-file=no --quiet --optimize={optimize} -e 'using DECAES; main()' -- @{input} --quiet"
        --prepare "{julia} --startup-file=no --project=$(@__DIR__) $(joinpath(@__DIR__, "prepare.jl")) {julia} {version}"
        --warmup $(args["warmup"])
        --min-runs $(args["min-runs"])
        --parameter-list version $(join(args["decaes-version"], ","))
        --parameter-list julia $(join(args["julia"], ","))
        --parameter-list threads $(join(args["threads"], ","))
        --parameter-list optimize $(join(args["optimize"], ","))
        --parameter-list input $(join(args["input"], ","))
        --export-markdown $(outpath("results.md"))
        --export-json $(outpath("results.json"))
    `
    # --show-output

    # Save benchmarking settings/files for future reference
    open(outpath("settings.txt"); write = true) do io
        for arg in ARGS
            startswith(arg, "@") ?
                println(io, readchomp(arg[2:end])) :
                println(io, arg)
        end
    end
    open(outpath("run_benchmarks.jl"); write = true) do io
        println(io, "mkpath(\"$(outpath())\")")
        println(io, "run($cmd)")
    end
    for jl in filter(endswith(".jl"), readdir(@__DIR__; join = true))
        cp(jl, outpath(basename(jl)); force = true)
    end

    # Run benchmarks
    @info "Benchmarking DECAES with settings:\n" * readchomp(outpath("settings.txt"))
    run(cmd)
    rm(joinpath(@__DIR__, ".bench.tmp"); force = true, recursive = true)

    return nothing
end

main()


text:
# Max function

"""
  Maximum(λ::Real=1.0)

Returns the function `g(x) = λ⋅maximum(x)`, for a nonnegative parameter `λ ⩾ 0`.
"""

Maximum{R <: Real}(lambda::R=1.0) = SumLargest(one(Int32), lambda)


text:
default_cddlib_backend(::Type{<:AbstractFloat}) = CDDLib.Library()
default_cddlib_backend(::Type{<:Rational}) = CDDLib.Library(:exact)


text:
#
Random.seed!(123)
using OrdinaryDiffEq, DiffEqOperators, LinearAlgebra

"""
GOALS

1. have defaultlinsolve use predefined ldiv! if available on that type
    - allow for `Diagonal`, `LinearAlgebra.I`

2. add misc shit in each Krylov iteration call:
    https://github.com/SciML/DiffEqBase.jl/issues/366

3. Make SEM.jl operators AbstractDiffEqOperators
4. Make SEM fields subsets of AbstractVector

"""

n = 10
x = Array(range(start=0,stop=1,length=n))
A = 0.01*Tridiagonal(-ones(n-1),2ones(n),-ones(n-1))
rn = (du, u, p, t) -> begin
    mul!(du, A, u)
end
u0 = rand(n)
prob = ODEProblem(ODEFunction(rn, jac_prototype=JacVecOperator{Float64}(rn, u0; autodiff=false)), u0, (0, 10.))
#=============================================#
# DiffEqBase/linear_nonlinear.jl
#=============================================#
"""
 OrdinaryDiffEq/test/interface/linear_nonlinear_tests.jl
"""

""" default """
#sol = solve(prob, QNDF(autodiff=false)) # linsolve=DefaultLinSolve()
#@show sol.retcode
""" GMRES """
#sol = solve(prob, QNDF(autodiff=false, linsolve=LinSolveGMRES()))
#@show sol.retcode
""" CG """
#sol = solve(prob, QNDF(autodiff=false, linsolve=LinSolveCG()))
#@show sol.retcode
#=============================================#
# LinearSolve.jl
#=============================================#

#=============================================#
"""
 OrdinaryDiffEq/test/interface/mass_matrix_tests.jl

 how to provide inverse to mass_matrix. where is mass_matrix being inverted???
 how to set preconditioner,
"""

#---------------------------------------------#
return


text:
using OMEinsum, Cracker
using Cracker: TrackedArrayType

for EC in [:DynamicEinCode, :StaticEinCode]
    @eval function OMEinsum.einsum(code::$EC, @nospecialize(xs::NTuple{N,TrackedArrayType} where N), size_dict::Dict)
        Cracker.trace(einsum, code, xs, size_dict)
    end
end


text:
## Usage: @assume x ~ Dist where x is a symbol to be used and Dist is
##   a valid distribution from the Distributions package
macro assume(ex)
  @assert ex.args[1] == symbol("@~")
  esc(quote
    $(ex.args[2]) = Turing.assume(Turing.sampler, $(ex.args[3]))
  end)
end

## Usage: @observe(x ~ Dist) where x is a value and Dist is a valid distribution
macro observe(ex)
  @assert ex.args[1] == symbol("@~")
  global TURING
  ex2 = Expr(:block, nothing)
  push!( ex2.args, :(Turing.observe(Turing.sampler,
                      logpdf($(ex.args[3]), $(ex.args[2])))) )
  esc(ex2)
end

## Usage: @predict x y z  where x, y, z are symbols
macro predict(ex...)
  ex_funcs = Expr(:block)
  for i = 1:length(ex)
    @assert typeof(ex[i]) == Symbol
    sym = string(ex[i])
    push!( ex_funcs.args, :(ct = current_task();
            Turing.predict(
              Turing.sampler, symbol($sym), get(ct, $(ex[i])))) )
  end
  # println(ex_funcs)
  esc(ex_funcs)
end

## Usage: @model f body
macro model(name, fbody)
  # functions defined via model macro have an implicit varinfo array.
  #   this varinfo array is useful is task cloning.

  fname = isa(name, Symbol) ? Expr(:call, name) : name # Turn f into f() if necessary.
  ex = Expr(:function, fname, fbody)

  push!(ex.args[2].args, :(produce(Val{:done})))
  TURING[:modelex] = ex
  return esc(ex)
end


text:
# Process generation
# Power spectrum
# Autocovariance function

function Base.rand(rng::AbstractRNG,proc::TimeSeriesModel,N::Integer,Δt_av=1.0,Tmin=0.01)
  t = randtimes(rng,Δt_av,N,Tmin)
  return rand(rng,proc,t)
end

function Base.rand(rng::AbstractRNG,e::SdeEngine,state::State,t::AbstractVector)
  eq = e.eq
  N = length(t)
  p = statedim(eq)
  if measdim(eq)≠1 error("Only 1-D observations are supported") end
  data = Vector{Obs}(undef,N)

  z = rand(rng,state)
  μy,σ2y = obsdist(e,State(z))
  yc = σ2y==0 ? μy : rand(rng,Normal(μy,√σ2y))
  data[1] = Obs(t[1],yc)

  for i in 2:N
    Δt = t[i]-t[i-1]
    state = predict(e,State(z),Δt)
    z = rand(rng,state)
    μy,σ2y = obsdist(e,State(z))
    yc = σ2y==0 ? μy : rand(rng,Normal(μy,√σ2y))
    data[i] = Obs(t[i],yc)
  end
  return data
end

function Base.rand(rng::AbstractRNG,proc::TimeSeriesModel,t::AbstractVector)
  Σ = datacov(proc,t)
  L = cholesky(Σ).L
  ϵ = randn(rng,length(t)) # Not using MvNormal because it does not accept `Σ::SymmetricToeplitz`
  y = L*ϵ
  data = [Obs(tc,yc) for (tc,yc) in zip(t,y)]
  return data
end

"
$SIGNATURES

Default Kalman engine for process `proc`
"
defaultengine(proc::SdeProcess{<:Any,<:Any}) = ElimQ(proc.eq)

defaultengine(proc::SdeProcess{<:Any,<:SdeRoots}) = SdeRootsEngine(proc.eq)

"""
    rand(rng,proc,t) -> Vector(VectorObservation)

Random draw from the stationary vector SDE process `proc` at times `t`.
"""
function Base.rand(rng::AbstractRNG,proc::SdeProcess,t::AbstractVector)
  proc′ = transform_ss(Sde,proc)
  e = ElimQ(proc′.eq)
  return rand(rng,e,proc′.state,t)
end

"""
$SIGNATURES

Autocovariance function of state `z` at lags `τ`
"""
function autocovz(eq::VectorSde,τ)
  Rz0 = statcovz(eq)
  d = statedim(eq)
  Rz = similar(Rz0,d,d,length(τ))
  for (i,tc) in enumerate(τ)
    Rz[:,:,i] = Rz0*exp(eq.A*tc)'
  end
  Rz
end

"""
$SIGNATURES

Autocovariance function of observations `y` at lags `τ`
"""
function autocov(eq::VectorSde,τ)
  d = measdim(eq)
  Ry = similar(eq.R,d,d,length(τ))
  Rz = autocovz(eq,τ)
  for i in 1:length(τ)
    Ry[:,:,i] = eq.C*Rz[:,:,i]*eq.C'+eq.R
  end
  Ry
end

autocov(eq::Sde,τ) = autocov(convert(VectorSde,eq),τ)[:]

autocov(eq::SdeRoots,τ) = real.(autocov(convert(VectorSde,eq),τ)[:])

autocov(proc::SdeProcess,τ) = autocov(proc.eq,τ)

autocov(m::WhiteNoise,τ) = [τc==0 ? m.σ^2 : zero(m.σ) for τc in τ]

powerspectrum(m::WhiteNoise,f) = zero(f)

"""
$SIGNATURES

Power spectrum of Stochastic Differential Equation `eq`
at frequencies `f`
"""
function powerspectrum(eq::VectorSde,f)
  N = length(f)
  p = statedim(eq)
  q = measdim(eq)
  S = zeros(Complex{Float64},q,q,N)
  for (i,fc) in enumerate(f)
    ω = 2π*fc
    Sz = (eq.A-im*ω*I)\eq.Q/(eq.A-im*ω*I)'
    S[:,:,i] = eq.C*Sz*eq.C'
  end
  S
end

function powerspectrum(eq::ScalarSde,f)
  mv = convert(VectorSde,eq)
  S = powerspectrum(mv,f)
  real.(S[:])
end

powerspectrum(proc::ScalarSdeProcess,f) = powerspectrum(proc.eq,f)

function Distributions.logpdf(proc::TimeSeriesModel,data)
  t = getfield.(data,:t); y = getfield.(data,:y)
  Σ = datacov(proc,t)
  N = length(data)
  lp = -1/2*logdet(Σ) - 1/2*dot(y,Σ,y) - N/2*log(2π)
  return lp
end


text:
using Interface
using Test

@testset "Interface.jl" begin
    # Write your tests here.
end


text:
function expected_hungarian!(server::RoutingAllocation, routing_sim::RoutingSimulator,
                              rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}

    available_drones = String[]
    for (drone_nm, dp) in server.agent_prop_set
        if server.agent_prop_set[drone_nm].at_depot == true
            push!(available_drones, drone_nm)
        end
    end

    available_pkgs = collect(keys(routing_sim.active_packages))

    n_drones = length(available_drones)
    n_pkgs = length(available_pkgs)

    if n_drones == 0 || n_pkgs == 0
        return
    end

    cost_matrix = ones(n_drones, n_pkgs)

    # Loop over drones X packages and assign cost matrix entries
    for (i, drone_nm) in enumerate(available_drones)
        for (j, pkg_nm) in enumerate(available_pkgs)

            depot_loc = server.agent_set[drone_nm].depot_loc
            pkg_loc = routing_sim.active_packages[pkg_nm].delivery

            if Distances.evaluate(EuclideanLatLongMetric(), convert_to_vector(depot_loc), convert_to_vector(pkg_loc)) <= routing_sim.distance_thresh

                # Compute the cost based on windows
                drone_pkg_window = server.agent_task_windows[(drone_nm, pkg_nm)]

                travel_time = drone_pkg_window[3] - drone_pkg_window[2] # SUCCESS - FINISH
                meanval = travel_time
                stdval = meanval/routing_sim.tt_est_std_scale
                tt_dist = Distributions.Epanechnikov(meanval, stdval)
                succ_prob = Distributions.cdf(tt_dist, drone_pkg_window[2])

                cost_matrix[i, j] = -succ_prob*routing_sim.delivery_reward

            end # if distance <= thresh
        end # for enumerate(available_pkgs)
    end # for enumerate(available_drones)

    if minimum(cost_matrix) == 1.0
        @warn "All matchings are failures"
        return
    end

    assignment, cost = hungarian(cost_matrix)

    for (i, assgn) in enumerate(assignment)

        if assgn != 0 && cost_matrix[i,assgn] < 1.0

            drone_nm = available_drones[i]
            pkg_nm = available_pkgs[assgn]

            server.agent_task_allocation[drone_nm] = (name=pkg_nm, time=Inf)

            routing_sim.true_delivery_return[(drone_nm, pkg_nm)] = sample_true_delivery_return_time(server.agent_task_windows[(drone_nm, pkg_nm)],
                                                                                server.current_time,
                                                                                routing_sim.tt_est_std_scale,
                                                                                rng)
            server.agent_prop_set[drone_nm].at_depot = false
            server.agent_prop_set[drone_nm].current_package = pkg_nm

            new_busy_package = routing_sim.active_packages[pkg_nm]
            routing_sim.busy_packages[pkg_nm] = new_busy_package
            delete!(routing_sim.active_packages, pkg_nm)
            routing_sim.num_active_packages -= 1

        else
            @info "Drone $(available_drones[i]) goes unassigned this round"
        end # assgn != 0
    end # for enumerate(assignment)

    # @infiltrate
end # end function

function earliest_due_date!(server::RoutingAllocation, routing_sim::RoutingSimulator,
                            rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}

    # First collect all available drones
    available_drones = Dict{Int,Vector{String}}()

    # Split out available drones in terms of depot
    # Since drones in same depot will use priority ordering
    for (drone_nm, dp) in server.agent_prop_set

        drone = server.agent_set[drone_nm]
        if dp.at_depot == true
            if ~(haskey(available_drones, drone.depot_number))
                available_drones[drone.depot_number] = [drone_nm]
            else
                push!(available_drones[drone.depot_number], drone_nm)
            end
        end
    end # for (drone_nm, dp)

    all_assigned_pkgs = Set{String}()

    # Assign all available drones
    for (depot_number, depot_drones) in available_drones

        depot_loc = server.agent_set[depot_drones[1]].depot_loc

        pkgs_in_range = Set{String}()
        for (pkg_nm, pp) in routing_sim.active_packages
            if Distances.evaluate(EuclideanLatLongMetric(), convert_to_vector(depot_loc), convert_to_vector(pp.delivery)) <= routing_sim.distance_thresh
                push!(pkgs_in_range, pkg_nm)
            end
        end # for (pkg_nm, pp)

        for drone_nm in depot_drones

            # Iterate over interaction events and add first one
            # that is not already assigned
            ie_idx = 0
            for (idx, ie) in enumerate(server.agent_prop_set[drone_nm].interaction_events)
                if ~(ie.task_name in all_assigned_pkgs) && (ie.task_name in pkgs_in_range) &&
                    ~(haskey(routing_sim.busy_packages, ie.task_name))
                    ie_idx = idx
                end
            end

            # Now assign and sample true deliv-return times
            if ie_idx > 0
                pkg_nm = server.agent_prop_set[drone_nm].interaction_events[ie_idx].task_name

                if ~(haskey(routing_sim.busy_packages, pkg_nm))

                    push!(all_assigned_pkgs, pkg_nm)

                    server.agent_task_allocation[drone_nm] = (name=pkg_nm, time=Inf)

                    routing_sim.true_delivery_return[(drone_nm, pkg_nm)] = sample_true_delivery_return_time(server.agent_task_windows[(drone_nm, pkg_nm)],
                                                                                server.current_time,
                                                                                routing_sim.tt_est_std_scale,
                                                                                rng)
                    server.agent_prop_set[drone_nm].at_depot = false

                    new_busy_package = routing_sim.active_packages[pkg_nm]
                    routing_sim.busy_packages[pkg_nm] = new_busy_package
                    delete!(routing_sim.active_packages, pkg_nm)
                    routing_sim.num_active_packages -= 1
                end
            end # if ie_idx != -1
        end # for drone_nm in depot_drones
    end # for (dn, dd) in available_drones
end # end function

## Implementations for POMDPs.jl
@with_kw struct RoutingDroneState
    drone_idx::Int64
    pkg_assignment::Int64 # -1 at return, 0 at depot, and > 0 for on the way
    current_time::Int64
    true_delivery_return_tuple::Tuple{Float64,Float64}
end

"""
Only MCTS for routing
"""
@with_kw mutable struct RoutingMCTSMDP <: POMDPs.MDP{RoutingDroneState,Int64}
    routing_sim::RoutingSimulator
    server::RoutingAllocation
    horizon::Int64
    drone_pkg_assignment::Vector{Int64}
    true_delivery_return::Vector{Tuple{Float64,Float64}}
end

function update_routing_mcts_fullstate!(mdp::RoutingMCTSMDP, rng::AbstractRNG)

    drone_pkg_assignment = Int64[]
    true_delivery_return = Tuple{Float64,Float64}[]

    for (d, drone_nm) in enumerate(mdp.server.agent_ordering)

        drone_props = mdp.server.agent_prop_set[drone_nm]

        if ~(isempty(drone_props.current_package))
            p = parse(Int64, SubString(drone_props.current_package, 4))
            del_return = sample_true_delivery_return_time(mdp.server.agent_task_windows[(drone_nm, drone_props.current_package)],
                                                          mdp.server.current_time,
                                                          mdp.routing_sim.tt_est_std_scale,
                                                          rng)
            push!(true_delivery_return, del_return)
            push!(drone_pkg_assignment, p)
        else
            if drone_props.at_depot == true
                push!(drone_pkg_assignment, 0)
                push!(true_delivery_return, (0.0, 0.0))
            else
                # Drone is returning - should still have an assignment in agent_task_allocation
                # We SAMPLE that for the true del return time
                push!(drone_pkg_assignment, -1)

                del_return = nothing
                for (k, v) in mdp.server.agent_task_allocation
                    if k == drone_nm
                        pkg_nm, _ = v
                        del_return = sample_true_delivery_return_time(mdp.server.agent_task_windows[(drone_nm, pkg_nm)],
                                                                      mdp.server.current_time,
                                                                      mdp.routing_sim.tt_est_std_scale,
                                                                      rng)
                        break
                    end
                end

                # We MUST get a del_return by construction
                @assert ~(isnothing(del_return))

                push!(true_delivery_return, del_return)
            end
        end
    end

    mdp.drone_pkg_assignment = drone_pkg_assignment
    mdp.true_delivery_return = true_delivery_return
end

function get_current_routing_mcts_state(mdp::RoutingMCTSMDP, idx::Int64)
    return RoutingDroneState(idx, mdp.drone_pkg_assignment[idx], convert(Int64, mdp.server.current_time), mdp.true_delivery_return[idx])
end

POMDPs.discount(mdp::RoutingMCTSMDP) = 1.0
POMDPs.isterminal(mdp::RoutingMCTSMDP, s::RoutingDroneState) = (s.current_time == mdp.horizon)

function POMDPs.actions(mdp::RoutingMCTSMDP, s::RoutingDroneState)

    drone_nm = mdp.server.agent_ordering[s.drone_idx]
    in_range_pkgs = [0]

    if s.pkg_assignment == 0 # At depot

        pkgs_to_exclude = Set{Int64}()

        # Iterate over drones up to that point and exclude packages
        for i = 1:s.drone_idx-1
            if mdp.drone_pkg_assignment[i] > 0
                push!(pkgs_to_exclude, mdp.drone_pkg_assignment[i])
            end
        end

        depot_loc = mdp.server.agent_set[drone_nm].depot_loc
        for (pkg_nm, pp) in mdp.routing_sim.active_packages
            pkg_id = parse(Int64, SubString(pkg_nm, 4)) # Extracting number from name
            if Distances.evaluate(EuclideanLatLongMetric(), convert_to_vector(depot_loc), convert_to_vector(pp.delivery)) <= mdp.routing_sim.distance_thresh &&
                ~(pkg_id in pkgs_to_exclude)
                push!(in_range_pkgs, pkg_id)
            end
        end

    end

    return in_range_pkgs
end

function POMDPs.gen(mdp::RoutingMCTSMDP, s::RoutingDroneState, a::Int64, rng::AbstractRNG)

    new_pkg_assignment = s.pkg_assignment
    new_true_del_return_tuple = s.true_delivery_return_tuple
    next_time = s.current_time + 1
    on_time_pkg = 0

     if s.pkg_assignment == 0

         if a > 0
            drone_nm = mdp.server.agent_ordering[s.drone_idx]
            pkg_nm = string("pkg", a)
            new_true_del_return_tuple = sample_true_delivery_return_time(mdp.server.agent_task_windows[(drone_nm, pkg_nm)],
                                                                          mdp.server.current_time,
                                                                          mdp.routing_sim.tt_est_std_scale,
                                                                          rng)
            new_pkg_assignment = a
        end
    else
        @assert s.true_delivery_return_tuple != (0.0, 0.0)
        del_time, ret_time = s.true_delivery_return_tuple

        if s.pkg_assignment == -1
            if ret_time > s.current_time && ret_time < next_time # Back to base; otherwise change nothing
                new_pkg_assignment = 0
                new_true_del_return_tuple = (0.0, 0.0)
            end
        else
            # On the way
            pkg_nm = string("pkg", s.pkg_assignment)
            if del_time > s.current_time && del_time < next_time # Delivered in sim

                if haskey(mdp.routing_sim.busy_packages, pkg_nm)
                    tw = mdp.routing_sim.busy_packages[pkg_nm].time_window[2]
                elseif haskey(mdp.routing_sim.active_packages, pkg_nm)
                    tw = mdp.routing_sim.active_packages[pkg_nm].time_window[2]
                else
                    throw(ErrorException("Neither busy nor active has pkg!"))
                end

                if del_time <= tw
                    on_time_pkg = 1
                end
                new_pkg_assignment = -1 # On the way back
            end
        end # s.pkg_assignment == -1
    end # s.pkg_assignment == 0
    return (sp=RoutingDroneState(s.drone_idx, new_pkg_assignment, next_time, new_true_del_return_tuple), r = on_time_pkg)
end # end function

function update_server_with_mdp_action!(mdp::RoutingMCTSMDP, action::Vector{Int64}, rng::AbstractRNG)

    for (agent_action, drone_nm) in zip(action, mdp.server.agent_ordering)
        if agent_action > 0 # Do nothing otherwise
            pkg_nm = string("pkg", agent_action)
            mdp.server.agent_task_allocation[drone_nm] = (name=pkg_nm, time=Inf)

            mdp.routing_sim.true_delivery_return[(drone_nm, pkg_nm)] =
                        sample_true_delivery_return_time(mdp.server.agent_task_windows[(drone_nm, pkg_nm)],
                                                         mdp.server.current_time,
                                                         mdp.routing_sim.tt_est_std_scale,
                                                         rng)
            mdp.server.agent_prop_set[drone_nm].at_depot = false
            mdp.server.agent_prop_set[drone_nm].current_package = pkg_nm

            new_busy_package = mdp.routing_sim.active_packages[pkg_nm]
            mdp.routing_sim.busy_packages[pkg_nm] = new_busy_package
            delete!(mdp.routing_sim.active_packages, pkg_nm)
            mdp.routing_sim.num_active_packages -= 1
        end
    end
end

struct NearestPkgPolicy <: Policy
    mdp::RoutingMCTSMDP
end

function POMDPs.action(pol::NearestPkgPolicy, s::RoutingDroneState)

    mdp = pol.mdp
    action = 0

    if s.pkg_assignment == 0 # At depot, assign to nearest available
        pkgs_to_exclude = Set{Int64}()

        # Iterate over drones up to that point and exclude packages
        for i = 1:s.drone_idx-1
            if mdp.drone_pkg_assignment[i] > 0
                push!(pkgs_to_exclude, mdp.drone_pkg_assignment[i])
            end
        end

        drone_nm = mdp.server.agent_ordering[s.drone_idx]
        depot_loc = mdp.server.agent_set[drone_nm].depot_loc
        min_pkg = 0
        min_dist = mdp.routing_sim.distance_thresh

        for (pkg_nm, pp) in mdp.routing_sim.active_packages
            pkg_id = parse(Int64, SubString(pkg_nm, 4)) # Extracting number from name
            dist = Distances.evaluate(EuclideanLatLongMetric(), convert_to_vector(depot_loc), convert_to_vector(pp.delivery))
            if dist <= min_dist && ~(pkg_id in pkgs_to_exclude)
                min_dist = dist
                min_pkg = pkg_id
            end
        end

        if min_pkg != 0
            action = min_pkg
        end
    end

    return action
end


text:
smooth_up_filter = r_type.( reshape((1/4) * [1 2 1;2 4.0 2;1 2 1],3,3,1,1))
smooth_down_filter =r_type.( reshape((1/16) * [1 2 1;2 4 2;1 2 1],3,3,1,1))
laplacian_filter = r_type.(reshape([0 -1 0;-1 4.0 -1;0 -1 0],3,3,1,1))

function block_filter!(filter_size, kernel, channels)
    w = zeros(r_type, filter_size, filter_size, channels, channels)
    for i in 1:channels
        w[:,:,i,i] = kernel
    end
    return w
end

function big_block_filter!(filter_size, kernel, channels)
    w = u_type.(zeros(r_type, filter_size, filter_size, channels, channels))|>cgpu
    for i in 1:channels
        w[:,:,:,i] = u_type.(kernel)|>cgpu
    end
    return w
end

block_laplacian_filter = block_filter!(3, laplacian_filter, 2)

up = ConvTranspose(smooth_up_filter, r_type.([0.0]), stride=2)|> pu
down = Conv(smooth_down_filter, r_type.([0.0]), stride=2)|> pu

block_up = ConvTranspose(block_filter!(3, smooth_up_filter, 2), r_type.([0.0]), stride=2)
block_down = Conv(block_filter!(3, smooth_down_filter, 2), r_type.([0.0]), stride=2)

i_conv = Conv(block_filter!(1, reshape([1.0],1,1,1,1), 2),r_type.([0.0]))|> pu

function laplacian_conv!(grid; h= 1)
    filter = r_type.((1.0 / (h^2)) * laplacian_filter)
    conv = Conv(filter, r_type.([0.0]), pad=(1,1))
    return conv(grid)
end

function helmholtz_chain!(grid, matrix; h = 1)
    filter = r_type.((1.0 / (h^2)) * laplacian_filter)
    conv = Conv(filter, r_type.([0.0]), pad=(1,1))|>pu
    y = conv(grid|>pu) - matrix .* grid
    return y
end

function helmholtz_chain_channels!(grid, matrix; h = 1)
    filter = r_type.((1.0 / (h^2)) * block_laplacian_filter)
    conv = Conv(filter, [0.0], pad=(1,1))|> pu
    y = conv(grid)-sum(matrix .* grid, dims=4)
    return y
end


text:
function TAIUTCDatFilesLoader(arg0::JString)
    return TAIUTCDatFilesLoader((JString,), arg0)
end

function TAIUTCDatFilesLoader(arg0::JString, arg1::DataProvidersManager)
    return TAIUTCDatFilesLoader((JString, DataProvidersManager), arg0, arg1)
end

function load_offsets(obj::TAIUTCDatFilesLoader)
    return jcall(obj, "loadOffsets", List, ())
end


text:
"""
Google Cloud Datastore API
"""
module _datastore

using Base64

export datastore

using ..api
using ...root

"""Enumerations"""
module types
    export ValueType, OperatorType, wrap, unwrap

    using Dates

    import JSON

    """Datastore value types"""
    @enum(ValueType,
        nullValue,
        booleanValue, integerValue, doubleValue, timestampValue,
        keyValue, stringValue, blobValue,
        geoPointValue, entityValue, arrayValue
    )
    value_type_map = Dict(Symbol(x) => x for x in instances(ValueType))
    Base.convert(::Type{ValueType}, x::Symbol) = haskey(value_type_map, x) ? value_type_map[x] : error("Unknown Datastore value type: $x")
    JSON.lower(x::ValueType) = string(x)
    for s in instances(ValueType)
        @eval export $(Symbol(s))
    end

    """Datastore comparison operator types"""
    @enum(OperatorType, LESS_THAN, LESS_THAN_OR_EQUAL, GREATER_THAN, GREATER_THAN_OR_EQUAL, EQUAL, HAS_ANCESTOR)
    operator_type_map = Dict(Symbol(x) => x for x in instances(OperatorType))
    Base.convert(::Type{OperatorType}, x::Symbol) = haskey(operator_type_map, x) ? operator_type_map[x] : error("Unknown Datastore comparison operator type: $x")
    JSON.lower(x::OperatorType) = string(x)
    for s in instances(OperatorType)
        @eval export $(Symbol(s))
    end

    value_types = Dict{Type, ValueType}(
        Bool => booleanValue,
        Integer => integerValue,
        AbstractFloat => doubleValue,
        TimeType => timestampValue,
        AbstractString => stringValue,
        Char => stringValue,
        Enum => stringValue,
        Nothing => nullValue,
    )
    function wrap(x)
        T = typeof(x)
        if T <: AbstractArray
            return Dict(arrayValue => Dict(:values => map(wrap, x)))
        end
        # move up type hierarchy until we get a hit
        while T !== Any
            if haskey(value_types, T)
                return Dict(value_types[T] => x)
            end
            T = supertype(T)
        end
        Dict(blobValue => base64encode(JSON.json(x)))
    end
    function unwrap(x)
        value_type, value = first(x)
        value_type = convert(ValueType, value_type)
        if value_type == timestampValue
            DateTime(chop(value), ISODateTimeFormat)
        elseif value_type == integerValue
            parse(Int64, value)
        elseif value_type == arrayValue
            map(unwrap, value[:values])
        elseif value_type == blobValue
            JSON.parse(String(base64decode(value)))
        else
            value
        end
    end
end

"""
Google Cloud Datastore API root.
"""
datastore = APIRoot(
    "https://datastore.googleapis.com/v1/projects/{project}",
    Dict(
        "cloud-platform" => "Full access to all resources and services in the specified Cloud Platform project.",
        "datastore" => "View and manage your Google Cloud Datastore data",
    );
    Project=APIResource("";
        allocateIds=APIMethod(:DELETE, ":allocateIds", "Allocates IDs for the given keys, which is useful for referencing an entity before it is inserted."),
        beginTransaction=APIMethod(:POST, ":beginTransaction", "Begins a new transaction."),
        commit=APIMethod(:POST, ":commit", "Commits a transaction, optionally creating, deleting or modifying some entities."),
        lookup=APIMethod(:POST, ":lookup", "Looks up entities by key."),
        rollback=APIMethod(:POST, ":rollback", "Rolls back a transaction."),
        runQuery=APIMethod(:POST, ":runQuery", "Queries for entities."),
    ),
)

end


text:
# Pizza Slices
# Let's try this out on a problem. Let's say you want to order pizzas for 10 people (people = 10) and each person wants 2.5 slices on average (avg = 2.5). A pizza has 8 slices per pizza (slices = 8). How many pizzas should you order (pizzas = ?)? So we have the following

# Meaning	Variable
# Number of people	people
# Average number of slices each person eats	avg
# Number of slices on a piece of pizza	slices

people = 10;
avg = 2.5;
slices = 8;

pizzas = people * avg / slices;

println(pizzas);

pizzasRounded = ceil(pizzas);

println(pizzasRounded);

max_stars_repo_name: fahimfarhan/legendary-coding-odyssey

max_issues_repo_name: fahimfarhan/legendary-coding-odyssey

max_forks_repo_name: fahimfarhan/legendary-coding-odyssey


text:
using Dispersal, Aqua, SafeTestsets

if VERSION >= v"1.5.0"
    # Amibiguities are not owned by Dispersal
    # Aqua.test_ambiguities([Dispersal, Base, Core])
    Aqua.test_unbound_args(Dispersal)
    Aqua.test_undefined_exports(Dispersal)
    Aqua.test_project_extras(Dispersal)
    Aqua.test_stale_deps(Dispersal)
    Aqua.test_deps_compat(Dispersal)
    Aqua.test_project_toml_formatting(Dispersal)
end

@time @safetestset "integration" begin include("integration.jl") end
@time @safetestset "downsampling" begin include("downsampling.jl") end
@time @safetestset "human dispersal" begin include("human.jl") end
@time @safetestset "jump dispersal" begin include("jump.jl") end
@time @safetestset "kernels" begin include("kernels.jl") end
@time @safetestset "allee effects" begin include("allee.jl") end
@time @safetestset "growth" begin include("growth.jl") end
@time @safetestset "mortality" begin include("mortality.jl") end
@time @safetestset "chain" begin include("chain.jl") end


text:
const BAR = (
            " ",
            LOWER_ONE_EIGHTH_BLOCK,
            LOWER_ONE_QUARTER_BLOCK,
            LOWER_THREE_EIGHTHS_BLOCK,
            LOWER_HALF_BLOCK,
            LOWER_FIVE_EIGHTHS_BLOCK,
            LOWER_THREE_QUARTERS_BLOCK,
            LOWER_SEVEN_EIGHTHS_BLOCK,
            FULL_BLOCK,
            )

Base.@kwdef struct BarChart
    block::Block
    data::Vector{Tuple{String, Int}}
    width::Int = 0
    gap::Int = 0
end

function draw(bar_chart::BarChart, area::Rect, buf::Buffer)

    draw(bar_chart.block, area, buf)
    chart_area = inner(bar_chart.block, area)

    w = bar_chart.width
    gap = bar_chart.gap
    if w == 0 || gap == 0
        w = ( width(chart_area) ÷ length(bar_chart.data) ) - 1
        gap = 1
    end

    height(chart_area) < 2 && return

    max_data = maximum([v for (_, v) in bar_chart.data])

    data = [Int(round(d * (height(chart_area) - 1) * 8 / max_data, RoundUp)) for (_, d) in bar_chart.data]

    for j in reverse(1:height(chart_area))
        for (i, d) in enumerate(data)

            symbol = d > 8 ? BAR[9] : BAR[Int(round(d / 8, RoundUp)) + 1]

            for x in 1:w

                set(
                    buf,
                    left(chart_area) + (i - 1) * (w + gap) + x,
                    top(chart_area) + j,
                    symbol,
                )

            end

            if d > 8
                d -= 8
            else
                d = 0
            end

            data[i] = d
        end
    end

    data = [Int(round(d * (height(chart_area) - 1) * 8 / max_data, RoundUp)) for (_, d) in bar_chart.data]

    for (i, d) in enumerate(data)
        s = bar_chart.data[i][1]
        if length(s) < w
            set(
                buf,
                left(chart_area) + (i - 1) * (w + gap) + (w ÷ 2) - (length(s) ÷ 2) + 1,
                top(chart_area) + height(chart_area) - d ÷ 8 - 1,
                s,
            )
        end

    end

    # for (i, &(label, value)) in self.data.iter().take(max_index).enumerate() {
    #     if value != 0 {
    #         let value_label = &self.values[i];
    #         let width = value_label.width() as u16;
    #         if width < self.bar_width {
    #             buf.set_string(
    #                 chart_area.left()
    #                     + i as u16 * (self.bar_width + self.bar_gap)
    #                     + (self.bar_width - width) / 2,
    #                 chart_area.bottom() - 2,
    #                 value_label,
    #                 self.value_style,
    #             );
    #         }
    #     }
    #     buf.set_stringn(
    #         chart_area.left() + i as u16 * (self.bar_width + self.bar_gap),
    #         chart_area.bottom() - 1,
    #         label,
    #         self.bar_width as usize,
    #         self.label_style,
    #     );
    # }

end


text:
function main()
    #load txt file
    values = readlines("input.txt")
    #parse all values to int
    values = parse.(Int, values)
    #detect increases and count them
    inc_count = 0
    noninc_count = 0
    curr = -Inf
    for v in values
        println(string(v)*" > "*string(curr)*"?")
        if curr == -Inf
            #first entry of the list, don't count increases
            curr = v
            #println("First")
        elseif v > curr
            inc_count += 1
            curr = v
            #println("Yes!")
        else
            noninc_count += 1
            curr = v
            #println("No...")
        end
    end
    println(inc_count)
    println(noninc_count)
    println(size(values, 1))
    println(1 + inc_count + noninc_count)
end

main()


text:
using Test

using Petri
using AlgebraicPetri
using Catlab.Theories
using Catlab.CategoricalAlgebra.ShapeDiagrams
using Catlab.CategoricalAlgebra.FinSets

function compare_petricospan(p₁::PetriCospan, p₂::PetriCospan)
    @test dom(p₁) == dom(p₂)
    @test codom(p₁) == codom(p₂)
    @test force(left(p₁)) == force(left(p₂))
    @test force(right(p₁)) == force(right(p₂))
    @test decoration(p₁) == decoration(p₂)
end

@testset "Core" begin
    include("core.jl")
end

@testset "Petri" begin
    include("petri.jl")
end


text:
include("abstract_trajectory.jl")
include("trajectory.jl")


text:
using Test, StatsModels
using DataFrames

using Unfold
include("test_utilities.jl")

data, evts = loadtestdata("test_case_3a") #
##

f_spl = @formula 0 ~ 1 + conditionA + spl(continuousA, 4) # 1
f = @formula 0 ~ 1 + conditionA + continuousA # 1
data_r = reshape(data, (1, :))
data_e, times = Unfold.epoch(data = data_r, tbl = evts, τ = (-1.0, 1.0), sfreq = 10) # cut the data into epochs

m_mul = coeftable(fit(UnfoldModel, f, evts, data_e, times))
m_mul_spl = coeftable(fit(UnfoldModel, f_spl, evts, data_e, times))

# asking for 4 splines should generate 4 splines
@test length(unique(m_mul_spl.coefname)) == 6 # XXX check back with Unfold whether this is the same! could be n-1 splines in Unfold. We should keep that comparable I guess

# test safe prediction
m = fit(UnfoldModel, f_spl, evts, data_e, times)
r = predict(m,DataFrame(conditionA=[0,0],continuousA=[0.9,1.9]))
@test all(ismissing.(r.yhat[r.continuousA.==1.9]))
@test !any(ismissing.(r.yhat[r.continuousA.==0.9]))

# test time expanded
basisfunction = firbasis(τ = (-1, 1), sfreq = 10, name = "A")
m_tul = coeftable(fit(UnfoldModel, f, evts, data_r, basisfunction))
m_tul_spl = coeftable(fit(UnfoldModel, f_spl, evts, data_r, basisfunction))

# test safe predict
m = fit(UnfoldModel, f_spl, evts, data_r, basisfunction)
@test_broken predict(m,DataFrame(conditionA=[0,0,0],continuousA=[0.9,0.9,1.9]))
#@test_broken all(ismissing.)

#evts_grid = gridexpand()
# results from timeexpanded and non should be equal
#yhat_tul  = predict(m_tul_spl,evts_grid)
#yhat_mul  = predict(m_mul_spl,evts_grid)
if 1 == 0
    using AlgebraOfGraphics
    yhat_mul.conditionA = categorical(yhat_mul.conditionA)
    yhat_mul.continuousA = categorical(yhat_mul.continuousA)
    m = mapping(:times, :yhat, color = :continuousA, linestyle = :conditionA)
    df = yhat_mul
    AlgebraOfGraphics.data(df) * visual(Lines) * m |> draw
end

# test much higher number of splines
f_spl_many = @formula 0 ~ 1 + spl(continuousA, 131) # 1
m_mul_spl_many = coeftable(fit(UnfoldModel, f_spl_many, evts, data_e, times))
@test length(unique(m_mul_spl_many.coefname)) == 132


text:
"""
    BuoyancyTracer <: AbstractBuoyancyModel{Nothing}

Type indicating that the tracer `b` represents buoyancy.
"""
struct BuoyancyTracer <: AbstractBuoyancyModel{Nothing} end

const BuoyancyTracerModel = Buoyancy{<:BuoyancyTracer}

required_tracers(::BuoyancyTracer) = (:b,)

@inline buoyancy_perturbation(i, j, k, grid, ::BuoyancyTracer, C) = @inbounds C.b[i, j, k]

@inline ∂x_b(i, j, k, grid, ::BuoyancyTracer, C) = ∂xᶠᵃᵃ(i, j, k, grid, C.b)
@inline ∂y_b(i, j, k, grid, ::BuoyancyTracer, C) = ∂yᵃᶠᵃ(i, j, k, grid, C.b)
@inline ∂z_b(i, j, k, grid, ::BuoyancyTracer, C) = ∂zᵃᵃᶠ(i, j, k, grid, C.b)

@inline    top_buoyancy_flux(i, j, grid, ::BuoyancyTracer, top_tracer_bcs, clock, fields) = getbc(top_tracer_bcs.b, i, j, grid, clock, fields)
@inline bottom_buoyancy_flux(i, j, grid, ::BuoyancyTracer, bottom_tracer_bcs, clock, fields) = getbc(bottom_tracer_bcs.b, i, j, grid, clock, fields)


text:
using Fleck
using Documenter
using Literate

# Set Literate.jl config if not being compiled on recognized service.
config = Dict{String,String}()
if !(haskey(ENV, "GITHUB_ACTIONS") || haskey(ENV, "GITLAB_CI"))
  config["nbviewer_root_url"] = "https://nbviewer.jupyter.org/github/AlgebraicJulia/AlgebraicPetri.jl/blob/gh-pages/dev"
  config["repo_root_url"] = "https://github.com/AlgebraicJulia/AlgebraicPetri.jl/blob/master/docs"
end

const literate_dir = joinpath(@__DIR__, "..", "examples")
const generated_dir = joinpath(@__DIR__, "src", "examples")

for (root, dirs, files) in walkdir(literate_dir)
  out_dir = joinpath(generated_dir, relpath(root, literate_dir))
  for file in files
    f,l = splitext(file)
    if l == ".jl" && !startswith(f, "_")
      Literate.markdown(joinpath(root, file), out_dir;
        config=config, documenter=true, credit=false)
    end
  end
end

makedocs(;
    modules=[Fleck],
    authors="Andrew Dolgert <adolgert@uw.edu>",
    repo="https://github.com/adolgert/Fleck.jl/blob/{commit}{path}#L{line}",
    sitename="Fleck.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://adolgert.github.io/Fleck.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Develop" => "develop.md"
    ],
)

deploydocs(;
    repo="github.com/adolgert/Fleck.jl",
)


text:
# add specific loop closure using Odo constraint
# yes, needs to be optimized and refactored

using IncrementalInference, CloudGraphs

## Uncomment out for command line operation
cloudGraph, addrdict = standardcloudgraphsetup()
session = addrdict["session"]

# interactive operation
# session = "SESSROX"
# include(joinpath(dirname(@__FILE__),"blandauthremote.jl"))

# Register protobuf types of interest in CloudGraphs layer
registerGeneralVariableTypes!(cloudGraph)

fg = emptyFactorGraph()
fg.cg = cloudGraph

fullLocalGraphCopy!(fg)

v1a = getVert(fg, :l9)
v1b = getVert(fg, :l16)

f1 = addFactor!(fg,[v1a;v1b],Odo(zeros(2,1),0.1*eye(2),[1.0]),ready=0)

setDBAllReady!(conn)

# x144, x109 and x15, 161 and x115, x137 and x122, 132 and x131, 123
# x37, x149 and x113, 139 and x112,141 and x17,x160, and x34,x154
# x119, x135 and x10, x175 and x38, x148 and x8, x181 and x21, x157
v2a = getVert(fg, :x21)
v2b = getVert(fg, :x157)

DX = zeros(3,1)
DX[3] = pi
pp = Pose2Pose2(DX, 0.1*eye(3), [1.0]) #[prev;v],
f = addFactor!(fg, [v2a;v2b], pp, ready=0)

setDBAllReady!(conn)

# DX = zeros(3)
# x164, x170

# 9/14 l4, l7 and l5, l6 and l3, l11 and l3, l13 and

# 8/22 l1, l6 and l1 and l10


text: abstract type UDWT_Filter_Biorthogonal{T<:Number} end


text:
module ViscosityDrop

using Oceananigans
using Optim
using JLD2
using DelimitedFiles

using Oceananigans.Grids
include("simulation/grid.jl")
export make_grid

include("simulation/boundary_conditions.jl")
export no_slip_bc

using Oceananigans.Advection
include("simulation/model.jl")
export make_model

include("simulation/buoyancy.jl")
export set_buoyancy!

using Oceananigans.AbstractOperations
using Oceananigans.Fields
using Oceananigans.OutputWriters
include("simulation/simulation.jl")
export make_simulation, run!

include("simulation/wrapper.jl")
export do_all

include("post_process/circle.jl")
include("post_process/fit_functions.jl")
include("post_process/get_results.jl")
export get_results

end # module


text:
#
# Turbulent kinetic energy wall models
#

# Fallbacks
TKEBoundaryConditions(T, wall_model) = DefaultBoundaryConditions(T)

#
# Prescribes the "near wall" turbulent kinetic energy
# at the upper grid cell at i = N.
#

Base.@kwdef struct PrescribedNearWallTKE{T} <: AbstractParameters
    Cʷu★ :: T = 3.75
end

@inline update_near_wall_tke!(m::Model{L, H, <:PrescribedNearWallTKE}) where {L, H} =
    m.solution.e[m.grid.N] = m.tke_wall_model.Cʷu★ * u★(m)^2

#
# Prescribes the surface value of turbulent kinetic energy
# via a ValueBoundaryCondition.
#

Base.@kwdef struct PrescribedSurfaceTKEValue{T} <: AbstractParameters
    Cʷu★ :: T = 3.75
end

@inline (boundary_tke::PrescribedSurfaceTKEValue)(model) = boundary_tke.Cʷu★ * u★(model)^2

TKEValueSurfaceConditions(T, wall_model::PrescribedSurfaceTKEValue) =
    FieldSurfaceConditions(GradientSurfaceCondition(-zero(T)), ValueSurfaceCondition(wall_model))

#
# Prescribes the surface flux of turbulent kinetic energy
# as propotional to the sum of u★^3 and wΔ^3 = Δz * max(0, Qb), where Qb is the buoyancy flux.
# The free parameters for both u★ and wΔ are multiplied by the TKE dissipation parameter
# to prevent their mutual correlation.
#

Base.@kwdef struct PrescribedSurfaceTKEFlux{T} <: AbstractParameters
    Cʷu★ :: T = 1.3717
    CʷwΔ :: T = 1.0
end

@inline prescribed_surface_tke_flux(model) =  - model.tke_equation.Cᴰ * (   model.tke_wall_model.Cʷu★ * u★(model)^3
                                                                          + model.tke_wall_model.CʷwΔ * wΔ³(model) )

TKEBoundaryConditions(T, wall_model::PrescribedSurfaceTKEFlux) =
    FieldBoundaryConditions(GradientBoundaryCondition(-zero(T)), FluxBoundaryCondition(prescribed_surface_tke_flux))


text:
"""
    basicblock(inplanes, outplanes, downsample = false)

Create a basic residual block
([reference](https://arxiv.org/abs/1512.03385v1)).

# Arguments:
- `inplanes`: the number of input feature maps
- `outplanes`: a list of the number of output feature maps for each convolution
               within the residual block
- `downsample`: set to `true` to downsample the input
"""
function basicblock(inplanes, outplanes, downsample = false)
  stride = downsample ? 2 : 1
  Chain(conv_bn((3, 3), inplanes, outplanes[1]; stride = stride, pad = 1, bias = false),
        conv_bn((3, 3), outplanes[1], outplanes[2], identity; stride = 1, pad = 1, bias = false))
end

"""
    bottleneck(inplanes, outplanes, downsample = false; stride = [1, (downsample ? 2 : 1), 1])

Create a bottleneck residual block
([reference](https://arxiv.org/abs/1512.03385v1)). The bottleneck is composed of
3 convolutional layers each with the given `stride`.
By default, `stride` implements ["ResNet v1.5"](https://catalog.ngc.nvidia.com/orgs/nvidia/resources/resnet_50_v1_5_for_pytorch)
which uses `stride == [1, 2, 1]` when `downsample == true`.
This version is standard across various ML frameworks.
The original paper uses `stride == [2, 1, 1]` when `downsample == true` instead.

# Arguments:
- `inplanes`: the number of input feature maps
- `outplanes`: a list of the number of output feature maps for each convolution
               within the residual block
- `downsample`: set to `true` to downsample the input
- `stride`: a list of the stride of the 3 convolutional layers
"""
function bottleneck(inplanes, outplanes, downsample = false;
                    stride = [1, (downsample ? 2 : 1), 1])
  Chain(conv_bn((1, 1), inplanes, outplanes[1]; stride = stride[1], bias = false),
        conv_bn((3, 3), outplanes[1], outplanes[2]; stride = stride[2], pad = 1, bias = false),
        conv_bn((1, 1), outplanes[2], outplanes[3], identity; stride = stride[3], bias = false))
end

"""
    bottleneck_v1(inplanes, outplanes, downsample = false)

Create a bottleneck residual block
([reference](https://arxiv.org/abs/1512.03385v1)). The bottleneck is composed of
3 convolutional layers with all a stride of 1 except the first convolutional
layer which has a stride of 2.

# Arguments:
- `inplanes`: the number of input feature maps
- `outplanes`: a list of the number of output feature maps for each convolution
               within the residual block
- `downsample`: set to `true` to downsample the input
"""
bottleneck_v1(inplanes, outplanes, downsample = false) =
    bottleneck(inplanes, outplanes, downsample; stride = [(downsample ? 2 : 1), 1, 1])

"""
    resnet(block, residuals::NTuple{2, Any}, connection = addrelu;
           channel_config, block_config, nclasses = 1000)

Create a ResNet model
([reference](https://arxiv.org/abs/1512.03385v1)).

# Arguments
- `block`: a function with input `(inplanes, outplanes, downsample=false)` that returns
           a new residual block (see [`Metalhead.basicblock`](#) and [`Metalhead.bottleneck`](#))
- `residuals`: a 2-tuple of functions with input `(inplanes, outplanes, downsample=false)`,
               each of which will return a function that will be used as a new "skip" path to match a residual block.
              [`Metalhead.skip_identity`](#) and [`Metalhead.skip_projection`](#) can be used here.
- `connection`: the binary function applied to the output of residual and skip paths in a block
- `channel_config`: the growth rate of the output feature maps within a residual block
- `block_config`: a list of the number of residual blocks at each stage
- `nclasses`: the number of output classes
"""
function resnet(block, residuals::AbstractVector{<:NTuple{2, Any}}, connection = addrelu;
                channel_config, block_config, nclasses = 1000)
  inplanes = 64
  baseplanes = 64
  layers = []
  push!(layers, conv_bn((7, 7), 3, inplanes; stride = 2, pad = 3, bias = false))
  push!(layers, MaxPool((3, 3), stride = (2, 2), pad = (1, 1)))
  for (i, nrepeats) in enumerate(block_config)
    # output planes within a block
    outplanes = baseplanes .* channel_config
    # push first skip connection on using first residual
    # downsample the residual path if this is the first repetition of a block
    push!(layers, Parallel(connection, block(inplanes, outplanes, i != 1),
                                       residuals[i][1](inplanes, outplanes[end], i != 1)))
    # push remaining skip connections on using second residual
    inplanes = outplanes[end]
    for _ in 2:nrepeats
      push!(layers, Parallel(connection, block(inplanes, outplanes, false),
                                         residuals[i][2](inplanes, outplanes[end], false)))
      inplanes = outplanes[end]
    end
    # next set of output plane base is doubled
    baseplanes *= 2
  end

  return Chain(Chain(layers),
               Chain(AdaptiveMeanPool((1, 1)), MLUtils.flatten, Dense(inplanes, nclasses)))
end

"""
    resnet(block, shortcut_config::Symbol, connection = addrelu;
           channel_config, block_config, nclasses = 1000)

Create a ResNet model
([reference](https://arxiv.org/abs/1512.03385v1)).

# Arguments
- `block`: a function with input `(inplanes, outplanes, downsample=false)` that returns
           a new residual block (see [`Metalhead.basicblock`](#) and [`Metalhead.bottleneck`](#))
- `shortcut_config`: the type of shortcut style (either `:A`, `:B`, or `:C`)
    - `:A`: uses a [`Metalhead.skip_identity`](#) for all residual blocks
    - `:B`: uses a [`Metalhead.skip_projection`](#) for the first residual block
            and [`Metalhead.skip_identity`](@) for the remaining residual blocks
    - `:C`: uses a [`Metalhead.skip_projection`](#) for all residual blocks
- `connection`: the binary function applied to the output of residual and skip paths in a block
- `channel_config`: the growth rate of the output feature maps within a residual block
- `block_config`: a list of the number of residual blocks at each stage
- `nclasses`: the number of output classes
"""
function resnet(block, shortcut_config::AbstractVector{<:Symbol}, args...; kwargs...)
  shortcut_dict = Dict(
    :A => (skip_identity, skip_identity),
    :B => (skip_projection, skip_identity),
    :C => (skip_projection, skip_projection))

  if any(sc -> !haskey(shortcut_dict,sc),shortcut_config)
    error("Unrecognized shortcut_config ($shortcut_config) passed to `resnet` (use only :A, :B, or :C).")
  end

  shortcut = [shortcut_dict[sc] for sc in shortcut_config]
  resnet(block, shortcut, args...; kwargs...)
end

function resnet(block, shortcut_config::Symbol, args...; block_config, kwargs...)
    resnet(block, fill(shortcut_config, length(block_config)), args...;
           block_config = block_config, kwargs...)
end

resnet(block, residuals::NTuple{2}, args...; kwargs...) = resnet(block, [residuals], args...; kwargs...)

const resnet_config =
  Dict(18 => (([1, 1], [2, 2, 2, 2], [:A, :B, :B, :B]), basicblock),
       34 => (([1, 1], [3, 4, 6, 3], [:A, :B, :B, :B]), basicblock),
       50 => (([1, 1, 4], [3, 4, 6, 3], [:B, :B, :B, :B]), bottleneck),
       101 => (([1, 1, 4], [3, 4, 23, 3], [:B, :B, :B, :B]), bottleneck),
       152 => (([1, 1, 4], [3, 8, 36, 3], [:B, :B, :B, :B]), bottleneck))

"""
    ResNet(channel_config, block_config, shortcut_config;
           block, connection = addrelu, nclasses = 1000)

Create a `ResNet` model
([reference](https://arxiv.org/abs/1512.03385v1)).
See also [`resnet`](#).

# Arguments
- `channel_config`: the growth rate of the output feature maps within a residual block
- `block_config`: a list of the number of residual blocks at each stage
- `shortcut_config`: the type of shortcut style (either `:A`, `:B`, or `:C`).
   `shortcut_config` can also be a vector of symbols if different shortcut styles are applied to
   different residual blocks.
- `block`: a function with input `(inplanes, outplanes, downsample=false)` that returns
           a new residual block (see [`Metalhead.basicblock`](#) and [`Metalhead.bottleneck`](#))
- `connection`: the binary function applied to the output of residual and skip paths in a block
- `nclasses`: the number of output classes
"""
struct ResNet
  layers
end

function ResNet(channel_config, block_config, shortcut_config;
                block, connection = addrelu, nclasses = 1000)
  layers = resnet(block,
                  shortcut_config,
                  connection;
                  channel_config = channel_config,
                  block_config = block_config,
                  nclasses = nclasses)

  ResNet(layers)
end

@functor ResNet

(m::ResNet)(x) = m.layers(x)

backbone(m::ResNet) = m.layers[1]
classifier(m::ResNet) = m.layers[2]

"""
    ResNet(depth = 50; pretrain = false, nclasses = 1000)

Create a ResNet model with a specified depth
([reference](https://arxiv.org/abs/1512.03385v1))
following [these modification](https://catalog.ngc.nvidia.com/orgs/nvidia/resources/resnet_50_v1_5_for_pytorch)
referred as ResNet v1.5.

See also [`Metalhead.resnet`](#).

# Arguments
- `depth`: depth of the ResNet model. Options include (18, 34, 50, 101, 152).
- `nclasses`: the number of output classes

!!! warning
    Only `ResNet(50)` currently supports pretrained weights.

For `ResNet(18)` and `ResNet(34)`, the parameter-free shortcut style (type `:A`)
is used in the first block and the three other blocks use type `:B` connection
(following the implementation in PyTorch). The published version of
`ResNet(18)` and `ResNet(34)` used type `:A` shortcuts for all four blocks. The
example below shows how to create a 18 or 34-layer `ResNet` using only type `:A`
shortcuts:

```julia
using Metalhead

resnet18 = ResNet([1, 1], [2, 2, 2, 2], :A; block = Metalhead.basicblock)

resnet34 = ResNet([1, 1], [3, 4, 6, 3], :A; block = Metalhead.basicblock)
```

The bottleneck of the orginal ResNet model has a stride of 2 on the first
convolutional layer when downsampling (instead of the second convolutional layers
as in ResNet v1.5). The architecture of the orignal ResNet model can be obtained
as shown below:

```julia
resnet50_v1 = ResNet([1, 1, 4], [3, 4, 6, 3], :B; block = Metalhead.bottleneck_v1)
```
"""
function ResNet(depth::Int = 50; pretrain = false, nclasses = 1000)
    @assert depth in keys(resnet_config) "`depth` must be one of $(sort(collect(keys(resnet_config))))"

    config, block = resnet_config[depth]
    model = ResNet(config...; block = block, nclasses = nclasses)
    pretrain && loadpretrain!(model, string("ResNet", depth))
    model
end

# Compat with Metalhead 0.6; remove in 0.7
@deprecate ResNet18(; kw...) ResNet(18; kw...)
@deprecate ResNet34(; kw...) ResNet(34; kw...)
@deprecate ResNet50(; kw...) ResNet(50; kw...)
@deprecate ResNet101(; kw...) ResNet(101; kw...)
@deprecate ResNet152(; kw...) ResNet(152; kw...)


text:
using Bijectors, Random

Random.seed!(123456)

include("transform.jl")


text:
"""
    CoefficientsQadicReversed(k::Int, z::T, q::Int, n::Int) where {T<:Integer}

Return an array ``[ v_1, v_2, \\ldots, v_n ]`` of length `n`
such that ``k*(q^n-1)/z`` has the form
``v_n q^0 + v_{n-1} q^1 + \\cdots + v_1 q^{n-1}``,
that is, the result contains the `q`-adic coefficients of ``k*(q^n-1)/z``
in reversed order.

We assume that ``0 \\leq k \\leq z``, ``q > 1``, and ``n > 0``.

The (perhaps large) number ``k*(q^n-1)/z`` need not be created
in the computations.
"""
function CoefficientsQadicReversed(k::Int, z::T, q::Int, n::Int) where {T<:Integer}
    v::Vector{Int} = zeros(Int, n)
    if k == z
      for i in 1:n
        v[i] = q-1
      end
    else
      r = k
      for i in 1:n
        d, r = divrem(r * q, z)
        v[i] = d
      end
    end

    return v
end

# The following function is needed only for the sake of efficiency
# comparisons.
function CoefficientsQadic(ke::T, q::Int, n::Int) where {T<:Integer}
    v::Vector{Int} = zeros(Int, n)
    for i in 1:n
        ke, r = divrem(ke, q)
        v[i] = r
    end

    return v
end

function islessorequal(mon1::Vector{Int}, mon2::Vector{Int}, n::Int)
  for i in 1:n
    if mon2[i] < mon1[i]
      return false
    end
  end

  return true
end

"""
    LoewyStructureInfo(q::Int, n::Int, z::T) where {T<:Integer}

Return a dictionary with the following keys.

`:monomials`
    the array of reversed `q`-adic expansions for multiples of
    `e = (q^n-1)/z`, each of length `n`,
    see [`CoefficientsQadicReversed`](@ref),

`:layers`
    the array of the Loewy layers to which the monomials belong,

`:chain`
    an array of positions of monomials of a longest ascending chain,

`:m`
    the value `m(q,e)`,

`:LL`
    the Loewy length of `A[q,n,z]`,
    equal to the length of the `:layers` value plus 1,

`:parameters`
    the array `[ q, n, z ]`.
"""
function LoewyStructureInfo(q::Int, n::Int, z::T) where {T<:Integer}
    if (z != 1 && powermod(q, n, z) != 1)
      error("<z> must divide <q>^<n> - 1")
    end

    monomials = [zeros(Int, n)]
    layers = Int[0]
    degrees = Int[0]    # just for speedup: avoid comparisons
    predecessors = Int[0]

    m = n * (q-1)
    for i in 1:z
      mon = CoefficientsQadicReversed(i, z, q, n)
      lambda = 0
      pred = 1
      mm = Base.sum(mon)
      for j in 2:i
        if lambda < layers[j] && degrees[j] < mm &&
                                 islessorequal(monomials[j], mon, n)
          lambda = layers[j]
          pred = j
        end
      end
      push!(monomials, mon)
      push!(layers, lambda + 1)
      push!(degrees, mm)
      push!(predecessors, pred)

      if lambda == 0
        if mm < m
          m = mm
        end
      end
    end

    # Extract information about one longest chain.
    i = length(monomials)
    pred = Int[]
    while i > 0
      push!(pred, i)
      i = predecessors[i]
    end

    return Dict(:monomials => monomials,
                :layers => layers,
                :chain => pred,
                :m => m,
                :LL => layers[z+1] + 1,
                :parameters => [q, n, z])
end;

"""
    LoewyVector(data::Dict)

Return an array that stores at position `i+1` the dimension of the `i`-th
Loewy layer of the Singer algebra described by the dictionary `data`.
"""
function LoewyVector(data::Dict) Vector{Int}
    v = zeros(Int, data[:LL])
    for i in data[:layers]
      v[i+1] = v[i+1] + 1
    end

    return v
end;

function SingerAlgE(q::Int, n::Int, z::T) where {T <: Integer}
    return div(BigInt(q)^n-1, z)
end;

"""
    LoewyLength(q::Int, z::T) where {T <: Integer}
    LoewyLength(q::Int, n::Int, z::T, m::Int) where {T <: Integer}
    LoewyLength(q::Int, n::Int, z::T, A::GAP.GapObj) where {T <: Integer}
    LoewyLength(q::Int, n::Int, z::T, m::Int, A::GAP.GapObj) where {T <: Integer}

Return the Loewy length of the Singer Algebra ``A[q, z]``.

If `n` and `m` are given then `n` must be the multiplicative order of `q`
modulo `z`, and `m` must be the value ``m(q, e)``.

If `A` is given then it must be the GAP object that represents ``A[q, z]``;
in this case, stored attributes of `A` are used.
"""
function LoewyLength(q::Int, z::T) where {T <: Integer}
    n = OrderMod(q, z)
    return LoewyLength(q, n, z,
               MinimalDegree(q, SingerAlgE(q, n, z)))
end;

function LoewyLength(q::Int, n::Int, z::Int, m::Int)
    l = LoewyLengthCheap(q, n, z, m)
    if l == 0
      l = LoewyLengthHard(q, n, z, m)
    end

    return l
end;

function LoewyLength(q::Int, n::Int, z::T, A::GAP.GapObj) where {T <: Integer}
    # First we need m; we store this attribute in `A`.
    if GAP.Globals.HasMinimalDegreeOfSingerAlgebraJulia(A)
      m = GAP.Globals.MinimalDegreeOfSingerAlgebraJulia(A)
    else
      m = MinimalDegree(q, n, SingerAlgE(q, n, z), A)
    end

    if GAP.Globals.HasLoewyStructureInfoJulia(A)
      # The Loewy info is now known.
      l = GAP.Globals.LoewyStructureInfoJulia(A)[:LL]
    else
      l = LoewyLength(q, n, z, m, A)
    end

    return l
end;

function LoewyLength(q::Int, n::Int, z::T, m::Int, A::GAP.GapObj) where {T <: Integer}
    # Compute the value.
    l = LoewyLengthCheap(q, n, z, m)
    if l == 0
      info = LoewyStructureInfo(q, n, z)
      # Store the attribute.
      GAP.Globals.SetLoewyStructureInfoJulia(A, info)
      l = info[:LL]
    end

    return l
end;

"""
    SufficientCriterionForLoewyBoundAttained(q::Int, n::Int, z::T, m::Int) where {T <: Integer}

Return a string that describes the criterion from
[BHHK1](https://doi.org/10.1016/j.jalgebra.2019.05.004) or
[BHHK2](http://arxiv.org/abs/1912.03065v1)
from which it follows that the Loewy length of the algebra ``A[q,z]``
is equal to the upper bound ``\\lfloor n (q-1) / m \\rfloor + 1``,
where ``n = \\ord_z(q)``, ``e = (q^n-1) / z``, and ``m = m(q, e)``;
if no such criterion was found then the returned string is empty.
"""
function SufficientCriterionForLoewyBoundAttained(q::Int, n::Int, z::T, m::Int) where {T <: Integer}

    if n <= 3
      return "Cor. I.7.1 (n <= 3)"
    elseif z < 70
      return "z < 70"
    elseif m == 2
      return "La. I.6.3"
    elseif mod(q - 1, m) == 0
      return "Thm. I.7.1"
    elseif n * (q-1) < 3 * m
      return "La. I.7.1 (iii)"
    end

    qm = mod(q, m)
    if 1 < qm && div(n * (qm-1), m) == 1
      return "Prop. II.3.15"
    end

    bq = BigInt(q)
    qnm1 = bq^n - 1
    e = div(qnm1, z)

    if e <= 32
      return "Prop. II.6.1 (e <= 32)"
    elseif e <= 3 * m
      return "Prop. II.6.4"
    elseif n <= 5 && mod(div(qnm1, q-1), e) == 0
      return "Prop. II.5.3, II.5.6 (e | (q^n-1)/(q-1), n <= 5)"
    elseif mod(n, div(m * OrderMod(q, e, n), gcd(m, q-1))) == 0
      return "La. II.5.2 (ii)"
    end

    if mod(div(qnm1, q-1), e) == 0
      # we know m <= n
      if m >= n-1
        return "Prop. II.5.1 (ii)"
      elseif m == n-2 && mod(q, m) in 1:div(m+1, 2)
        return "Prop. II.5.1 (iii)"
      elseif mod(div((n-m) * (q-1), m), n-m) == 0
        return "Prop. II.5.1 (i)"
      end
    end

    if mod(n, m) == 0
      if mod(Base.sum([bq^div(i*n, m) for i in 0:(m-1)]), e) == 0
        return "La. II.4.1 for r = 1"
      end
      for r in DivisorsInt(n)
        if mod(div(n, r), m) == 0 &&
           mod(Base.sum([bq^div(n*i, m*r) for i in 0:(m-1)]), e) == 0
          return "La. II.4.1"
        end
      end
    end

    p = SmallestRootInt(e)

    if p == 2
      return "Thm. II.4.3 (iii)"
    elseif PrimesGAP.isprime(p)
      # (Here we know that 'q mod p <> 1' holds.)
      # Check whether 'p' is a Pierpont prime.
      p = p - 1
      while iseven(p)
        p = div(p, 2)
      end
      while mod(p, 3) == 0
        p = div(p, 3)
      end
      if p == 1
        return "Thm. II.4.3 (ii)"
      end
    end

    # Concerning Remark II.5.4,
    # we need not check \Phi_2 (since m = e if e divides q-1,
    # and we have checked whether m divides q-1),
    # \Phi_4, \Phi_8 (because then m = 2).
    phi6 = bq^2 - q + 1     # this is \Phi_6(q)
    if mod(phi6, e) == 0 ||
       mod(bq^3*(bq^3+1) + 1, e) == 0 ||  # this is \Phi_9(q)
       mod(bq^3*(q-1) + phi6, e) == 0     # this is \Phi_{10}(q)
      return "Rem. II.5.4"
    end

    # We give up.
    return ""
end;

function LoewyLengthCheap(q::Int, n::Int, z::T, m::Int) where {T <: Integer}
    if SufficientCriterionForLoewyBoundAttained(q, n, z, m) != ""
      # The bound is attained.
      return div(n*(q-1), m) + 1
    else
      return 0
    end
end;

function LoewyLengthHard(q::Int, n::Int, z::T, m::Int) where {T <: Integer}
    return LoewyStructureInfo(q, n, z)[:LL]
end;

function OverflowTable(data::Dict)
    if ! haskey(data, :overflowTable)
      monoms = data[ :monomials ]
      d = length(monoms)
      mat = falses(d, d)
      q = data[:parameters][1]
      n = data[:parameters][2]
      for i in 1:d
        for j in 1:(i-1)
          res = all(x -> x < q, monoms[i] + monoms[j])
          mat[i, j] = res
          mat[j, i] = res
        end
        mat[i, i] = all(x -> x < q, monoms[i] + monoms[i])
      end

      data[:overflowTable] = mat
    end

    return data[:overflowTable]
end

function MultTable(data::Dict)
    if ! haskey(data, :multtable)
      monoms = data[:monomials]
      n = length(monoms)
      nn = div(n+1, 2)
      mat = zeros(Int, n, n)
      q = data[:parameters][1]
      if isodd(q)
        qhalf = div(q, 2) + 1
      else
        qhalf = div(q, 2)
      end
      for i in 1:nn
        for j in 1:(i-1)
          if all(x -> x < q, monoms[i] + monoms[j])
            mat[i, j] = i+j-1
            mat[j, i] = i+j-1
          end
        end
        if all(x -> x < qhalf, monoms[i])
          mat[i, i] = i+i-1
        end
      end
      for i in (nn+1):n
        # The [i,j] entry can be nonzero only if i+j <= z+2 holds.
        for j in 1:(n+1-i)
          if all(x -> x < q, monoms[i] + monoms[j])
            mat[i,j] = i+j-1
            mat[j,i] = i+j-1
          end
        end
      end
      data[:multtable] = mat
    end

    return data[:multtable]
end

function OverflowCount(data::Dict)
    if ! haskey(data, :overflowCount)
      mat = OverflowTable(data)
      list = [count(mat[i,:]) for i in 1:size(mat, 1)]
      data[:overflowCount] = list
    end

    return data[:overflowCount]
end

function BasisOfSum(data::Dict, list1::Vector{Int}, list2::Vector{Int})
    return sort(union(list1, list2))
end

function BasisOfIntersection(data::Dict, list1::Vector{Int}, list2::Vector{Int})
    return intersect(list1, list2)
end

function BasisOfProductSpace(data::Dict, list1::Vector{Int}, list2::Vector{Int})
    mat = MultTable(data)
    basis = Int[]

    for i in list1
      for j in list2
        if mat[i,j] != 0
          push!(basis, i+j-1)
        end
      end
    end

    return unique(sort(basis))
end

function BasisOfIdeal(data::Dict, list::Vector{Int})
    return BasisOfProductSpace(data, list, collect(1:length(data[:monomials])))
end

function BasisOfAnnihilator(data::Dict, list::Vector{Int})
    mat = MultTable(data)
    return filter(i -> all(j -> mat[i,j] == 0, list), 1:size(mat, 1))
end;

function BasesOfRadicalSeries(data::Dict)
    lay = data[:layers]
    ll = data[:LL]
    jbylayer = map(i -> findall(isequal(i), lay), 1:(ll-1))
    J = [jbylayer[ll-1]]
    for i in 2:(ll-1)
      push!(J, sort(union(J[i-1], jbylayer[ll-i])))
    end

    return reverse(J)
end

function BasesOfSocleSeries(data::Dict)
    lay = reverse(data[:layers])
    ll = data[:LL]
    sbylayer = map(i -> findall(isequal(i), lay), 0:(ll-2))
    S = [sbylayer[1]]
    for i in 2:(ll-1)
      push!(S, sort(union(S[i-1], sbylayer[i])))
    end

    return S
end

function BasisOfPowers(data::Dict, I::Vector{Int}, p::Int, m::Int)
    result = Int[]
    q = data[:parameters][1]
    pm = BigInt(p)^m
    quo = div(q-1, pm)
    for k in I
      if maximum(data[:monomials][k]) <= quo
        push!(result, pm * (k-1) + 1)
      end
    end

    return result
end

function BasisOfPMRoots(data::Dict, I::Vector{Int}, p::Int, m::Int)
    q= data[:parameters][1]
    pm = BigInt(p)^m
    quo = div(q-1, pm)

    # The first condition means that B_k^{p^m} is zero
    # and hence inside the given subspace.
    # For the second condition, we may assume that B_k^{p^m}
    # is nonzero and hence equal to B_{p^m (k-1)+1}.
    return filter(k -> maximum(data[:monomials][k]) > quo ||
                        pm * (k-1) + 1 in I,
                  1:length(data[:monomials]))
end

IsSubset(I::Vector{Int}, J::Vector{Int}) = all(x -> x in I, J)

function BasisOfPC(data::Dict, I::Vector{Int}, J::Vector{Int})
    # Note that we are interested in subspaces of the radical.
    return filter(k -> IsSubset(J, BasisOfProductSpace(data, [k], I)),
                  2:length(data[:monomials]))
end

max_stars_repo_name: oscar-system/SingerAlg

max_stars_repo_licenses: Naumen Condor-1.1 MS-PL

max_issues_repo_name: oscar-system/SingerAlg

max_issues_repo_licenses: Naumen Condor-1.1 MS-PL

max_forks_repo_name: oscar-system/SingerAlg

max_forks_repo_licenses: Naumen Condor-1.1 MS-PL