text:
@testset "Testing forward on trivial QP" begin
    # using example on https://osqp.org/docs/examples/setup-and-solve.html
    Q = [
        4.0 1.0
        1.0 2.0
    ]
    q = [1.0, 1.0]
    G = [
         1.0 1.0
         1.0 0.0
         0.0 1.0
        -1.0 -1.0
        -1.0 0.0
        0.0 -1.0
    ]
    h = [1, 0.7, 0.7, -1, 0, 0]

    model = diff_optimizer(Ipopt.Optimizer)
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variables(model, 2)

    # define objective
    quad_terms = MOI.ScalarQuadraticTerm{Float64}[]
    for i in 1:2
        for j in i:2 # indexes (i,j), (j,i) will be mirrored. specify only one kind
            push!(
                quad_terms,
                MOI.ScalarQuadraticTerm(Q[i,j],x[i],x[j])
            )
        end
    end

    objective_function = MOI.ScalarQuadraticFunction(
                            MOI.ScalarAffineTerm.(q, x),
                            quad_terms,
                            0.0
                        )
    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(), objective_function)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    # add constraints
    for i in 1:6
        MOI.add_constraint(
            model,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(G[i,:], x), 0.0),
            MOI.LessThan(h[i])
        )
    end

    MOI.optimize!(model)

    @test model.primal_optimal ≈ [0.3, 0.7] atol=ATOL rtol=RTOL
end

@testset "Differentiating trivial QP 1" begin
    Q = [
        4.0 1.0
        1.0 2.0
    ]
    q = [1.0, 1.0]
    G = [1.0 1.0]
    h = [-1.0]

    model = diff_optimizer(OSQP.Optimizer)
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variables(model, 2)

    # define objective
    quad_terms = MOI.ScalarQuadraticTerm{Float64}[]
    for i in 1:2
        for j in i:2 # indexes (i,j), (j,i) will be mirrored. specify only one kind
            push!(
                quad_terms,
                MOI.ScalarQuadraticTerm(Q[i,j], x[i], x[j])
            )
        end
    end

    objective_function = MOI.ScalarQuadraticFunction(
                            MOI.ScalarAffineTerm.(q, x),
                            quad_terms,
                            0.0,
                        )
    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(), objective_function)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    # add constraint
    MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(G[1, :], x), 0.0),
        MOI.LessThan(h[1])
    )

    MOI.optimize!(model)

    @test model.primal_optimal ≈ [-0.25; -0.75] atol=ATOL rtol=RTOL

    grad_wrt_h = backward!(model, ["h"], ones(2))[1]

    @test grad_wrt_h ≈ [1.0] atol=2ATOL rtol=RTOL
end

# @testset "Differentiating a non-convex QP" begin
#     Q = [0.0 0.0; 1.0 2.0]
#     q = [1.0; 1.0]
#     G = [1.0 1.0;]
#     h = [-1.0;]

#     model = diff_optimizer(OSQP.Optimizer)
#     x = MOI.add_variables(model, 2)

#     # define objective
#     quad_terms = MOI.ScalarQuadraticTerm{Float64}[]
#     for i in 1:2
#         for j in i:2 # indexes (i,j), (j,i) will be mirrored. specify only one kind
#             push!(
#                 quad_terms,
#                 MOI.ScalarQuadraticTerm(Q[i,j], x[i], x[j]),
#             )
#         end
#     end

#     objective_function = MOI.ScalarQuadraticFunction(
#                             MOI.ScalarAffineTerm.(q, x),
#                             quad_terms,
#                             0.0,
#                         )
#     MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(), objective_function)
#     MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

#     # add constraint
#     MOI.add_constraint(
#         model,
#         MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(G[1, :], x), 0.0),
#         MOI.LessThan(h[1]),
#     )

#     @test_throws ErrorException MOI.optimize!(model) # should break
# end

@testset "Differentiating QP with inequality and equality constraints" begin
    # refered from: https://www.mathworks.com/help/optim/ug/quadprog.html#d120e113424
    # Find equivalent qpth program here - https://github.com/AKS1996/jump-gsoc-2020/blob/master/DiffOpt_tests_4_py.ipynb

    Q = [
         1.0 -1.0 1.0;
        -1.0  2.0 -2.0;
         1.0 -2.0 4.0
    ]
    q = [2.0, -3.0, 1.0]
    G = [0.0 0.0 1.0
         0.0 1.0 0.0
         1.0 0.0 0.0
         0.0 0.0 -1.0
         0.0 -1.0 0.0
         -1.0 0.0 0.0
    ]
    h = [1.0, 1.0, 1.0, 0.0, 0.0, 0.0,]
    A = [1.0 1.0 1.0;]
    b = [0.5]

    model = diff_optimizer(Ipopt.Optimizer)
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variables(model, 3)

    # define objective
    quad_terms = MOI.ScalarQuadraticTerm{Float64}[]
    for i in 1:3
        for j in i:3 # indexes (i,j), (j,i) will be mirrored. specify only one kind
            push!(
                quad_terms,
                MOI.ScalarQuadraticTerm(Q[i,j], x[i], x[j])
            )
        end
    end

    objective_function = MOI.ScalarQuadraticFunction(
                            MOI.ScalarAffineTerm.(q, x),
                            quad_terms,
                            0.0
                        )
    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(), objective_function)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    # add constraint
    for i in 1:6
        MOI.add_constraint(
            model,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(G[i, :], x), 0.0),
            MOI.LessThan(h[i])
        )
    end

    for i in 1:1
        MOI.add_constraint(
            model,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(A[i,:], x), 0.0),
            MOI.EqualTo(b[i])
        )
    end

    MOI.optimize!(model)

    z = model.primal_optimal

    @test z ≈ [0.0, 0.5, 0.0] atol=ATOL rtol=RTOL

    grads = backward!(model, ["Q","q","G","h","A","b"], ones(3))

    dl_dQ = grads[1]
    dl_dq = grads[2]
    dl_dG = grads[3]
    dl_dh = grads[4]
    dl_dA = grads[5]
    dl_db = grads[6]

    @test dl_dQ ≈ zeros(3,3)  atol=ATOL rtol=RTOL

    @test dl_dq ≈ zeros(3,1) atol=ATOL rtol=RTOL

    @test dl_dG ≈ zeros(6,3) atol=ATOL rtol=RTOL

    @test dl_dh ≈ zeros(6,1) atol=ATOL rtol=RTOL

    @test dl_dA ≈ [0.0 -0.5 0.0] atol=ATOL rtol=RTOL

    @test dl_db ≈ [1.0] atol=ATOL rtol=RTOL
end

# refered from https://github.com/jump-dev/MathOptInterface.jl/blob/master/src/Test/contquadratic.jl#L3
# Find equivalent CVXPYLayers and QPTH code here:
#               https://github.com/AKS1996/jump-gsoc-2020/blob/master/DiffOpt_tests_1_py.ipynb
@testset "Differentiating MOI examples 1" begin
    # homogeneous quadratic objective
    # Min x^2 + xy + y^2 + yz + z^2
    # st  x + 2y + 3z >= 4 (c1)
    #     x +  y      >= 1 (c2)
    #     x, y, z \in R

    model = diff_optimizer(OSQP.Optimizer)
    MOI.set(model, MOI.Silent(), true)
    v = MOI.add_variables(model, 3)
    @test MOI.get(model, MOI.NumberOfVariables()) == 3

    c1 = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([-1.0, -2.0, -3.0], v),
            0.0),
        MOI.LessThan(-4.0)
    )
    c2 = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([-1.0, -1.0, 0.0], v),
            0.0),
        MOI.LessThan(-1.0)
    )

    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MIN_SENSE

    obj = MOI.ScalarQuadraticFunction(
        MOI.ScalarAffineTerm{Float64}[],
        MOI.ScalarQuadraticTerm.(
            [2.0, 1.0, 2.0, 1.0, 2.0],
            v[[1, 1, 2, 2, 3]],
            v[[1, 2, 2, 3, 3]]
        ),
        0.0,
    )
    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(), obj)

    MOI.optimize!(model)

    z = model.primal_optimal

    @test z ≈ [4/7, 3/7, 6/7] atol=ATOL rtol=RTOL

    # obtain gradients
    grads = backward!(model, ["Q","q","G","h"], ones(3))

    dl_dQ = grads[1]
    dl_dq = grads[2]
    dl_dG = grads[3]
    dl_dh = grads[4]

    @test dl_dQ ≈ [-0.12244895  0.01530609 -0.11224488;
                    0.01530609  0.09183674  0.07653058;
                   -0.11224488  0.07653058 -0.06122449]  atol=ATOL rtol=RTOL

    @test dl_dq ≈ [-0.2142857;  0.21428567; -0.07142857] atol=ATOL rtol=RTOL

    @test dl_dG ≈ [0.05102035   0.30612245  0.255102;
                   0.06122443   0.36734694  0.3061224] atol=ATOL rtol=RTOL

    @test dl_dh ≈ [-0.35714284; -0.4285714] atol=ATOL rtol=RTOL
end

# refered from https://github.com/jump-dev/MathOptInterface.jl/blob/master/src/Test/contquadratic.jl#L3
# Find equivalent CVXPYLayers and QPTH code here:
#               https://github.com/AKS1996/jump-gsoc-2020/blob/master/DiffOpt_tests_2_py.ipynb
@testset "Differentiating MOI examples 2 - non trivial backward pass vector" begin
    # non-homogeneous quadratic objective
    #    minimize 2 x^2 + y^2 + xy + x + y
    #       s.t.  x, y >= 0
    #             x + y = 1

    model = diff_optimizer(Ipopt.Optimizer)
    MOI.set(model, MOI.Silent(), true)

    x = MOI.add_variable(model)
    y = MOI.add_variable(model)

    c1 = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0,1.0], [x,y]), 0.0),
        MOI.EqualTo(1.0)
    )

    vc1 = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-1.0,0.0], [x,y]), 0.0),
        MOI.LessThan(0.0)
    )
    @test vc1.value ≈ x.value atol=ATOL rtol=RTOL

    vc2 = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([0.0,-1.0], [x,y]), 0.0),
        MOI.LessThan(0.0)
    )
    @test vc2.value ≈ y.value atol=ATOL rtol=RTOL

    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    obj = MOI.ScalarQuadraticFunction(
        MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]),
        MOI.ScalarQuadraticTerm.([4.0, 2.0, 1.0], [x, y, x], [x, y, y]),
        0.0
    )
    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(), obj)

    MOI.optimize!(model)

    z = model.primal_optimal
    ν = model.dual_optimal[1]   # can be accessed in the order constraints were added
    λ = model.dual_optimal[2:3]

    @test z ≈ [0.25, 0.75] atol=ATOL rtol=RTOL
    @test λ ≈ [0.0, 0.0]   atol=ATOL rtol=RTOL
    @test ν ≈ 11/4       atol=ATOL rtol=RTOL

    # obtain gradients
    dl_dz = [1.3, 0.5]   # choosing a non trivial backward pass vector
    grads = backward!(model, ["Q", "q", "G", "h", "A", "b"], dl_dz)

    dl_dQ = grads[1]
    dl_dq = grads[2]
    dl_dG = grads[3]
    dl_dh = grads[4]
    dl_dA = grads[5]
    dl_db = grads[6]

    @test dl_dQ ≈ [-0.05   -0.05;
                   -0.05    0.15]  atol=ATOL rtol=RTOL

    @test dl_dq ≈ [-0.2; 0.2] atol=ATOL rtol=RTOL

    @test dl_dG ≈ [1e-8  1e-8; 1e-8 1e-8] atol=ATOL rtol=RTOL

    @test dl_dh ≈ [1e-8; 1e-8] atol=ATOL rtol=RTOL

    @test dl_dA ≈ [0.375 -1.075] atol=ATOL rtol=RTOL

    @test dl_db ≈ [0.7] atol=ATOL rtol=RTOL
end

@testset "Differentiating non trivial convex QP MOI" begin
    nz = 10
    nineq = 25
    neq = 10

    # read matrices from files
    names = ["P", "q", "G", "h", "A", "b"]
    matrices = []

    for name in names
        push!(matrices, readdlm(Base.Filesystem.abspath(Base.Filesystem.joinpath("data",name*".txt")), ' ', Float64, '\n'))
    end

    Q, q, G, h, A, b = matrices
    q = vec(q)
    h = vec(h)
    b = vec(b)

    model = diff_optimizer(Ipopt.Optimizer)
    MOI.set(model, MOI.Silent(), true)

    x = MOI.add_variables(model, nz)

    # define objective
    quadratic_terms = MOI.ScalarQuadraticTerm{Float64}[]
    for i in 1:nz
        for j in i:nz # indexes (i,j), (j,i) will be mirrored. specify only one kind
            push!(quadratic_terms, MOI.ScalarQuadraticTerm(Q[i,j], x[i], x[j]))
        end
    end

    objective_function = MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm.(q, x), quadratic_terms, 0.0)
    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(), objective_function)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    # set constraints
    for i in 1:nineq
        MOI.add_constraint(
            model,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(G[i,:], x), 0.0),MOI.LessThan(h[i])
        )
    end

    for i in 1:neq
        MOI.add_constraint(
            model,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(A[i,:], x), 0.0), MOI.EqualTo(b[i]),
        )
    end

    MOI.optimize!(model)

    # obtain gradients
    grads = backward!(model, ["Q", "q", "G", "h", "A", "b"], ones(nz))  # using dl_dz=[1,1,1,1,1,....]

    # read gradients from files
    names = ["dP", "dq", "dG", "dh", "dA", "db"]
    grads_actual = []

    for name in names
        push!(grads_actual, readdlm(Base.Filesystem.abspath(Base.Filesystem.joinpath("data",name*".txt")), ' ', Float64, '\n'))
    end

    grads_actual[2] = vec(grads_actual[2])
    grads_actual[4] = vec(grads_actual[4])
    grads_actual[6] = vec(grads_actual[6])

    # testing differences
    for i in 1:size(grads)[1]
        @test grads[i] ≈  grads_actual[i] atol=1e-2 rtol=1e-2
    end
end

@testset "Differentiating LP; checking gradients for non-active contraints" begin
    # Issue #40 from Gurobi.jl
    # min  x
    # s.t. x >= 0
    #      x >= 3

    model = diff_optimizer(Clp.Optimizer)
    MOI.set(model, MOI.Silent(), true)

    x = MOI.add_variables(model, 1)

    # define objective
    objective_function = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], x), 0.0)
    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), objective_function)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    # set constraints
    MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-1.0], x), 0.),
        MOI.LessThan(0.0)
    )
    MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-1.0], x), 0.),
        MOI.LessThan(-3.0)
    )

    MOI.optimize!(model)

    # obtain gradients
    grads = backward!(model, ["G", "h"], [1.0])

    @test grads[1] ≈ [0.0, 3.0] atol=ATOL rtol=RTOL
    @test grads[2] ≈ [0.0, -1.0] atol=ATOL rtol=RTOL
end

@testset "Differentiating LP; checking gradients for non-active contraints" begin
    # refered from - https://en.wikipedia.org/wiki/Simplex_algorithm#Example

    # max 2x + 3y + 4z
    # s.t. 3x+2y+z <= 10
    #      2x+5y+3z <= 15
    #      x,y,z >= 0

    model = diff_optimizer(SCS.Optimizer)
    MOI.set(model, MOI.Silent(), true)
    v = MOI.add_variables(model, 3)

    # define objective
    objective_function = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-2.0, -3.0, -4.0], v), 0.0)
    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), objective_function)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    # set constraints
    MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([3.0, 2.0, 1.0], v), 0.),
        MOI.LessThan(10.0)
    )
    MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 5.0, 3.0], v), 0.),
        MOI.LessThan(15.0)
    )
    MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-1.0, 0.0, 0.0], v), 0.),
        MOI.LessThan(0.0)
    )
    MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([0.0, -1.0, 0.0], v), 0.),
        MOI.LessThan(0.0)
    )
    MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([0.0, 0.0, -1.0], v), 0.),
        MOI.LessThan(0.0)
    )

    MOI.optimize!(model)

    # obtain gradients
    grads = backward!(model, ["Q", "q", "G", "h"], ones(3))  # using dl_dz=[1,1,1,1,1,....]

    @test grads[1] ≈ zeros(3,3) atol=ATOL rtol=RTOL
    @test grads[2] ≈ zeros(3) atol=ATOL rtol=RTOL
    @test grads[3] ≈ [0.0 0.0 0.0;
                    0.0 0.0 -5/3;
                    0.0 0.0 5/3;
                    0.0 0.0 -10/3;
                    0.0 0.0 0.0]   atol=ATOL rtol=RTOL
    @test grads[4] ≈ [0.0; 1/3; -1/3; 2/3; 0.0]   atol=ATOL rtol=RTOL
end

@testset "Differentiating simple SOCP" begin
    # referred from https://github.com/jump-dev/MathOptInterface.jl/blob/master/src/Test/contconic.jl#L789
    # find equivalent diffcp python program here: https://github.com/AKS1996/jump-gsoc-2020/blob/master/diffcp_socp_1_py.ipynb

    model = diff_optimizer(SCS.Optimizer)
    MOI.set(model, MOI.Silent(), true)
    x,y,t = MOI.add_variables(model, 3)

    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    ceq  = MOI.add_constraint(model, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(-1.0, t))], [1.0]), MOI.Zeros(1))
    cnon = MOI.add_constraint(model, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, y))], [-1/√2]), MOI.Nonnegatives(1))
    csoc = MOI.add_constraint(model, MOI.VectorAffineFunction(MOI.VectorAffineTerm.([1,2,3], MOI.ScalarAffineTerm.(1.0, [t,x,y])), zeros(3)), MOI.SecondOrderCone(3))

    MOI.optimize!(model)

    x = model.primal_optimal
    s = MOI.get(model, MOI.ConstraintPrimal(), model.con_idx)
    y = model.dual_optimal

    # these matrices are benchmarked with the output generated by diffcp
    # refer the python file mentioned above to get equivalent python source code
    @test x ≈ [-0.707107; 0.707107; 1.0] atol=ATOL rtol=RTOL
    @test s ≈ [[3.58469e-18], [2.62755e-17], [1.0, -0.707107, 0.707107]] atol=ATOL rtol=RTOL
    @test y ≈ [[1.41421], [1.0], [1.41421, 1.0, -1.0]] atol=ATOL rtol=RTOL

    dA = Matrix{Float64}(I, 5, 3)
    db = zeros(5)
    dc = zeros(3)

    dx, dy, ds = backward_conic!(model, dA, db, dc)

    @test dx ≈ [1.12132144; 0.707107; 0.70710656] atol=ATOL rtol=RTOL
    @test ds ≈ [0.0; 0.0; -2.92893438e-01;  1.12132144e+00; 7.07106999e-01]  atol=ATOL rtol=RTOL
    @test dy ≈ [2.4142175;   5.00000557;  3.8284315;   1.414214;   -4.00000495] atol=ATOL rtol=RTOL
end

@testset "Differentiating simple PSD program" begin
    # refered from https://github.com/jump-dev/MathOptInterface.jl/blob/master/src/Test/contconic.jl#L2339
    # find equivalent diffcp program here: https://github.com/AKS1996/jump-gsoc-2020/blob/master/diffcp_sdp_1_py.ipynb

    model = diff_optimizer(SCS.Optimizer)
    MOI.set(model, MOI.Silent(), true)
    X = MOI.add_variables(model, 3)
    vov = MOI.VectorOfVariables(X)
    cX = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction{Float64}(vov),
        MOI.PositiveSemidefiniteConeTriangle(2)
    )

    c  = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, X[2]))],
            [-1.0]
        ),
        MOI.Zeros(1)
    )

    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, [X[1], X[end]]), 0.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    MOI.optimize!(model)

    x = model.primal_optimal
    s = MOI.get(model, MOI.ConstraintPrimal(), model.con_idx)
    y = model.dual_optimal

    @test x ≈ ones(3) atol=ATOL rtol=RTOL
    @test s ≈ [ones(3), [0.0]] atol=ATOL rtol=RTOL
    @test y ≈ [[1.0, -1.0, 1.0], [2.0]]  atol=ATOL rtol=RTOL

    # test1: changing the constant in `c`, i.e. changing value of X[2]
    dA = zeros(4, 3)
    db = zeros(4)
    db[4] = 1.0
    dc = zeros(3)

    dx, dy, ds = backward_conic!(model, dA, db, dc)

    @test dx ≈ -ones(3) atol=ATOL rtol=RTOL  # will change the value of other 2 variables
    @test ds[1:3] ≈ -ones(3)  atol=ATOL rtol=RTOL  # will affect PSD constraint too

    # test2: changing X[1], X[3] but keeping the objective (their sum) same
    dA = zeros(4, 3)
    db = zeros(4)
    dc = zeros(3)
    dc[1] = -1.0
    dc[3] = 1.0

    dx, dy, ds = backward_conic!(model, dA, db, dc)

    @test dx ≈ [1.0, 0.0, -1.0] atol=ATOL rtol=RTOL  # note: no effect on X[2]
end

@testset "Differentiating conic with PSD and SOC constraints" begin
    # similar to https://github.com/jump-dev/MathOptInterface.jl/blob/master/src/Test/contconic.jl#L2417
    # find equivalent diffcp example here - https://github.com/AKS1996/jump-gsoc-2020/blob/master/diffcp_sdp_2_py.ipynb

    model = diff_optimizer(SCS.Optimizer)
    MOI.set(model, MOI.Silent(), true)

    δ = √(1 + (3*√2+2)*√(-116*√2+166) / 14) / 2
    ε = √((1 - 2*(√2-1)*δ^2) / (2-√2))
    y2 = 1 - ε*δ
    y1 = 1 - √2*y2
    obj = y1 + y2/2
    k = -2*δ/ε
    x2 = ((3-2obj)*(2+k^2)-4) / (4*(2+k^2)-4*√2)
    α = √(3-2obj-4x2)/2
    β = k*α

    X = MOI.add_variables(model, 6)
    x = MOI.add_variables(model, 3)

    vov = MOI.VectorOfVariables(X)
    cX = MOI.add_constraint(model, MOI.VectorAffineFunction{Float64}(vov), MOI.PositiveSemidefiniteConeTriangle(3))
    cx = MOI.add_constraint(model, MOI.VectorAffineFunction{Float64}(MOI.VectorOfVariables(x)), MOI.SecondOrderCone(3))

    c1 = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(1:1, MOI.ScalarAffineTerm.([1., 1., 1., 1.], [X[1], X[3], X[end], x[1]])),
            [-1.0]
        ),
        MOI.Zeros(1)
    )
    c2 = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(1:1, MOI.ScalarAffineTerm.([1., 2, 1, 2, 2, 1, 1, 1], [X; x[2]; x[3]])),
            [-0.5]
        ),
        MOI.Zeros(1),
    )

    # this is a useless constraint - refer the tests below
    # even if we comment this, it won't affect the optimal values
    c_extra = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(1:1, MOI.ScalarAffineTerm.(ones(3), x)),
            [100.0]
        ),
        MOI.Nonnegatives(1)
    )

    objXidx = [1:3; 5:6]
    objXcoefs = 2*ones(5)
    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
    MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([objXcoefs; 1.0], [X[objXidx]; x[1]]), 0.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    MOI.optimize!(model)

    x = model.primal_optimal
    s = MOI.get(model, MOI.ConstraintPrimal(), model.con_idx)
    y = model.dual_optimal

    @test x ≈ [ 0.21725121; -0.25996907;  0.31108582;  0.21725009; -0.25996907;  0.21725121;
                0.2544097;   0.17989425;  0.17989425] atol=ATOL rtol=RTOL
    @test s ≈ [[0.217251, -0.25997, 0.31109, 0.217251, -0.25997, 0.217251],
                [0.254408, 0.179894, 0.179894], [0.0], [0.0], [100.614]] atol=ATOL rtol=RTOL   # TODO: it should be 100, not 100.614, its surely a residual error
    @test y ≈ [[1.13334, 0.678095, 1.13334, -0.321905, 0.678095, 1.13334],
                [0.455242, -0.321905, -0.321905], [0.544758], [0.321905], [0.0]]  atol=ATOL rtol=RTOL

    # test c_extra
    dA = zeros(12, 9)
    db = zeros(12)
    db[12] = 1.0
    dc = zeros(9)

    dx, dy, ds = backward_conic!(model, dA, db, dc)

    # a small change in the constant in c_extra should not affect any other variable or constraint other than c_extra itself
    @test dx ≈ zeros(9) atol=1e-2
    @test dy ≈ zeros(12) atol=1e-2
    @test ds[1:end-1] ≈ zeros(11) atol=1e-2
    @test ds[end] ≈ 1.0 atol=1e-2   # except c_extra itself
end

@testset "Differentiating conic with PSD and POS constraints" begin
    # refer https://github.com/jump-dev/MathOptInterface.jl/blob/master/src/Test/contconic.jl#L2575
    # find equivalent diffcp program here - https://github.com/AKS1996/jump-gsoc-2020/blob/master/diffcp_sdp_3_py.ipynb

    model = diff_optimizer(SCS.Optimizer)
    MOI.set(model, MOI.Silent(), true)

    x = MOI.add_variables(model, 7)
    @test MOI.get(model, MOI.NumberOfVariables()) == 7

    η = 10.0

    c1  = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(1, MOI.ScalarAffineTerm.(-1.0, x[1:6])),
            [η]
        ),
        MOI.Nonnegatives(1)
    )
    c2 = MOI.add_constraint(model, MOI.VectorAffineFunction(MOI.VectorAffineTerm.(1:6, MOI.ScalarAffineTerm.(1.0, x[1:6])), zeros(6)), MOI.Nonnegatives(6))
    α = 0.8
    δ = 0.9
    c3 = MOI.add_constraint(model, MOI.VectorAffineFunction(MOI.VectorAffineTerm.([fill(1, 7); fill(2, 5);     fill(3, 6)],
                                                            MOI.ScalarAffineTerm.(
                                                            [ δ/2,       α,   δ, δ/4, δ/8,      0.0, -1.0,
                                                                -δ/(2*√2), -δ/4, 0,     -δ/(8*√2), 0.0,
                                                                δ/2,     δ-α,   0,      δ/8,      δ/4, -1.0],
                                                            [x[1:7];     x[1:3]; x[5:6]; x[1:3]; x[5:7]])),
                                                            zeros(3)), MOI.PositiveSemidefiniteConeTriangle(2))
    c4 = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(1, MOI.ScalarAffineTerm.(0.0, [x[1:3]; x[5:6]])),
            [0.0]
        ),
        MOI.Zeros(1)
    )

    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x[7])], 0.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    MOI.optimize!(model)

    x = model.primal_optimal
    s = MOI.get(model, MOI.ConstraintPrimal(), model.con_idx)
    y = model.dual_optimal

    @test x' ≈ [20/3. 0.0 10/3. 0.0 0.0 0.0 1.90192379] atol=ATOL rtol=RTOL
    @test s ≈ [[0.0],  [20/3.0,  0.0,  10/3.0,  0.0,  0.0,  0.0],[4.09807621, -2.12132,  1.09807621], [0.0]] atol=ATOL rtol=RTOL
    @test y ≈ [[0.19019238], [0., 0.12597667, 0., 0.14264428, 0.14264428, 0.01274047],[0.21132487, 0.408248, 0.78867513] ,[0.0]] atol=ATOL rtol=RTOL

    dA = ones(11, 7)
    db = ones(11)
    dc = ones(7)

    dx, dy, ds = backward_conic!(model, dA, db, dc)

    atol = 0.3
    rtol = 0.01

    # compare these with https://github.com/AKS1996/jump-gsoc-2020/blob/master/diffcp_sdp_3_py.ipynb
    # results are not exactly as: 1. there is some residual error   2. diffcp results are SCS specific, hence scaled
    @test dx ≈ [-39.6066, 10.8953, -14.9189, 10.9054, 10.883, 10.9118, -21.7508] atol=atol rtol=rtol
    @test dy ≈ [-3.56905, 0.0, -0.380035, 0.0, -0.41398, -0.385321, -0.00743119, -0.644986, -0.550542, -2.36765, 0.0] atol=atol rtol=rtol
    @test ds ≈ [0.0, -50.4973, 0.0, -25.8066, 0.0, 0.0, 0.0, -7.96528, -1.62968, -2.18925, 0.0] atol=atol rtol=rtol

    # TODO: future example, how to differentiate wrt a specific constraint/variable, refer QPLib article for more
    dA = zeros(11, 7)
    dA[2:7, 1:6] = Matrix{Float64}(LinearAlgebra.I, 6, 6)  # differentiating only wrt POS constraint c2
    db = zeros(11)
    dc = zeros(7)

    dx, dy, ds = backward_conic!(model, dA, db, dc)

    # note that there's no change in the PSD slack values or dual optimas
    @test dy ≈ [0.0, 0.0, 0.125978, 0.0, 0.142644, 0.142641, 0.0127401, 0.0, 0.0, 0.0, 0.0] atol=atol rtol=RTOL
    @test ds ≈ [0.0, -6.66672, 0.0, -3.33336, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] atol=atol rtol=RTOL
end

@testset "Differentiating a simple PSD" begin
    # refer https://github.com/jump-dev/MathOptInterface.jl/blob/master/src/Test/contconic.jl#L2643
    # find equivalent diffcp program here - https://github.com/AKS1996/jump-gsoc-2020/blob/master/diffcp_sdp_0_py.ipynb

    model = DiffOpt.diff_optimizer(SCS.Optimizer)
    MOI.set(model, MOI.Silent(), true)

    x = MOI.add_variable(model)
    fx = MOI.SingleVariable(x)

    func = MOIU.operate(vcat, Float64, fx, one(Float64), fx, one(Float64), one(Float64), fx)

    c = MOI.add_constraint(model, func, MOI.PositiveSemidefiniteConeTriangle(3))

    # MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(), MOI.SingleVariable(x))
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, [x]), 0.0)
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    MOI.optimize!(model)

    x = model.primal_optimal
    s = MOI.get(model, MOI.ConstraintPrimal(), model.con_idx)
    y = model.dual_optimal

    @test x ≈ [1.0] atol=ATOL rtol=RTOL
    @test s ≈ [ones(6)] atol=ATOL rtol=RTOL
    @test y ≈ [[1/3.,  -1/6.,  1/3.,  -1/6.,  -1/6.,  1/3.]]  atol=ATOL rtol=RTOL

    # SCS/Mosek specific
    # @test s' ≈ [1.         1.41421356 1.41421356 1.         1.41421356 1.        ] atol=ATOL rtol=RTOL
    # @test y' ≈ [ 0.33333333 -0.23570226 -0.23570226  0.33333333 -0.23570226  0.33333333]  atol=ATOL rtol=RTOL

    dA = zeros(6, 1)
    db = ones(6)
    dc = zeros(1)

    dx, dy, ds = backward_conic!(model, dA, db, dc)

    @test dx ≈ [-0.5] atol=ATOL rtol=RTOL
    @test dy ≈ zeros(6) atol=ATOL rtol=RTOL
    @test ds ≈ [0.5, 1.0, 0.5, 1.0, 1.0, 0.5] atol=ATOL rtol=RTOL

    # test 2
    dA = zeros(6, 1)
    db = zeros(6)
    dc = ones(1)

    dx, dy, ds = backward_conic!(model, dA, db, dc)

    @test dx ≈ zeros(1) atol=ATOL rtol=RTOL
    @test dy ≈ [0.333333, -0.333333, 0.333333, -0.333333, -0.333333, 0.333333] atol=ATOL rtol=RTOL
    @test ds ≈ zeros(6) atol=ATOL rtol=RTOL
end


text:
export wcpair

wcpair = UInt8[
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, '-', 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 'T', 'V', 'G', 'H', 0, 0, 'C', 'D', 0, 0, 'M', 0, 'K', 'N', 0,
    0, 0, 'Y', 'S', 'A', 0, 'B', 'W', 0, 'R', 'Z', 0, 0, 0, 0, 0,
    0, 't', 'v', 'g', 'h', 0, 0, 'c', 'd', 0, 0, 'm', 0, 'k', 'n', 0,
    0, 0, 'y', 's', 'a', 0, 'b', 'w', 0, 'r', 'z', 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
]

ext: jl

max_stars_repo_name: neherlab/pangraph

max_issues_repo_name: neherlab/pangraph

max_forks_repo_name: neherlab/pangraph


text:
# Autogenerated wrapper script for xtrx_jll for armv6l-linux-gnueabihf-cxx11
export libXTRXSupport, liblms7compact, libxtrx, libxtrxdsp, libxtrxll, test_xtrx, test_xtrxflash, test_xtrxll

using soapysdr_jll
JLLWrappers.@generate_wrapper_header("xtrx")
JLLWrappers.@declare_library_product(libXTRXSupport, "libXTRXSupport.so")
JLLWrappers.@declare_library_product(liblms7compact, "liblms7compact.so")
JLLWrappers.@declare_library_product(libxtrx, "libxtrx.so.0")
JLLWrappers.@declare_library_product(libxtrxdsp, "libxtrxdsp.so.0")
JLLWrappers.@declare_library_product(libxtrxll, "libxtrxll.so.0")
JLLWrappers.@declare_executable_product(test_xtrx)
JLLWrappers.@declare_executable_product(test_xtrxflash)
JLLWrappers.@declare_executable_product(test_xtrxll)
function __init__()
    JLLWrappers.@generate_init_header(soapysdr_jll)
    JLLWrappers.@init_library_product(
        libXTRXSupport,
        "lib/SoapySDR/modules0.8/libXTRXSupport.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

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

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

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

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

    JLLWrappers.@init_executable_product(
        test_xtrx,
        "bin/test_xtrx",
    )

    JLLWrappers.@init_executable_product(
        test_xtrxflash,
        "bin/test_xtrxflash",
    )

    JLLWrappers.@init_executable_product(
        test_xtrxll,
        "bin/test_xtrxll",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using MakieLayout
using Makie

begin
    scene = Scene(resolution = (1000, 1000));
    screen = display(scene)
    campixel!(scene);

    maingl = GridLayout(scene, alignmode = Outside(30))

    las = Array{LAxis, 2}(undef, 4, 4)

    for i in 1:4, j in 1:4
        las[i, j] = maingl[i, j] = LAxis(scene)
    end

    lt = maingl[0, :] = LText(scene, text="Suptitle", textsize=50)
    lt2 = maingl[2:5, 5] = LText(scene, text="Side Title", textsize=50,
        rotation = -pi/2)
    lt3 = maingl[6, 1:4] = LText(scene, text="Sub Title", textsize=50)
    lt3 = maingl[2:5, 0] = LText(scene, text="Left Title", textsize=50,
        rotation=pi/2)
end


text:
# TODO: choose method based on problem size
defaultmethod(x) = :Logarithmic

type PWLData
    counter::Int
    PWLData() = new(0)
end

function initPWL!(m::JuMP.Model)
    if !haskey(m.ext, :PWL)
        m.ext[:PWL] = PWLData()
    end
    nothing
end

function piecewiselinear(m::JuMP.Model, x::JuMP.Variable, d, f::Function; method=defaultmethod(d))
    initPWL!(m)
    fd = [f(xx) for xx in d]
    piecewiselinear(m, x, d, fd; method=method)
end

piecewiselinear(m::JuMP.Model, x::JuMP.Variable, d, fd; method=defaultmethod(d)) =
    piecewiselinear(m, x, UnivariatePWLFunction(d, fd); method=method)

function piecewiselinear(m::JuMP.Model, x::JuMP.Variable, pwl::UnivariatePWLFunction; method=defaultmethod(pwl))
    initPWL!(m)
    counter = m.ext[:PWL].counter + 1
    d = [_x[1] for _x in pwl.x]
    fd = pwl.z
    n = length(d)
    z = JuMP.@variable(m, lowerbound=minimum(fd), upperbound=maximum(fd), basename="z_$counter")
    if method == :Incremental
        δ = JuMP.@variable(m, [1:n], lowerbound=0, upperbound=1, basename="δ_$counter")
        y = JuMP.@variable(m, [1:n-1], Bin, basename="y_$counter")
        JuMP.@constraint(m, x ==  d[1] + sum(δ[i]*( d[i+1]- d[i]) for i in 1:n-1))
        JuMP.@constraint(m, z == fd[1] + sum(δ[i]*(fd[i+1]-fd[i]) for i in 1:n-1))
        for i in 1:n-1
            JuMP.@constraint(m, δ[i+1] ≤ y[i])
            JuMP.@constraint(m, y[i] ≤ δ[i])
        end
    else # V-formulation method
        λ = JuMP.@variable(m, [1:n], lowerbound=0, upperbound=1, basename="λ_$counter")
        JuMP.@constraint(m, sum(λ) == 1)
        JuMP.@constraint(m, sum(λ[i]* d[i] for i in 1:n) == x)
        JuMP.@constraint(m, sum(λ[i]*fd[i] for i in 1:n) == z)
        if method == :Logarithmic
            sos2_logarthmic_formulation!(m, λ)
        elseif method == :CC
            sos2_cc_formulation!(m, λ)
        elseif method == :MC
            sos2_mc_formulation!(m, λ)
        end
    end
    m.ext[:PWL].counter = counter + 1
    z
end

function piecewiselinear(m::JuMP.Model, x::JuMP.Variable, y::JuMP.Variable, pwl::BivariatePWLFunction; method=defaultmethod(pwl))
    d = pwl.x
    fd = pwl.z
    n = length(d)
    JuMP.@variable(m, minimum(fd) ≤ z ≤ maximum(fd))
    error("Bivariate PWL functions currently not implemented")
end

function sos2_cc_formulation!(m::JuMP.Model, λ)
    counter = m.ext[:PWL].counter
    n = length(λ)-1
    y = JuMP.@variable(m, [1:n], Bin, basename="y_$counter")
    JuMP.@constraint(m, sum(y) == 1)
    for i in 1:n-1
        JuMP.@constraint(m, λ[i] ≤ y[i] + y[i+1])
    end
    nothing
end

function sos2_mc_formulation!(m::JuMP.Model, λ)
    counter = m.ext[:PWL].counter
    n = length(λ)-1
    γ = JuMP.@variable(m, [1:n,1:n+1], basename="γ_$counter")
    y = JuMP.@variable(m, [1:n], Bin,  basename="y_$counter")
    JuMP.@constraint(m, sum(γ[i,:] for i in 1:n) .== λ)
    for i in 1:n
        JuMP.@constraint(m, γ[i,i] + γ[i,i+1] ≥ y[i])
    end
    nothing
end

function reflected_gray(k::Int)
    if k == 2
        codes = [[0,0],[0,1],[1,1],[1,0]]
    elseif k < 2
        error()
    else
        codes′ = reflected_gray(k-1)
        codes = vcat([vcat(code,0) for code in codes′],
                     [vcat(code,1) for code in reverse(codes′)])
        return codes
    end
    codes
end

function sos2_logarthmic_formulation!(m::JuMP.Model, λ)
    counter = m.ext[:PWL].counter
    n = length(λ)-1
    k = ceil(Int,log2(n))
    H = reflected_gray(k)
    y = JuMP.@variable(m, [1:k], Bin, basename="y_$counter")
    for i in 1:k
        JuMP.@constraints(m, begin
            H[1][i]*λ[1] + sum(min(H[v][i],H[v-1][i])*λ[v] for v in 2:n-1) + H[n][i]*λ[n] ≤ y[i]
            H[1][i]*λ[1] + sum(max(H[v][i],H[v-1][i])*λ[v] for v in 2:n-1) + H[n][i]*λ[n] ≥ y[i]
        end)
    end
    nothing
end


text:
# This script automates the genenration of the Artifacts.toml file given a
# version of gmsh. It is essentially run once per release when the gmsh version
# is updated.

using Artifacts
using ArtifactUtils
using Pkg.BinaryPlatforms: Linux, MacOS, Windows, FreeBSD

const version = "4.8.4"

platform_to_dist = Dict(
        # glibc Linuces
        Linux(:i686,libc=:glibc)        => "Linux32",
        Linux(:armv7l,libc=:glibc)      => "Linux32",
        Linux(:x86_64,libc=:glibc)      => "Linux64",
        Linux(:aarch64,libc=:glibc)     => "Linux64",
        Linux(:powerpc64le,libc=:glibc) => "Linux64",

        # # musl Linuces
        Linux(:i686, libc=:musl)        => "Linux32",
        Linux(:x86_64, libc=:musl)      => "Linux64",
        Linux(:aarch64, libc=:musl)     => "Linux64",
        Linux(:armv7l, libc=:musl)      => "Linux32",

        # BSDs
        MacOS(:x86_64)                  => "MacOSX",
        FreeBSD(:x86_64)                => "Linux64",

        # # Windows
        Windows(:i686)                    => "Windows32",
        Windows(:x86_64)                  => "Windows64"
)

dist_to_url = Dict(
    "Linux32"   => "https://gmsh.info/bin/Linux/gmsh-$version-Linux32-sdk.tgz",
    "Linux64"   => "https://gmsh.info/bin/Linux/gmsh-$version-Linux64-sdk.tgz",
    "MacOSX"     => "https://gmsh.info/bin/MacOSX/gmsh-$version-MacOSX-sdk.tgz",
    # NOTE: the url for windows is currently set to a "personal" file on the
    # WaveProp organization. The reason is that Artifacts do not play well with
    # zipped files, and the official version of gmsh provided for windows on the
    # gmsh website is zipped. The files below are just the decompressed and
    # tarred version of the official files.
    "Windows32" => "https://github.com/WaveProp/.github/raw/main/gmsh-4.8.4-Windows32-sdk.tgz",
    "Windows64" => "https://github.com/WaveProp/.github/raw/main/gmsh-4.8.4-Windows64-sdk.tgz"
)

artifact_toml = joinpath(@__DIR__,"..","Artifacts.toml") |> normpath

artifact_name = "gmsh$version"

for (platform,dist) in platform_to_dist
    url                  = dist_to_url[dist]
    add_artifact!(artifact_toml,artifact_name,url;clear=true,platform=platform,force=true,lazy=false)
end

max_stars_repo_name: maltezfaria/GmshSDK

max_issues_repo_name: maltezfaria/GmshSDK

max_forks_repo_name: WaveProp/GmshSDK


text:
@doc """
	apply!(state::ExactState, gate::QGate, qubits::Vector{Int}; kwarg...)
Apply `gate` to `state` on given `qubits`.
# Keyword Arguments
* All following keyword arguments are passed to svd operation 1)to control the truncation accuracy 2) to specify the source of svd to call (BLAS or ITensors build-in)
- `maxdim::Int`: Maximal bond dimension for links between sites.
- `mindim::Int`: Minimal bond dimension for links between sites.
- `cutoff::Float64`: Truncation criterion.
- `absoluteCutoff::Bool`: Truncate singular values smaller than `cutoff`.
- `doRelCutoff::Bool`: Truncate the bottom `cutoff` fraction of singular values by norm.
- `fastSVD::Bool`: Calling BLAS svd if set to `false`; calling ITensors' recursive svd subroutine if set to `true`.
""" ->
apply!(state::ExactState, gate::QGate, qubits::Vector{Int}; kwarg...)


text:
using Documenter, PlotlyBase

makedocs(
    # options
    modules=[PlotlyBase],
    format=:html,
    sitename="PlotlyBase.jl"
)


text:
# Adapted from https://github.com/terasakisatoshi/sysimage_creator/
import Pkg

nthreads = Threads.nthreads()

using IJulia
kernelpath = IJulia.installkernel(
    "Julia",
    "--project=@.";
    env=Dict("JULIA_NUM_THREADS" => "$(nthreads)",)
)

max_stars_repo_name: ww-juliabook/template

max_issues_repo_name: ww-juliabook/template

max_forks_repo_name: ww-juliabook/template


text:
using LinearAlgebra
using SparseArrays

lchol(A) = LowerTriangular(Matrix(LinearAlgebra._chol!(copy(A), UpperTriangular)[1])')
lchol!(A) = LowerTriangular(Matrix(LinearAlgebra._chol!(A, UpperTriangular)[1])')

"""
    cholinverse!(L, x)
Solve L L'y =x using two backsolves,
L should be lower triangular
"""
function cholinverse!(L, x)
    naivesub!(L, x) # triangular backsolves
    naivesub!(UpperTriangular(L'), x)
    x
end

function naivesub!(At::Adjoint{<:Any,<:LowerTriangular}, b::AbstractVector, x::AbstractVector = b)
    A = At.parent
    n = size(A, 2)
    if !(n == length(b) == length(x))
        throw(DimensionMismatch("second dimension of left hand side A, $n, length of output x, $(length(x)), and length of right hand side b, $(length(b)), must be equal"))
    end
    @inbounds for j in n:-1:1
        iszero(A.data[j,j]) && throw(SingularException(j))
        xj = x[j] = A.data[j,j] \ b[j]
        for i in j-1:-1:1 # counterintuitively 1:j-1 performs slightly better
            b[i] -= A.data[j,i] * xj
        end
    end
    x
end

function naivesub!(At::Adjoint{<:Any,<:LowerTriangular}, B::AbstractMatrix, X::AbstractMatrix = B)
    A = At.parent
    n = size(A, 2)
    if !(n == size(B,1) == size(B,2) == size(X,1) == size(X,2))
        throw(DimensionMismatch())
    end
    @inbounds for k in 1:n
        for j in n:-1:1
            iszero(A.data[j,j]) && throw(SingularException(j))
            xjk = X[j,k] = A.data[j,j] \ B[j,k]
            for i in j-1:-1:1 # counterintuitively 1:j-1 performs slightly better
                B[i,k] -= A.data[j,i] * xjk
            end
        end
    end
    X
end

function naivesub!(A::LowerTriangular, B::AbstractMatrix, X::AbstractMatrix = B)
    n = size(A,2)
    if !(n == size(B,1) == size(X,1))
        throw(DimensionMismatch())
    end
    if !(size(B,2) == size(X,2))
        throw(DimensionMismatch())
    end

    @inbounds for k in 1:size(B,2)
        for j in 1:n
            iszero(A.data[j,j]) && throw(SingularException(j))
            xjk = X[j,k] = A.data[j,j] \ B[j,k]
            for i in j+1:n
                B[i,k] -= A.data[i,j] * xjk
            end
        end
    end
    X
end

function naivesub!(A::UpperTriangular, B::AbstractMatrix, X::AbstractMatrix = B)
    n = size(A, 2)
    if !(n == size(B,1) == size(X,1))
        throw(DimensionMismatch())
    end
    if !(size(B,2) == size(X,2))
        throw(DimensionMismatch())
    end

    @inbounds for k in 1:size(B, 2)
        for j in n:-1:1
            iszero(A.data[j,j]) && throw(SingularException(j))
            xjk = X[j,k] = A.data[j,j] \ B[j,k]
            for i in j-1:-1:1 # counterintuitively 1:j-1 performs slightly better
                B[i,k] -= A.data[i,j] * xjk
            end
        end
    end
    X
end


text:
#### Pattern matching
### Matching procedures
# A matcher is a function which takes 3 arguments
# 1. Expression
# 2. Dictionary
# 3. Callback: takes arguments Dictionary × Number of elements matched
#
function matcher(val::Any)
    function literal_matcher(next, data, bindings)
        islist(data) && isequal(car(data), val) ? next(bindings, 1) : nothing
    end
end

function matcher(slot::Slot)
    function slot_matcher(next, data, bindings)
        !islist(data) && return
        val = get(bindings, slot.name, nothing)
        if val !== nothing
            if isequal(val, car(data))
                return next(bindings, 1)
            end
        else
            if slot.predicate(car(data))
                next(assoc(bindings, slot.name, car(data)), 1)
            end
        end
    end
end

# returns n == offset, 0 if failed
function trymatchexpr(data, value, n)
    if !islist(value)
        return n
    elseif islist(value) && islist(data)
        if !islist(data)
            # didn't fully match
            return nothing
        end

        while isequal(car(value), car(data))
            n += 1
            value = cdr(value)
            data = cdr(data)

            if !islist(value)
                return n
            elseif !islist(data)
                return nothing
            end
        end

        return !islist(value) ? n : nothing
    elseif isequal(value, data)
        return n + 1
    end
end

function matcher(segment::Segment)
    function segment_matcher(success, data, bindings)
        val = get(bindings, segment.name, nothing)

        if val !== nothing
            n = trymatchexpr(data, val, 0)
            if n !== nothing
                success(bindings, n)
            end
        else
            res = nothing

            for i=length(data):-1:0
                subexpr = take_n(data, i)

                if segment.predicate(subexpr)
                    res = success(assoc(bindings, segment.name, subexpr), i)
                    if res !== nothing
                        break
                    end
                end
            end

            return res
        end
    end
end

function matcher(term::Term)
    matchers = (matcher(operation(term)), map(matcher, arguments(term))...,)
    function term_matcher(success, data, bindings)

        !islist(data) && return nothing
        !(istree(car(data))) && return nothing

        function loop(term, bindings′, matchers′) # Get it to compile faster
            if !islist(matchers′)
                if  !islist(term)
                    return success(bindings′, 1)
                end
                return nothing
            end
            car(matchers′)(term, bindings′) do b, n
                loop(drop_n(term, n), b, cdr(matchers′))
            end
        end

        loop(car(data), bindings, matchers) # Try to eat exactly one term
    end
end


text:
@inline function _time_cuda_elapsed(kernel::F, mem_dst, mem_src) where {F}
    t = CUDA.context!(context(mem_src)) do
        CUDA.@elapsed begin
            NVTX.@range "p2p: kernel call" begin
                kernel(mem_dst, mem_src)
            end
        end
    end
    return t
end

"""
    p2p_bandwidth([memsize::UnitPrefixedBytes]; kwargs...)

Performs a peer-to-peer memory copy benchmark (time measurement) and
returns an inter-gpu memory bandwidth estimate (in GiB/s) derived from it.

**Keyword arguments:**
* `src` (default: `0`): source device
* `dst` (default: `1`): destination device
* `nbench` (default: `5`): number of time measurements (i.e. p2p memcopies)
* `verbose` (default: `true`): set to false to turn off any printing.
* `hist` (default: `false`): when `true`, a UnicodePlots-based histogram is printed.
* `times` (default: `false`): toggle printing of measured times.
* `alternate` (default: `false`): alternate `src` and `dst`, i.e. copy data back and forth.
* `dtype` (default: `Float32`): see [`alloc_mem`](@ref).
* `io` (default: `stdout`): set the stream where the results should be printed.

**Examples:**
```julia
p2p_bandwidth()
p2p_bandwidth(MiB(1024))
p2p_bandwidth(KiB(20_000); dtype=Int32)
```
"""
function p2p_bandwidth(
    memsize::UnitPrefixedBytes=B(40_000_000);
    nbench=5,
    verbose=true,
    hist=false,
    times=false,
    alternate=false,
    dtype=Float32,
    src=0,
    dst=1,
    io::IO=stdout
)
    mem_src, mem_dst = alloc_mem(memsize; devs=(src, dst), dtype)
    actual_memsize = sizeof(mem_src)
    ts = zeros(nbench)

    NVTX.@range "p2p: nbench loop" begin
        @inbounds for i in 1:nbench
            if mod(i, alternate ? 2 : 1) == 0
                # ts[i] = _time_cuda_elapsed(copyto!, mem_dst, mem_src)
                ts[i] = _time_cuda_elapsed(
                    (x, y) -> copyto!(x, y, CUDA.stream()), mem_dst, mem_src
                )
            else
                # ts[i] = _time_cuda_elapsed(copyto!, mem_src, mem_dst)
                ts[i] = _time_cuda_elapsed(
                    (x, y) -> copyto!(x, y, CUDA.stream()), mem_dst, mem_src
                )
            end
        end
    end

    t_min = minimum(ts)
    t_max = maximum(ts)
    t_avg = mean(ts)

    actual_memsize_GiB = actual_memsize * 2^(-30)
    bws = actual_memsize_GiB ./ ts
    bw_min = minimum(bws)
    bw_max = maximum(bws)
    bw_avg = mean(bws)

    if verbose
        # println("Memsize: $(Base.format_bytes(actual_memsize))\n")
        if times
            println(io,"t_min: $t_min")
            println(io,"t_max: $t_max")
            println(io,"t_avg: $t_avg")
        end
        printstyled(io,"Bandwidth (GiB/s):\n"; bold=true)
        print(io," ├ max: ")
        printstyled(io,round(bw_max; digits=2), "\n"; color=:green, bold=true)
        println(io," ├ min: ", round(bw_min; digits=2))
        println(io," ├ avg: ", round(bw_avg; digits=2))
        print(io," └ std_dev: ")
        printstyled(io,round(std(bws); digits=2), "\n"; color=:yellow, bold=true)
        if hist
            println(io,UnicodePlots.histogram(bws; title="Bandwidths (GiB/s)", nbins=5))
        end
    end

    return bw_max
end

"""
    p2p_bandwidth_all(args...; kwargs...)

Run [`p2p_bandwidth`](@ref) for all combinations of devices.
Returns a matrix with the p2p memory bandwidth estimates.
"""
function p2p_bandwidth_all(args...;io::IO=stdout, verbose=false, kwargs...)
    ngpus = length(CUDA.devices())
    return [
        if src == dst
            nothing
        else
            p2p_bandwidth(args...; src=src, dst=dst, io=io, verbose=verbose, kwargs...)
        end for src in 0:(ngpus - 1), dst in 0:(ngpus - 1)
    ]
end

function Base.copyto!(
    dest::DenseCuArray{T}, src::DenseCuArray{T}, stream::CuStream; async=true
) where {T}
    # based on https://github.com/JuliaGPU/CUDA.jl/blob/5413070c13be8809916900205376f6c062c55a67/src/array.jl#L421
    doffs = 1
    soffs = 1
    n = length(src)
    context!(context(src)) do
        GC.@preserve src dest begin
            unsafe_copyto!(pointer(dest, doffs), pointer(src, soffs), n; async, stream)
            if Base.isbitsunion(T)
                unsafe_copyto!(
                    typetagdata(dest, doffs), typetagdata(src, soffs), n; async, stream
                )
            end
        end
    end
    return dest
end

"""
Same as [`p2p_bandwidth`](@ref) but measures the bidirectional bandwidth (copying data back and forth).
"""
function p2p_bandwidth_bidirectional(
    memsize::UnitPrefixedBytes=B(40_000_000);
    nbench=20,
    verbose=true,
    hist=false,
    times=false,
    dtype=Float32,
    dev1=0,
    dev2=1,
    repeat=100,
    io::IO=stdout
)
    mem_dev1, mem_dev2 = alloc_mem(memsize; dtype, devs=(dev1, dev2))
    actual_memsize = sizeof(mem_dev1)
    ts = zeros(nbench)

    NVTX.@range "p2p: nbench loop" begin
        @inbounds for i in 1:nbench
            ts[i] = _perform_p2p_memcpy_bidirectional(
                dev1, mem_dev1, dev2, mem_dev2; repeat
            )
        end
    end

    t_min = minimum(ts)
    t_max = maximum(ts)
    t_avg = mean(ts)

    actual_memsize_GiB = (actual_memsize) / (1024^3)
    bws = (2 * actual_memsize_GiB * repeat) ./ ts
    bw_min = minimum(bws)
    bw_max = maximum(bws)
    bw_avg = mean(bws)

    if verbose
        # println("Memsize: $(Base.format_bytes(actual_memsize))\n")
        if times
            println(io,"t_min: $t_min")
            println(io,"t_max: $t_max")
            println(io,"t_avg: $t_avg")
        end
        printstyled(io,"Bandwidth (GiB/s):\n"; bold=true)
        print(io," ├ max: ")
        printstyled(io,round(bw_max; digits=2), "\n"; color=:green, bold=true)
        println(io," ├ min: ", round(bw_min; digits=2))
        println(io," ├ avg: ", round(bw_avg; digits=2))
        print(io," └ std_dev: ")
        printstyled(io,round(std(bws); digits=2), "\n"; color=:yellow, bold=true)
        if hist
            println(io,UnicodePlots.histogram(bws; title="Bandwidths (GiB/s)", nbins=5))
        end
    end

    return bw_max
end

function _perform_p2p_memcpy_bidirectional(dev1, mem_dev1, dev2, mem_dev2; repeat=5)
    NVTX.@range "p2p: perform" begin
        start = Vector{CuEvent}(undef, 2)
        stop = Vector{CuEvent}(undef, 2)

        device!(dev1)
        stream1 = CuStream()
        device!(dev2)
        stream2 = CuStream()

        for (i, devid) in pairs((dev1, dev2))
            device!(devid)
            start[i] = CuEvent()
            stop[i] = CuEvent()
        end

        device!(dev1)
        synchronize(stream1)
        synchronize(stream2)

        # Block the stream until all the work is queued up
        # delay_flag = Ref(false)
        # device!(dev1)
        # @cuda stream=stream1 delay(delay_flag)

        # Force stream2 not to start until stream1 does, in order to ensure
        # the events on stream1 fully encompass the time needed for all
        # operations
        record(start[1], stream1)
        CUDA.wait(start[1], stream2)

        NVTX.@range "p2p: copyto!'s" begin
            for k in 1:repeat
                copyto!(mem_dev2, mem_dev1, stream2)
                copyto!(mem_dev1, mem_dev2, stream1)
            end
        end

        # Notify stream1 that stream2 is complete and record the time of
        # the total transaction
        record(stop[2], stream2)
        CUDA.wait(stop[2], stream1)
        record(stop[1], stream1)

        # Release the queued operations
        # delay_flag[] = true
        synchronize(stream1)
        synchronize(stream2)

        t = elapsed(start[1], stop[1])
    end
    return t
end

"""
Same as [`p2p_bandwidth_all`](@ref) but measures the bidirectional bandwidth (copying data back and forth).
"""
function p2p_bandwidth_bidirectional_all(args...; kwargs...)
    ngpus = length(CUDA.devices())
    return [
        if src == dst
            nothing
        else
            p2p_bandwidth_bidirectional(
                args...; dev1=src, dev2=dst, verbose=false, kwargs...
            )
        end for src in 0:(ngpus - 1), dst in 0:(ngpus - 1)
    ]
end

## delay kernel (Warning: not sure the ref flag is working as intended!)
function delay(delay_flag, timeout_clocks=10_000_000)
    # Wait until the application notifies us that it has completed queuing up the
    # experiment, or timeout and exit, allowing the application to make progress
    start_clock = CUDA.clock(UInt64)

    while !delay_flag[]
        sample_clock = CUDA.clock(UInt64)

        if (sample_clock - start_clock) > timeout_clocks
            break
        end
    end
end


text:
module AVUtil
  include(joinpath(dirname(@__FILE__),"..","..","..","init.jl"))
  w(f) = joinpath(avutil_dir, f)

  include(w("LIBAVUTIL.jl"))

  Base.zero(::Type{AVRational}) = AVRational(0, 1)

end

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
using PyPlot # would be nice to not need to load this on every proccessor
@everywhere using Shapes, RayTrace # PyPlot is slow and in the ultimate version will be removed
@everywhere begin
# constants:
const ε = 1e-4 # the ε size with which the ray advances every iteration
const lens_r = 1. # the lens radius
const maxiter = int(2π*lens_r/ε) # maximum iteration allowed. This value is based on the arc-length of an imaginary ray that travels along the surface of the lens, times two! So should be safe.
const c = [100.,-20.,330.] # location of the center of the lens. It has no real use here, but good to have it.
const retina_r = 1.2*lens_r # the retina's radius, good practice to make it a function of the lens radius
const L = Inf # the distance between the center of the lens and the source light
const nrays = 1000 # the number of discrete rays
const n_medium = 1. # the refractive index of the medium surrounding the lens
const lens_r2 = lens_r*lens_r # just to make the refractive index gradient function faster
const goal = c .- [0.,0.,lens_r] # this is where all rays that start at L=Inf should end up at.
ri(r::Float64) = sqrt(2. - r*r/lens_r2) # this is the refractive index gradient function. I set this RIG to be equal to the Luneburg lens, so that we can check the results real quick.
const lens = Lens(c,lens_r,ri,ε) # initiate the lens
const retina = Retina(c,retina_r) # initiate the retina
#const procsToUse = 2

# actual calculations

end

@everywhere function oneIter(i)
# iterate through all the rays, ploting them
    p = zeros(3,maxiter) # ray locations
    diagnose = [:Δ => 0., :len => 0, :goal => goal] # diagnosys, including the distance to the goal, the number of refraction points, the goal
    trace!(p,lens,retina,L,n_medium,maxiter,diagnose) # main function that ray-traces the system
    return (p,diagnose)
    # plotting:
    #ind = round(linspace(2,diagnose[:len],10)) # let's plot just 10 points inside the lens, no need to go crazy
    #plot3D(vec(p[1,ind]),vec(p[2,ind]),vec(p[3,ind]),color="red") # plotting it in red...
end

x = pmap(oneIter,1:nrays)

Δ = zeros(nrays) # the matrix of euclidean distance between where the ray exsits the lens and the goal
for i = 1:nrays
    p = x[i][1]
    diagnose = x[i][2]
    Δ[i] = diagnose[:Δ] # put that distance in the matrix
    ind = round(linspace(2,diagnose[:len],10)) # let's plot just 10 points inside the lens, no need to go crazy
    plot3D(vec(p[1,ind]),vec(p[2,ind]),vec(p[3,ind]),color="red") # plotting it in red...
end

# set this up ONCE for a parallel environment
#=np = nprocs()
if np >= procsToUse
	println("Already running on ",np," processors.")
else
	deltaproc = procsToUse-np
	println("Starting ",deltaproc," more processors.")
	addprocs(deltaproc)
	println("Now running on ",nprocs()," processors.")
end=#

# sphere plot
n = 100 # number of points to plot=#
u = linspace(0., π, n) # the two angles that define the location of the points on the sphere
v = linspace(0., 2π, n)' # notice that I transpose this one with '
u .+= 0*v # the transposition is "carried" and expanded through the .+ operator to make u a n by n matrix
v .+= 0*u # same here
x = lens.c.x + lens.r*cos(u).*sin(v) # then standard vector math
y = lens.c.y + lens.r*sin(u).*sin(v)
z = lens.c.z + lens.r*cos(v)
plot_surface(x, y, z, color="blue",alpha=.5,linewidth=0,rstride=4,cstride=4) # and finally plotting the sphere with some transparency so that we'll be able to see the ray

# plotting a part of the retina, like a hemisphere, but not hemi... This code is exactly the same as for the initiation of the rays, when they are gonna hit the lens
θ = .5 # this is the (arbitrary) angle of the sphere-retina that will be plotted

# z and h are uniformly sampled from these ranges to generate the uniformly spaced points on the sphere
z = linspace(cos(θ),1.,n)
h = linspace(0.,2π,n)'

# now come the x, y, and z that lie on the retina's surface
x = retina.c.x + retina.r*sqrt(1. - z.^2).*cos(h)
y = retina.c.y + retina.r*sqrt(1. - z.^2).*sin(h)
z = retina.c.z - retina.r*z .+ 0.*h # I flip the z so it'll point downwards

plot_surface(x, y, z, color="green",alpha=.25,linewidth=0,rstride=4,cstride=4) # and finally plotting the sphere with some transparency

σ = mean(Δ) # mean error

max_stars_repo_name: JuliaTagBot/borat

max_issues_repo_name: JuliaTagBot/borat

max_forks_repo_name: JuliaTagBot/borat


text:
import Splittings:compute_interpolants, interpolate, CubicSpline

@testset "CubicSpline interpolation" begin

    function interpolation_test(n1 = 128)

        x1min, x1max = 0.0, 1.0
        x = collect(range(x1min,stop=x1max,length=n1))
        y = sin.(2π*x)
        x_new = zeros(Float64,n1)
        y_new = zeros(Float64,n1)
        coeffs = compute_interpolants(n1, y)
        for (i, xi) in enumerate(x)
            x_new[i] = xi - 0.1
            x_new[i] = x1min + mod(x_new[i]-x1min,x1max-x1min)
            y_new[i] = interpolate(coeffs, n1, x1min, x1max, x_new[i])
        end
        maximum(abs.(sin.(2π*(x.-0.1)) - y_new))

    end

    @test ≈(interpolation_test(), 0.0, atol=1e-7)

end

import LinearAlgebra: transpose
import Splittings: UniformMesh, advection!

@testset " BSL advections with cubic splines " begin

  n1, n2 = 128, 128
  x1min, x1max = -5, 10
  x2min, x2max = -5, 10
  mesh1 = UniformMesh(x1min, x1max, n1)
  mesh2 = UniformMesh(x2min, x2max, n2)

  f  = zeros(Float64,(n1,n2))
  f .= exp.(-mesh1.points.^2) * transpose(exp.(-mesh2.points.^2))
  fᵗ = zeros(Float64,(n2,n1))

  dt =  0.5

  e = ones(Float64, n1)
  v = ones(Float64, n2)

  advection!(f, mesh1,  v, dt, CubicSpline(), 1)
  advection!(f, mesh2,  e, dt, CubicSpline(), 2)
  advection!(f, mesh1, -v, dt, CubicSpline(), 1)
  advection!(f, mesh2, -e, dt, CubicSpline(), 2)

  f0 =  exp.(-mesh1.points.^2) * transpose(exp.(-mesh2.points.^2))
  println( maximum( abs.(f .- f0)))

  @test f ≈ f0 atol=1e-3

end


text:
for SHT in [SH, GSH]
    precompile(cache, (SHT, Type{Float64}, Int))
    precompile(biposh_flippoints, (SHT, Float64, Float64, Float64, Float64, Int, Int, Int, Int))
    precompile(biposh_flippoints, (SHT, Float64, Float64, Float64, Float64, Colon, Colon, Int, Int))
    precompile(biposh_flippoints, (SHT, Float64, Float64, Float64, Float64, Int, Int, L2L1Triangle))
    precompile(biposh_flippoints, (SHT, Float64, Float64, Float64, Float64, Colon, Colon, L2L1Triangle))
end
for YT in [PB, Hansen, Irreducible], B in [Cartesian, Polar, SphericalCovariant, HelicityCovariant]
    precompile(cache, (VSH{YT,B}, Type{Float64}, Int))
    precompile(biposh_flippoints, (VSH{YT,B}, Float64, Float64, Float64, Float64, Int, Int, Int, Int))
    precompile(biposh_flippoints, (VSH{YT,B}, Float64, Float64, Float64, Float64, Colon, Colon, Int, Int))
    precompile(biposh_flippoints, (VSH{YT,B}, Float64, Float64, Float64, Float64, Int, Int, L2L1Triangle))
    precompile(biposh_flippoints, (VSH{YT,B}, Float64, Float64, Float64, Float64, Colon, Colon, L2L1Triangle))
end


text:
using Genie
using Genie.Router
import Genie.Renderer.Json: json
using Genie.Renderer.Json, Genie.Requests

using Flux
using Flux: Chain, Dense, Conv, MaxPool, flatten, relu, cpu
using BSON
using Flux: onehotbatch, onecold, normalise, Dropout
using DataStructures
using StatsBase
#using GLMakie
#using Images

#img = load("./../datasetspectro2/down#3ae5c04f_nohash_0.wav.png")
#fig = Figure(size = (1000, 700), title = "Object Tracker")
#ax = GLMakie.Axis(
   # fig[1, 1],aspect = DataAspect())
#node = Node(rotr90(img))
#makieimg = image!(ax, node)

#display(fig)

function ConvNet5(; imgsize = (28, 28, 1), nclasses = 4)
    return Chain(
        normalise,
        Conv((3, 3), imgsize[end] => 32, relu),
        Conv((3, 3), 32 => 64, relu),
        MaxPool((2, 2)),
        Dropout(0.25),
        flatten,
        Dense(9216, 120, relu),
        Dropout(0.5),
        Dense(120, nclasses),
    )
end

predictdict = Dict(
    0 => "left",
    1 => "right",
    2 => "up",
    3 => "down",
)

device = cpu
model = ConvNet5() |> device

BSON.@load "./src/model_4.0.bson" model
queueofdata = Queue{Any}();
xtrain = zeros(Float32, (28, 28, 1, 15))
function launchServer(port)

    Genie.config.run_as_server = true
    Genie.config.server_host = "0.0.0.0"
    Genie.config.server_port = port

    println("port set to $(port)")

    #route("/") do
    #    "Hi there! This is server 1"
    #end

    route("/adddata", method = POST) do
        data1 = eval(Meta.parse(rawpayload()))
        enqueue!(queueofdata, data1);
        data1 = nothing
        #println("Current queue after post: ",length(queueofdata))
    end

    route("/digitreg", method = POST) do
        data = eval(Meta.parse(rawpayload()))
        data = reshape(data, (28, 28, 1))
        xtrain = zeros(Float32, (28, 28, 1, 2))
        xtrain[:, :, :, 1] = data
        results = model(xtrain)
        final = predictdict[findmax(results[:, 1])[2]-1]
        return final
    end

    route("/digitreg", method = GET) do
    	#println("Current queue: ",length(queueofdata))
    	if length(queueofdata) >= 1
	    	#node[] = rotr90(RGB{N0f8}.(dequeue!(queueofdata)))
	    	l = length(queueofdata)
		map(x-> xtrain[:,:,:,x] =  reshape(Float32.(dequeue!(queueofdata)), (28, 28, 1)), 1:l)

		results = model(xtrain)
		results = map(x-> predictdict[findmax(results[:, x])[2]-1],1:l)
		println("This here:",results)
		dictfinal = countmap(results)
		a = findmax(dictfinal)
		#if dictfinal[a[2]] == 3
			#return	"Idk what"
		#end
		# println(a)
		return a[2]
	end
    end

    Genie.AppServer.startup(async = false)
end

# launchServer(parse(Int, ARGS[1]))
launchServer(8003)

max_stars_repo_name: ashwani-rathee/Sardor

max_issues_repo_name: ashwani-rathee/Sardor

max_forks_repo_name: ashwani-rathee/Sardor


text:
using SparseDiffTools
using LightGraphs
using LightGraphs: SimpleGraph
using Test

using Random
Random.seed!(100)

# println("Starting acyclic coloring test...")
#= Test data =#
test_graphs = Vector{SimpleGraph}(undef, 0)
test_graphs_dir = Vector{SimpleDiGraph}(undef, 0)

for _ in 1:5
    nv = rand(5:20)
    ne = rand(1:100)
    graph = SimpleGraph(nv)
    for e in 1:ne
        v1 = rand(1:nv)
        v2 = rand(1:nv)
        while v1 == v2
            v2 = rand(1:nv)
        end
        add_edge!(graph, v1, v2)
    end
    push!(test_graphs, copy(graph))
end

#=
 Coloring needs to satisfy two conditions:

1. every pair of adjacent vertices receives distinct colors
(a distance-1 coloring)

2. For any vertex v, any color that leads to a two-colored path
involving v and three other vertices  is  impermissible  for  v.
In other words, every path on four vertices uses at least three
colors.
=#

#Sample graph from Gebremedhin AH, Manne F, Pothen A. **New Acyclic and Star Coloring Algorithms with Application to Computing Hessians**

#=

(1) ----- (2) ----- (3) ---- (4)
 | \      /          |       |
 |   \   /           |       |
 |     \/            |       |
 (5)---(6) -------- (7) ----(8)
 |      |         /  |       |
 |      |     (11)   |       |
 |      |    /       |       |
 (9) ---(10) ----- (12)-----(13)
  \                 /
   \              /
     ------------
=#

gx = SimpleGraph(13)

add_edge!(gx,1,2)
add_edge!(gx,1,5)
add_edge!(gx,1,6)
add_edge!(gx,2,3)
add_edge!(gx,2,6)
add_edge!(gx,3,4)
add_edge!(gx,3,7)
add_edge!(gx,4,8)
add_edge!(gx,5,6)
add_edge!(gx,5,9)
add_edge!(gx,6,7)
add_edge!(gx,6,10)
add_edge!(gx,7,8)
add_edge!(gx,7,11)
add_edge!(gx,7,12)
add_edge!(gx,8,13)
add_edge!(gx,9,10)
add_edge!(gx,9,12)
add_edge!(gx,10,11)
add_edge!(gx,10,12)
add_edge!(gx,12,13)

push!(test_graphs, gx)

#create directed copies
for g in test_graphs
    dg = SimpleDiGraph(nv(g))
    for e in edges(g)
        src1 = src(e)
        dst1 = dst(e)
        add_edge!(dg, src1, dst1)
        add_edge!(dg, dst1, src1)
    end
    push!(test_graphs_dir, dg)
end

for i in 1:5
    g = test_graphs[i]
    dg = test_graphs_dir[i]

    out_colors = SparseDiffTools.color_graph(g, SparseDiffTools.AcyclicColoring())

    #test condition 1
    for v in vertices(g)
        color = out_colors[v]
        for j in inneighbors(g, v)
            @test out_colors[j] != color
        end
    end
end

for i in 3:4
    g = test_graphs[i]
    dg = test_graphs_dir[i]

    out_colors = SparseDiffTools.color_graph(g, SparseDiffTools.AcyclicColoring())

    #test condition 2
    cycles = simplecycles(dg)
    for c in cycles
        colors = zeros(Int, 0)
        if length(c) > 2
            for v in c
                push!(colors, out_colors[v])
            end
            @test length(unique(colors)) >= 3
        end
    end
    # println("finished testing graph $i")
end

# println("finished testing...")


text:
# This file contains code that was formerly a part of Julia. License is MIT: http://julialang.org/license

# PriorityQueue
# -------------

"""
    PriorityQueue{K, V}([ord])

Construct a new [`PriorityQueue`](@ref), with keys of type
`K` and values/priorites of type `V`.
If an order is not given, the priority queue is min-ordered using
the default comparison for `V`.

A `PriorityQueue` acts like a `Dict`, mapping values to their
priorities, with the addition of a `dequeue!` function to remove the
lowest priority element.

```jldoctest
julia> a = PriorityQueue(["a","b","c"],[2,3,1],Base.Order.Forward)
PriorityQueue{String,Int64,Base.Order.ForwardOrdering} with 3 entries:
  "c" => 1
  "b" => 3
  "a" => 2
```
"""
mutable struct PriorityQueue{K,V,O<:Ordering} <: AbstractDict{K,V}
    # Binary heap of (element, priority) pairs.
    xs::Array{Pair{K,V}, 1}
    o::O

    # Map elements to their index in xs
    index::Dict{K, Int}

    function PriorityQueue{K,V,O}(o::O) where {K,V,O<:Ordering}
        new{K,V,O}(Vector{Pair{K,V}}(), o, Dict{K, Int}())
    end

    function PriorityQueue{K,V,O}(o::O, itr) where {K,V,O<:Ordering}
        xs = Vector{Pair{K,V}}(undef, length(itr))
        index = Dict{K, Int}()
        for (i, (k, v)) in enumerate(itr)
            xs[i] = Pair{K,V}(k, v)
            if haskey(index, k)
                throw(ArgumentError("PriorityQueue keys must be unique"))
            end
            index[k] = i
        end
        pq = new{K,V,O}(xs, o, index)

        # heapify
        for i in heapparent(length(pq.xs)):-1:1
            percolate_down!(pq, i)
        end

        pq
    end
end

# Any-Any constructors
PriorityQueue(o::Ordering=Forward) = PriorityQueue{Any,Any,typeof(o)}(o)

# Construction from Pairs
PriorityQueue(ps::Pair...) = PriorityQueue(Forward, ps)
PriorityQueue(o::Ordering, ps::Pair...) = PriorityQueue(o, ps)
PriorityQueue{K,V}(ps::Pair...) where {K,V} = PriorityQueue{K,V,ForwardOrdering}(Forward, ps)
PriorityQueue{K,V}(o::Ord, ps::Pair...) where {K,V,Ord<:Ordering} = PriorityQueue{K,V,Ord}(o, ps)

# Construction specifying Key/Value types
# e.g., PriorityQueue{Int,Float64}([1=>1, 2=>2.0])
PriorityQueue{K,V}(kv) where {K,V} = PriorityQueue{K,V}(Forward, kv)
function PriorityQueue{K,V}(o::Ord, kv) where {K,V,Ord<:Ordering}
    try
        PriorityQueue{K,V,Ord}(o, kv)
    catch e
        if not_iterator_of_pairs(kv)
            throw(ArgumentError("PriorityQueue(kv): kv needs to be an iterator of tuples or pairs"))
        else
            rethrow(e)
        end
    end
end

# Construction inferring Key/Value types from input
# e.g. PriorityQueue{}

PriorityQueue(o1::Ordering, o2::Ordering) = throw(ArgumentError("PriorityQueue with two parameters must be called with an Ordering and an interable of pairs"))
PriorityQueue(kv, o::Ordering=Forward) = PriorityQueue(o, kv)
function PriorityQueue(o::Ordering, kv)
    try
        _priority_queue_with_eltype(o, kv, eltype(kv))
    catch e
        if not_iterator_of_pairs(kv)
            throw(ArgumentError("PriorityQueue(kv): kv needs to be an iterator of tuples or pairs"))
        else
            rethrow(e)
        end
    end
end

_priority_queue_with_eltype(o::Ord, ps, ::Type{Pair{K,V}} ) where {K,V,Ord} = PriorityQueue{  K,  V,Ord}(o, ps)
_priority_queue_with_eltype(o::Ord, kv, ::Type{Tuple{K,V}}) where {K,V,Ord} = PriorityQueue{  K,  V,Ord}(o, kv)
_priority_queue_with_eltype(o::Ord, ps, ::Type{Pair{K}}   ) where {K,  Ord} = PriorityQueue{  K,Any,Ord}(o, ps)
_priority_queue_with_eltype(o::Ord, kv, ::Type            ) where {    Ord} = PriorityQueue{Any,Any,Ord}(o, kv)

## TODO: It seems impossible (or at least very challenging) to create the eltype below.
##       If deemed possible, please create a test and uncomment this definition.
# _priority_queue_with_eltype{  D,Ord}(o::Ord, ps, ::Type{Pair{K,V} where K}) = PriorityQueue{Any,  D,Ord}(o, ps)

length(pq::PriorityQueue) = length(pq.xs)
isempty(pq::PriorityQueue) = isempty(pq.xs)
haskey(pq::PriorityQueue, key) = haskey(pq.index, key)

"""
    peek(pq)

Return the lowest priority key from a priority queue without removing that
key from the queue.
"""
peek(pq::PriorityQueue) = pq.xs[1]

function percolate_down!(pq::PriorityQueue, i::Integer)
    x = pq.xs[i]
    @inbounds while (l = heapleft(i)) <= length(pq)
        r = heapright(i)
        j = r > length(pq) || lt(pq.o, pq.xs[l].second, pq.xs[r].second) ? l : r
        if lt(pq.o, pq.xs[j].second, x.second)
            pq.index[pq.xs[j].first] = i
            pq.xs[i] = pq.xs[j]
            i = j
        else
            break
        end
    end
    pq.index[x.first] = i
    pq.xs[i] = x
end

function percolate_up!(pq::PriorityQueue, i::Integer)
    x = pq.xs[i]
    @inbounds while i > 1
        j = heapparent(i)
        if lt(pq.o, x.second, pq.xs[j].second)
            pq.index[pq.xs[j].first] = i
            pq.xs[i] = pq.xs[j]
            i = j
        else
            break
        end
    end
    pq.index[x.first] = i
    pq.xs[i] = x
end

# Equivalent to percolate_up! with an element having lower priority than any other
function force_up!(pq::PriorityQueue, i::Integer)
    x = pq.xs[i]
    @inbounds while i > 1
        j = heapparent(i)
        pq.index[pq.xs[j].first] = i
        pq.xs[i] = pq.xs[j]
        i = j
    end
    pq.index[x.first] = i
    pq.xs[i] = x
end

function getindex(pq::PriorityQueue{K,V}, key) where {K,V}
    pq.xs[pq.index[key]].second
end

function get(pq::PriorityQueue{K,V}, key, deflt) where {K,V}
    i = get(pq.index, key, 0)
    i == 0 ? deflt : pq.xs[i].second
end

# Change the priority of an existing element, or equeue it if it isn't present.
function setindex!(pq::PriorityQueue{K, V}, value, key) where {K,V}
    if haskey(pq, key)
        i = pq.index[key]
        oldvalue = pq.xs[i].second
        pq.xs[i] = Pair{K,V}(key, value)
        if lt(pq.o, oldvalue, value)
            percolate_down!(pq, i)
        else
            percolate_up!(pq, i)
        end
    else
        enqueue!(pq, key, value)
    end
    value
end

"""
    enqueue!(pq, k=>v)

Insert the a key `k` into a priority queue `pq` with priority `v`.

```jldoctest
julia> a = PriorityQueue(PriorityQueue("a"=>1, "b"=>2, "c"=>3))
PriorityQueue{String,Int64,Base.Order.ForwardOrdering} with 3 entries:
  "c" => 3
  "b" => 2
  "a" => 1

julia> enqueue!(a, "d"=>4)
PriorityQueue{String,Int64,Base.Order.ForwardOrdering} with 4 entries:
  "c" => 3
  "b" => 2
  "a" => 1
  "d" => 4
```
"""
function enqueue!(pq::PriorityQueue{K,V}, pair::Pair{K,V}) where {K,V}
    key = pair.first
    if haskey(pq, key)
        throw(ArgumentError("PriorityQueue keys must be unique"))
    end
    push!(pq.xs, pair)
    pq.index[key] = length(pq)
    percolate_up!(pq, length(pq))

    return pq
end

"""
enqueue!(pq, k, v)

Insert the a key `k` into a priority queue `pq` with priority `v`.

"""
enqueue!(pq::PriorityQueue, key, value) = enqueue!(pq, key=>value)
enqueue!(pq::PriorityQueue{K,V}, kv) where {K,V} = enqueue!(pq, Pair{K,V}(kv.first, kv.second))

"""
    dequeue!(pq)

Remove and return the lowest priority key from a priority queue.

```jldoctest
julia> a = PriorityQueue(["a","b","c"],[2,3,1],Base.Order.Forward)
PriorityQueue{String,Int64,Base.Order.ForwardOrdering} with 3 entries:
  "c" => 1
  "b" => 3
  "a" => 2

julia> dequeue!(a)
"c"

julia> a
PriorityQueue{String,Int64,Base.Order.ForwardOrdering} with 2 entries:
  "b" => 3
  "a" => 2
```
"""
function dequeue!(pq::PriorityQueue)
    x = pq.xs[1]
    y = pop!(pq.xs)
    if !isempty(pq)
        pq.xs[1] = y
        pq.index[y.first] = 1
        percolate_down!(pq, 1)
    end
    delete!(pq.index, x.first)
    x.first
end

function dequeue!(pq::PriorityQueue, key)
    idx = pq.index[key]
    force_up!(pq, idx)
    dequeue!(pq)
    key
end

"""
    dequeue_pair!(pq)

Remove and return a the lowest priority key and value from a priority queue as a pair.

```jldoctest
julia> a = PriorityQueue(["a","b","c"],[2,3,1],Base.Order.Forward)
PriorityQueue{String,Int64,Base.Order.ForwardOrdering} with 3 entries:
  "c" => 1
  "b" => 3
  "a" => 2

julia> dequeue_pair!(a)
"c" => 1

julia> a
PriorityQueue{String,Int64,Base.Order.ForwardOrdering} with 2 entries:
  "b" => 3
  "a" => 2
```
"""
function dequeue_pair!(pq::PriorityQueue)
    x = pq.xs[1]
    y = pop!(pq.xs)
    if !isempty(pq)
        pq.xs[1] = y
        pq.index[y.first] = 1
        percolate_down!(pq, 1)
    end
    delete!(pq.index, x.first)
    x
end

function dequeue_pair!(pq::PriorityQueue, key)
    idx = pq.index[key]
    force_up!(pq, idx)
    dequeue_pair!(pq)
end

"""
    delete!(pq, key)
Delete the mapping for the given key in a priority queue, and return the priority queue.
# Examples
```jldoctest
julia> q = PriorityQueue(Base.Order.Forward, "a"=>2, "b"=>3, "c"=>1)
PriorityQueue{String,Int64,Base.Order.ForwardOrdering} with 3 entries:
  "c" => 1
  "b" => 3
  "a" => 2
julia> delete!(q, "b")
DataStructures.PriorityQueue{String,Int64,Base.Order.ForwardOrdering} with 2 entries:
  "c" => 1
  "a" => 2
```
"""
function delete!(pq::PriorityQueue, key)
    dequeue_pair!(pq, key)
    pq
end

# Unordered iteration through key value pairs in a PriorityQueue
function _iterate(pq::PriorityQueue, state)
    state == nothing && return nothing
    (k, idx), i = state
    return (pq.xs[idx], i)
end

iterate(pq::PriorityQueue) = _iterate(pq, iterate(pq.index))

iterate(pq::PriorityQueue, i) = _iterate(pq, iterate(pq.index, i))


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

mutable struct IoK8sApiAppsV1beta1StatefulSetSpec <: SwaggerModel
    podManagementPolicy::Any # spec type: Union{ Nothing, String } # spec name: podManagementPolicy
    replicas::Any # spec type: Union{ Nothing, Int32 } # spec name: replicas
    revisionHistoryLimit::Any # spec type: Union{ Nothing, Int32 } # spec name: revisionHistoryLimit
    selector::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1LabelSelector } # spec name: selector
    serviceName::Any # spec type: Union{ Nothing, String } # spec name: serviceName
    template::Any # spec type: Union{ Nothing, IoK8sApiCoreV1PodTemplateSpec } # spec name: template
    updateStrategy::Any # spec type: Union{ Nothing, IoK8sApiAppsV1beta1StatefulSetUpdateStrategy } # spec name: updateStrategy
    volumeClaimTemplates::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1PersistentVolumeClaim} } # spec name: volumeClaimTemplates

    function IoK8sApiAppsV1beta1StatefulSetSpec(;podManagementPolicy=nothing, replicas=nothing, revisionHistoryLimit=nothing, selector=nothing, serviceName=nothing, template=nothing, updateStrategy=nothing, volumeClaimTemplates=nothing)
        o = new()
        validate_property(IoK8sApiAppsV1beta1StatefulSetSpec, Symbol("podManagementPolicy"), podManagementPolicy)
        setfield!(o, Symbol("podManagementPolicy"), podManagementPolicy)
        validate_property(IoK8sApiAppsV1beta1StatefulSetSpec, Symbol("replicas"), replicas)
        setfield!(o, Symbol("replicas"), replicas)
        validate_property(IoK8sApiAppsV1beta1StatefulSetSpec, Symbol("revisionHistoryLimit"), revisionHistoryLimit)
        setfield!(o, Symbol("revisionHistoryLimit"), revisionHistoryLimit)
        validate_property(IoK8sApiAppsV1beta1StatefulSetSpec, Symbol("selector"), selector)
        setfield!(o, Symbol("selector"), selector)
        validate_property(IoK8sApiAppsV1beta1StatefulSetSpec, Symbol("serviceName"), serviceName)
        setfield!(o, Symbol("serviceName"), serviceName)
        validate_property(IoK8sApiAppsV1beta1StatefulSetSpec, Symbol("template"), template)
        setfield!(o, Symbol("template"), template)
        validate_property(IoK8sApiAppsV1beta1StatefulSetSpec, Symbol("updateStrategy"), updateStrategy)
        setfield!(o, Symbol("updateStrategy"), updateStrategy)
        validate_property(IoK8sApiAppsV1beta1StatefulSetSpec, Symbol("volumeClaimTemplates"), volumeClaimTemplates)
        setfield!(o, Symbol("volumeClaimTemplates"), volumeClaimTemplates)
        o
    end
end # type IoK8sApiAppsV1beta1StatefulSetSpec

const _property_map_IoK8sApiAppsV1beta1StatefulSetSpec = Dict{Symbol,Symbol}(Symbol("podManagementPolicy")=>Symbol("podManagementPolicy"), Symbol("replicas")=>Symbol("replicas"), Symbol("revisionHistoryLimit")=>Symbol("revisionHistoryLimit"), Symbol("selector")=>Symbol("selector"), Symbol("serviceName")=>Symbol("serviceName"), Symbol("template")=>Symbol("template"), Symbol("updateStrategy")=>Symbol("updateStrategy"), Symbol("volumeClaimTemplates")=>Symbol("volumeClaimTemplates"))
const _property_types_IoK8sApiAppsV1beta1StatefulSetSpec = Dict{Symbol,String}(Symbol("podManagementPolicy")=>"String", Symbol("replicas")=>"Int32", Symbol("revisionHistoryLimit")=>"Int32", Symbol("selector")=>"IoK8sApimachineryPkgApisMetaV1LabelSelector", Symbol("serviceName")=>"String", Symbol("template")=>"IoK8sApiCoreV1PodTemplateSpec", Symbol("updateStrategy")=>"IoK8sApiAppsV1beta1StatefulSetUpdateStrategy", Symbol("volumeClaimTemplates")=>"Vector{IoK8sApiCoreV1PersistentVolumeClaim}")
Base.propertynames(::Type{ IoK8sApiAppsV1beta1StatefulSetSpec }) = collect(keys(_property_map_IoK8sApiAppsV1beta1StatefulSetSpec))
Swagger.property_type(::Type{ IoK8sApiAppsV1beta1StatefulSetSpec }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiAppsV1beta1StatefulSetSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiAppsV1beta1StatefulSetSpec }, property_name::Symbol) =  _property_map_IoK8sApiAppsV1beta1StatefulSetSpec[property_name]

function check_required(o::IoK8sApiAppsV1beta1StatefulSetSpec)
    (getproperty(o, Symbol("serviceName")) === nothing) && (return false)
    (getproperty(o, Symbol("template")) === nothing) && (return false)
    true
end

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


text:
using HighPrecisionNonsymmetricEigenproblem
using Test

@testset "HighPrecisionNonsymmetricEigenproblem.jl" begin

end


text:
#=
reverser.jl
based on
https://stackoverflow.com/questions/27411401/julia-reverse-n-dimensional-arrays
=#

using MIRT: reverser

using Test: @test

@test reverser(1:3) == 3:-1:1
@test reverser(1:3, 1) == 3:-1:1
@test reverser((1:3)', 1) == (1:3)'
@test reverser((1:3)', 2) == (3:-1:1)'


text:
@symbol_func function P_kappa(cur_reactor::AbstractReactor)
  cur_P = K_kappa(cur_reactor)

  cur_P *= cur_reactor.R_0 ^ 3

  cur_P *= cur_reactor.n_bar

  cur_P *= cur_reactor.T_bar

  cur_P /= tau_E(cur_reactor)

  cur_P
end


text:
@testset "Plastic" begin
    # constructor
    m = Plastic(E=200e3, ν=0.3, σ_y=200., H=50., r=0.5, κ_∞=13., α_∞=13.)
    cache = get_cache(m)

    # initial state
    state = initial_material_state(m)

    # strain at yield point for uniaxial stress
    ε11_yield = m.σ_y / m.E

    # elastic branch
    Δε = SymmetricTensor{2,3,Float64}((i,j) -> i==1 && j==1 ? 0.5ε11_yield : (i == j ? -0.5ε11_yield*m.ν : 0.0))
    σ, ∂σ∂ε, temp_state = material_response(m, Δε, state; cache=cache)
    @test σ ≈ SymmetricTensor{2,3,Float64}((i,j) -> i==1 && j==1 ? 0.5ε11_yield*m.E : 0.0)
    @test ∂σ∂ε == m.Eᵉ
    @test temp_state.κ == 0.0
    @test temp_state.α == zero(SymmetricTensor{2,3})
    @test temp_state.μ == 0.0

    # yield point
    Δε = Δε = SymmetricTensor{2,3,Float64}((i,j) -> i==1 && j==1 ? ε11_yield : (i == j ? -ε11_yield*m.ν : 0.0))
    σ, ∂σ∂ε, temp_state = material_response(m, Δε, state; cache=cache)
    @test sqrt(3/2)*norm(dev(σ-temp_state.α)) ≈ m.σ_y

    # plastic branch
    Δε = Δε = SymmetricTensor{2,3,Float64}((i,j) -> i==1 && j==1 ? 2ε11_yield : (i == j ? -2ε11_yield*m.ν : 0.0))
    σ, ∂σ∂ε, temp_state = material_response(m, Δε, state; cache=cache)
    @test sqrt(3/2)*norm(dev(σ-temp_state.α)) > m.σ_y
end

function get_Plastic_loading()
    strain1 = range(0.0,  0.005, length=5)
    strain2 = range(0.005, 0.001, length=5)
    strain3 = range(0.001, 0.007, length=5)

    _C = [strain1[2:end]..., strain2[2:end]..., strain3[2:end]...]
    ε = [SymmetricTensor{2,3}((x, x/10, 0.0, 0.0, 0.0, 0.0)) for x in _C]

    return ε
end

@testset "Plastic jld2" begin
    m = Plastic(E=200e3, ν=0.3, σ_y=200., H=50., r=0.5, κ_∞=13., α_∞=13.)

    loading = get_Plastic_loading()
    check_jld2(m, loading, "Plastic1")#, debug_print=true, OVERWRITE_JLD2=true)
end


text:
#=
utilities for HSL
only works with Float64 and Float32
=#

import LinearAlgebra.BlasReal

mutable struct HSLSymCache{T <: BlasReal} <: SparseSymCache{T}
    analyzed::Bool
    ma57::HSL.Ma57{T}
    function HSLSymCache{T}() where {T <: BlasReal}
        cache = new{T}()
        cache.analyzed = false
        return cache
    end
end

int_type(::HSLSymCache) = Int

function update_fact(
    cache::HSLSymCache{T},
    A::SparseMatrixCSC{T, Int}
    ) where {T <: BlasReal}
    if !cache.analyzed
        cache.ma57 = HSL.Ma57(A)
        cache.analyzed = true
    else
        copyto!(cache.ma57.vals, A.nzval)
    end
    HSL.ma57_factorize(cache.ma57)
    return
end

function inv_prod(
    cache::HSLSymCache{T},
    x::Vector{T},
    A::SparseMatrixCSC{T, Int},
    b::Vector{T},
    ) where {T <: BlasReal}
    # MA57 only has the option to take iterative refinement steps for a single-column RHS
    ma57 = cache.ma57
    copyto!(x, b)
    HSL.ma57_solve!(ma57, x)

    # info we could print
    # ma57.info.backward_error1
    # ma57.info.backward_error2
    # ma57.info.matrix_inf_norm
    # ma57.info.scaled_residuals
    # ma57.info.cond1
    # ma57.info.cond2
    # ma57.info.error_inf_norm

    return x
end

free_memory(cache::HSLSymCache) = nothing


text:
__precompile__()

module RoamesGeometry

using StaticArrays
using CoordinateTransformations
using Rotations

import Base: getindex, setindex!, size, length, lastindex, in, isempty, eltype, IndexStyle,
             push!, pop!, show
using Base: @propagate_inbounds, RefValue
import CoordinateTransformations: Transformation, IdentityTransformation,
                                  Translation, LinearMap, AffineMap

if VERSION < v"0.7"
    import Base: endof
else
    using LinearAlgebra, Pkg

    import Base: lastindex
    linspace(a,b,c) = range(a, stop = b, length = c)
end

export AbstractRegion, BoundingBox, TriangularPrism, Sphere, Triangle, Circle, Cylinder,
       Line, LineString, Quadratic, Catenary, Polygon

export boundingbox, pad, intersects, wireframe, distance, powerline_distances, volume,
       area, might_intersect, isclosed, issimple, closest_point, closest_points, convert2d,
       convert3d, winding_number, containsMoreThanN

export wkt, load_wkt, save_wkt, read_wkt, write_wkt

# Any affine transformation type
const AffineTransformation = Union{IdentityTransformation, Translation, LinearMap, AffineMap}

# Zero-dimensional entities
#  - Points are AbstractVector{<:Real} of length 2 or 3, usually a StaticVector{N, <:Real}
convert2d(p::StaticVector{2, <:Real}) = p
convert2d(p::StaticVector{3, <:Real}) = similar_type(typeof(p), Size(2))((p[1], p[2]))
convert3d(p::StaticVector{2, T}, z::Real = zero(T)) where {T<:Real} = similar_type(typeof(p), Size(3))((p[1], p[2], T(z)))
convert3d(p::StaticVector{3, T}, z::Real) where {T<:Real} = similar_type(typeof(p), Size(3))((p[1], p[2], T(z)))
convert3d(p::StaticVector{3, <:Real}) = p

# 1D geometries
include("Line.jl")
include("LineString.jl")
include("Quadratic.jl")
include("Catenary.jl")

# 2D geometries
# - triangle?
# - mesh?
# - polygon?

# 3D geometries
include("AbstractRegion.jl")
include("TriangularPrism.jl")
include("Cylinder.jl")
include("Sphere.jl")
include("Triangle.jl")
include("Circle.jl")
include("Polygon.jl")
include("BoundingBox.jl")

include("might_intersect.jl") # included last as uses unions of geometries
include("distances.jl") # calculates the distances between e.g. points and catenaries
include("wkt.jl") # read and write WKT files
include("displaz.jl") # optionally loads methods if the Displaz package is installed

if VERSION >= v"0.7.0"
    using AcceleratedArrays
    using TypedTables
    using FileIO
    using LasIO
    using HDF5
    using Colors
    using FixedPointNumbers
    using Dates

    export GridIndex
    export load_pointcloud, save_pointcloud

    include("GridIndex.jl")
    include("pointcloud_io.jl")
end

end # module


text:
export Localization, Locgaspari, periodicmetric!, periodicmetric, cartesianmetric, gaspari

struct Localization
    L::Float64
    Gxy::Function
    Gyy::Function
end

# function Locgaspari(d, L, periodic::Bool)
#
#     G = zeros(d,d)
#     @inbounds for i=1:d
#         # Check if the domain is periodic
#         if periodic == true
#             @inbounds for j=i:d
#                 rdom = min(abs(j - i), abs(-d + (j-1) - i), abs(d + j - i))
#                 G[i,j] = gaspari(rdom/L)
#             end
#         else
#             @inbounds for j=i:d
#             rdom = abs(j - i)
#             G[i,j] = gaspari(rdom/L)
#             end
#         end
#
#     end
#     return Symmetric(G)
# end

# Some metric for the distance between variables

periodicmetric!(i,j,d) = min(abs(j - i), abs(-d + j  - i), abs(d + j - i))
periodicmetric(d) = (i,j) -> periodicmetric!(i,j,d)

cartesianmetric(i,j) = abs(i-j)

# Construct a possibly non-square localisation matrix using
# the Gaspari-Cohn kernel
function Locgaspari(d::Tuple{Int64, Int64}, L, metric::Function)
    dx, dy = d
    G = zeros(dx,dy)
    @inbounds for i=1:dx
                 for j=1:dy
                 G[i,j] = gaspari(metric(i,j)/L)
                 end
               end
    return G
end

# Caspari-Cohn kernel, Formula found in Data assimilation in the geosciences:
# An overview of methods, issues, and perspectives
g1(r) = 1 - (5/3)*r^2 +(5/8)*r^3 +(1/2)*r^4 -0.25*r^5
g2(r) = 4 - 5r +(5/3)*r^2 + (5/8)*r^3 -(1/2)*r^4 +(1/12)*r^5 -(2/3)*r^(-1)
gaspari(r) = r>2.0 ? 0.0 : r<1.0 ? g1(r) : g2(r)