text:
#############################################################################
# MathProgBase.jl
#############################################################################
# test/conininterface.jl
# Test the MathProgBase.jl interface for a provided conic solver.
#############################################################################

using Compat.Test
using Compat.LinearAlgebra
using Compat.SparseArrays
using MathProgBase

function coniclineartest(solver::MathProgBase.AbstractMathProgSolver;duals=false, tol=1e-6)
    @testset "Testing linear problems through conic interface with $solver" begin
        # Problem LIN1 - all vars in nonneg cone
        # min -3x - 2y - 4z
        # st    x +  y +  z == 3 '4' -> z=2, x->2, obj -> -14
        #            y +  z == 2
        #       x>=0 y>=0 z>=0
        # Opt solution = -11
        # x = 1, y = 0, z = 2
        @testset "LIN1" begin
            c = [-3.0, -2.0, -4.0]
            A = [ 1.0   1.0   1.0;
                  0.0   1.0   1.0]
            b = [ 3.0,  2.0]
            m = MathProgBase.ConicModel(solver)
            MathProgBase.loadproblem!(m, c, A, b, [(:Zero,1:2)], [(:NonNeg, 1:3)])
            MathProgBase.optimize!(m)
            @test MathProgBase.status(m) == :Optimal
            @test isapprox(MathProgBase.getobjval(m), -11, atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[1], 1.0, atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[2], 0.0, atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[3], 2.0, atol=tol)
            if duals
                d = MathProgBase.getdual(m)
                @test isapprox(d[1], 3.0, atol=tol)
                @test isapprox(d[2], 1.0, atol=tol)
                vardual = c + A'd
                var = MathProgBase.getvardual(m)
                @test isapprox(vardual, var, atol=tol)
            end
        end

        # Problem LIN1A - same as Problem LIN1, but with variable bounds
        #              as constraints instead
        # min   -3x - 2y - 4z
        # st  [3] - [x + y + z] ZERO
        #     [2] - [    y + z] ZERO
        #     [0] - [x        ] NONPOS   x>=0 -> 0<=x -> 0-x<=0
        #     [0] - [    y    ] NONPOS   y>=0 -> 0<=y -> 0-y<=0
        #     [0] - [       -z] NONNEG   z>=0 -> 0<=z -> 0-z<=0 -> 0+z>=0
        # Opt solution = -11
        # x = 1, y = 0, z = 2
        @testset "LIN1A" begin
            c = [-3.0, -2.0, -4.0]
            A = [ 1.0   1.0   1.0;
                  0.0   1.0   1.0;
                  1.0   0.0   0.0;
                  0.0   1.0   0.0;
                  0.0   0.0  -1.0]
            b = [ 3.0,  2.0,  0.0,  0.0,  0.0]
            m = MathProgBase.ConicModel(solver)
            MathProgBase.loadproblem!(m,c, A, b, [(:Zero,1:2),(:NonPos,3:4),(:NonNeg,5)],
            [(:Free, 1:3)])
            MathProgBase.optimize!(m)
            @test MathProgBase.status(m) == :Optimal
            @test isapprox(MathProgBase.getobjval(m), -11, atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[1], 1.0, atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[2], 0.0, atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[3], 2.0, atol=tol)
            if duals
                d = MathProgBase.getdual(m)
                @test isapprox(d[1], 3.0, atol=tol)
                @test isapprox(d[2], 1.0, atol=tol)
                @test isapprox(d[3], 0.0, atol=tol)
                @test isapprox(d[4], -2.0, atol=tol)
                @test isapprox(d[5], 0.0, atol=tol)
                vardual = c + A'd
                var = MathProgBase.getvardual(m)
                @test isapprox(vardual, var, atol=tol)
            end
        end

        # Problem LIN2 - mixed free, nonneg, nonpos, zero, shuffled cones
        # min  3x + 2y - 4z + 0s
        # st    x           -  s  == -4    (i.e. x >= -4)
        #            y            == -3
        #       x      +  z       == 12
        #       x free
        #       y <= 0
        #       z >= 0
        #       s zero
        # Opt solution = -82
        # x = -4, y = -3, z = 16, s == 0
        @testset "LIN2" begin
            m = MathProgBase.ConicModel(solver)
            MathProgBase.loadproblem!(m,
            [ 3.0,  2.0, -4.0,  0.0],
            [ 1.0   0.0   0.0  -1.0;
              0.0   1.0   0.0   0.0;
              1.0   0.0   1.0   0.0],
            [-4.0, -3.0, 12.0],
            [(:Zero,1:3)],
            [(:Free,1), (:NonNeg,3), (:Zero,4), (:NonPos,2)])
            MathProgBase.optimize!(m)
            @test MathProgBase.status(m) == :Optimal
            @test isapprox(MathProgBase.getobjval(m), -82, atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[1], -4.0, atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[2], -3.0, atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[3], 16.0, atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[4], 0.0, atol=tol)
            if duals
                d = MathProgBase.getdual(m)
                @test isapprox(d[1], -7.0, atol=tol)
                @test isapprox(d[2], -2.0, atol=tol)
                @test isapprox(d[3], 4.0, atol=tol)
            end
        end

        # Problem LIN2A - Problem LIN2 but with y,z variable bounds as constraints
        # min       3x + 2y - 4z + 0s
        # st  [-4]-[ x           -  s] ZERO    (i.e. x >= -4)
        #     [-3]-[      y          ] ZERO
        #     [12]-[ x      +  z     ] ZERO
        #     [ 0]-[      y          ] NONEG
        #     [ 0]-[           z     ] NONPOS
        #       x free
        #       s zero
        # Opt solution = -82
        # x = -4, y = -3, z = 16, s == 0
        @testset "LIN2A" begin
            m = MathProgBase.ConicModel(solver)
            MathProgBase.loadproblem!(m,
            [ 3.0,  2.0, -4.0,  0.0],
            [ 1.0   0.0   0.0  -1.0;
              0.0   1.0   0.0   0.0;
              1.0   0.0   1.0   0.0;
              0.0   1.0   0.0   0.0;
              0.0   0.0   1.0   0.0],
            [-4.0, -3.0, 12.0,  0.0,  0.0],
            [(:Zero,1:3),(:NonNeg,4),(:NonPos,5)],
            [(:Free,1), (:NonNeg,3), (:Zero,4), (:NonPos,2)])
            MathProgBase.optimize!(m)
            @test MathProgBase.status(m) == :Optimal
            @test isapprox(MathProgBase.getobjval(m), -82, atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[1], -4.0, atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[2], -3.0, atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[3], 16.0, atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[4], 0.0, atol=tol)
            if duals
                d = MathProgBase.getdual(m)
                @test isapprox(d[1], -7.0, atol=tol)
                @test isapprox(d[2], -2.0, atol=tol)
                @test isapprox(d[3], 4.0, atol=tol)
                @test isapprox(d[4], 0.0, atol=tol)
                @test isapprox(d[5], 0.0, atol=tol)
            end
        end

        # Problem LIN3 - Infeasible LP
        # min  0
        # s.t. x ≥ 1
        #      x ≤ -1
        # in conic form:
        # min 0
        # s.t. -1 + x ∈ R₊
        #       1 + x ∈ R₋
        @testset "LIN3 infeasible" begin
            b = [-1, 1]
            A = reshape([-1, -1], 2, 1)
            c = [0]
            constr_cones = [(:NonNeg,1:1),(:NonPos,2:2)]
            var_cones = [(:Free,1:1)]

            m = MathProgBase.ConicModel(solver)
            MathProgBase.loadproblem!(m, c, A, b, constr_cones, var_cones)
            MathProgBase.optimize!(m)
            @test MathProgBase.status(m) == :Infeasible
            if duals
                y = MathProgBase.getdual(m)
                @test y[1] > 0
                @test y[2] < 0
                @test isapprox((A'y)[1], 0.0, atol=tol)
                @test -dot(b,y) > 0
            end
        end

        # Problem LIN4 - Infeasible LP
        # min  0
        # s.t. x ≥ 1
        #      x ≤ 0
        # in conic form:
        # min 0
        # s.t. -1 + x ∈ R₊
        #           x ∈ R₋
        @testset "LIN4 infeasible" begin
            b = [-1]
            A = reshape([-1], 1, 1)
            c = [0]
            constr_cones = [(:NonNeg,1:1)]
            var_cones = [(:NonPos,1:1)]

            m = MathProgBase.ConicModel(solver)
            MathProgBase.loadproblem!(m, c, A, b, constr_cones, var_cones)
            MathProgBase.optimize!(m)
            @test MathProgBase.status(m) == :Infeasible
            if duals
                y = MathProgBase.getdual(m)
                @test y[1] > 0
                @test -dot(b,y) > 0
            end
        end
    end
end

function conicSOCtest(solver::MathProgBase.AbstractMathProgSolver;duals=false, tol=1e-6)
    @testset "Testing SOC problems through conic interface with $solver" begin
        # Problem SOC1
        # min 0x - 1y - 1z
        #  st  x            == 1
        #      x >= ||(y,z)||
        @testset "SOC1" begin
            m = MathProgBase.ConicModel(solver)
            MathProgBase.loadproblem!(m,
            [ 0.0, -1.0, -1.0],
            [ 1.0   0.0   0.0],
            [ 1.0],
            [(:Zero,1)],
            [(:SOC,1:3)])
            MathProgBase.optimize!(m)
            @test MathProgBase.status(m) == :Optimal
            @test isapprox(MathProgBase.getobjval(m), -sqrt(2.0), atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[1], 1.0, atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[2], 1.0/sqrt(2.0), atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[3], 1.0/sqrt(2.0), atol=tol)
            if duals
                d = MathProgBase.getdual(m)
                @test isapprox(d[1], sqrt(2.0), atol=tol)
            end
        end

        # Problem SOC1A - Problem SOC1 but in ECOS form
        # min      0x - 1y - 1z
        #  st [1]-[ x          ] ZERO
        #     [0]-[-x          ] SOC
        #     [0]-[     -y     ] SOC
        #     [0]-[          -z] SOC
        @testset "SOC1A" begin
            m = MathProgBase.ConicModel(solver)
            MathProgBase.loadproblem!(m,
            [ 0.0, -1.0, -1.0],
            [ 1.0   0.0   0.0;
             -1.0   0.0   0.0;
              0.0  -1.0   0.0;
              0.0   0.0  -1.0],
            [ 1.0, 0.0, 0.0, 0.0],
            Any[(:Zero,1),(:SOC,2:4)],
            [(:Free,1:3)])
            MathProgBase.optimize!(m)
            @test MathProgBase.status(m) == :Optimal
            @test isapprox(MathProgBase.getobjval(m), -sqrt(2.0), atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[1], 1.0, atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[2], 1.0/sqrt(2.0), atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[3], 1.0/sqrt(2.0), atol=tol)
            if duals
                d = MathProgBase.getdual(m)
                @test isapprox(d[1], sqrt(2.0), atol=tol)
                @test isapprox(d[2], sqrt(2.0), atol=tol)
                @test isapprox(d[3], -1.0, atol=tol)
                @test isapprox(d[4], -1.0, atol=tol)
            end
            # Permute indices - used to cause Mosek failure
            m = MathProgBase.ConicModel(solver)
            MathProgBase.loadproblem!(m,
            [ 0.0, -1.0, -1.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],
            [ 1.0, 0.0, 0.0, 0.0],
            Any[(:Zero,1),(:SOC,[4,3,2])],
            [(:Free,1:3)])
            MathProgBase.optimize!(m)
            @test MathProgBase.status(m) == :Optimal
            @test isapprox(MathProgBase.getobjval(m), -sqrt(2.0), atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[1], 1.0, atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[2], 1.0/sqrt(2.0), atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[3], 1.0/sqrt(2.0), atol=tol)
            if duals
                d = MathProgBase.getdual(m)
                @test isapprox(d[1], sqrt(2.0), atol=tol)
                @test isapprox(d[4], sqrt(2.0), atol=tol)
                @test isapprox(d[3], -1.0, atol=tol)
                @test isapprox(d[2], -1.0, atol=tol)
            end

        end

        # Problem SOC2
        # min  x
        # s.t. y ≥ 1/√2
        #      x² + y² ≤ 1
        # in conic form:
        # min  x
        # s.t.  -1/√2 + y ∈ R₊
        #        1 - t ∈ {0}
        #      (t,x,y) ∈ SOC₃
        @testset "SOC2" begin
            b = [-1/sqrt(2), 1, 0, 0, 0]
            A = [ 0 -1 0
                  0 0 1
                  0 0 -1
                 -1 0 0
                  0 -1 0 ]
            c = [ 1, 0, 0 ]
            constr_cones = [(:NonNeg,1:1),(:Zero,2:2),(:SOC,3:5)]
            var_cones = [(:Free,1:3)]

            m = MathProgBase.ConicModel(solver)
            MathProgBase.loadproblem!(m, c, A, b, constr_cones, var_cones)
            MathProgBase.optimize!(m)
            @test MathProgBase.status(m) == :Optimal
            @test isapprox(MathProgBase.getobjval(m), -1/sqrt(2), atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[1], -1/sqrt(2), atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[2], 1/sqrt(2), atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[3], 1.0, atol=tol)

            if duals
                y = MathProgBase.getdual(m)
                @test isapprox(-dot(b,y), -1/sqrt(2), atol=tol)
                @test isapprox(norm(c + A'y), 0.0, atol=tol)
                s = MathProgBase.getvardual(m)
                @test isapprox(norm(s), 0.0, atol=tol)
            end
        end

        # Problem SOC2A
        # Same as above but with NonPos instead of NonNeg
        # min  x
        # s.t.  1/√2 - y ∈ R₋
        #        1 - t ∈ {0}
        #      (t,x,y) ∈ SOC₃
        @testset "SOC2A" begin
            b = [1/sqrt(2), 1, 0, 0, 0]
            A = [ 0 1 0
                  0 0 1
                  0 0 -1
                  -1 0 0
                  0 -1 0 ]
            c = [ 1, 0, 0 ]
            constr_cones = [(:NonPos,1:1),(:Zero,2:2),(:SOC,3:5)]
            var_cones = [(:Free,1:3)]

            m = MathProgBase.ConicModel(solver)
            MathProgBase.loadproblem!(m, c, A, b, constr_cones, var_cones)
            MathProgBase.optimize!(m)
            @test MathProgBase.status(m) == :Optimal
            @test isapprox(MathProgBase.getobjval(m), -1/sqrt(2), atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[1], -1/sqrt(2), atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[2], 1/sqrt(2), atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[3], 1.0, atol=tol)

            if duals
                y = MathProgBase.getdual(m)
                @test isapprox(-dot(b,y), -1/sqrt(2), atol=tol)
                @test isapprox(norm(c + A'y), 0.0, atol=tol)
                s = MathProgBase.getvardual(m)
                @test isapprox(norm(s), 0.0, atol=tol)
            end
        end

        # Problem SOC3 - Infeasible
        # min 0
        # s.t. y ≥ 2
        #      x ≤ 1
        #      |y| ≤ x
        # in conic form:
        # min 0
        # s.t. -2 + y ∈ R₊
        #      -1 + x ∈ R₋
        #       (x,y) ∈ SOC₂
        @testset "SOC3 infeasible" begin
            b = [-2, -1]
            A = [0 -1; -1 0]
            c = [0,0]
            constr_cones = [(:NonNeg,1:1),(:NonPos,2:2)]
            var_cones = [(:SOC,1:2)]

            m = MathProgBase.ConicModel(solver)
            MathProgBase.loadproblem!(m, c, A, b, constr_cones, var_cones)
            MathProgBase.optimize!(m)
            @test MathProgBase.status(m) == :Infeasible

            if duals
                y = MathProgBase.getdual(m)
                @test y[1] > 0
                @test y[2] < 0
                @test (A'y)[1]/abs((A'y)[2]) ≥ 1-tol
                @test -dot(b,y) > 0
            end
        end

        if duals
            # Problem SOC4
            # Like SOCINT1 but with copies of variables and integrality relaxed
            # Tests out-of-order indices in cones
            @testset "SOC4" begin
                b = [1.0,0.0,0.0]
                A = [1.0 0.0 0.0 0.0 0.0
                     0.0 1.0 0.0 -1.0 0.0
                     0.0 0.0 1.0 0.0 -1.0]
                c = [0.0,-2.0,-1.0,0.0,0.0]
                constr_cones = [(:Zero,[1]),(:Zero,[2,3])]
                var_cones = [(:SOC,[1,4,5]),(:Free,[2,3])]

                m = MathProgBase.ConicModel(solver)
                MathProgBase.loadproblem!(m, c, A, b, constr_cones, var_cones)
                MathProgBase.optimize!(m)
                @test MathProgBase.status(m) == :Optimal

                x = MathProgBase.getsolution(m)
                y = MathProgBase.getdual(m)
                s = MathProgBase.getvardual(m)

                @test isapprox(dot(c,x), -dot(y,b), atol=tol)
                @test x[1]^2 ≥ x[4]^2 + x[5]^2 - tol
                @test s[1]^2 ≥ s[4]^2 + s[5]^2 - tol
                @test isapprox(s[2], 0.0, atol=tol)
                @test isapprox(s[3], 0.0, atol=tol)
                @test isapprox(norm((c+A'y) - s), 0.0, atol=tol)
            end
        end
    end
end

function conicSOCRotatedtest(solver::MathProgBase.AbstractMathProgSolver;duals=false, tol=1e-6)
    @testset "Testing SOCRotated problems through conic interface with $solver" begin
        # Problem SOCRotated1
        # min 0a + 0b - 1x - 1y
        #  st  a            == 1/2
        #  st  b            == 1
        #      2a*b >= y^2+z^2
        @testset "SOCRotated1" begin
            c = [ 0.0, 0.0, -1.0, -1.0]
            A = [ 1.0  0.0   0.0   0.0
                  0.0  1.0   0.0   0.0]
            b = [ 0.5, 1.0]
            m = MathProgBase.ConicModel(solver)
            MathProgBase.loadproblem!(m,c,A,b,[(:Zero,1:2)],[(:SOCRotated,1:4)])
            MathProgBase.optimize!(m)
            @test MathProgBase.status(m) == :Optimal
            @test isapprox(MathProgBase.getobjval(m), -sqrt(2.0), atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[1], 0.5, atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[2], 1.0, atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[3], 1.0/sqrt(2.0), atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[4], 1.0/sqrt(2.0), atol=tol)
            if duals
                d = MathProgBase.getdual(m)
                dualobj = -dot(b,d)
                @test isapprox(dualobj, -sqrt(2.0), atol=tol)
                @test all(d .>= 0)
                vardual = MathProgBase.getvardual(m)
                @test isapprox(norm(vardual - (c+ A'd)), 0.0, atol=tol)
                @test 2*vardual[1]*vardual[2] ≥ vardual[3]^2 + vardual[4]^2 - tol
            end
        end

        # Problem SOCRotated1A - Problem SOCRotated1 with a and b substituted
        # min          -y - z
        #  st [0.5] - [      ] SOCRotated
        #     [1.0] - [      ] SOCRotated
        #     [0.0] - [-y    ] SOCRotated
        #     [0.0] - [    -z] SOCRotated
        @testset "SOCRotated1A" begin
            c = [-1.0,-1.0]
            A = [0.0 0.0; 0.0 0.0; -1.0 0.0; 0.0 -1.0]
            b = [0.5, 1.0, 0.0, 0.0]
            m = MathProgBase.ConicModel(solver)
            MathProgBase.loadproblem!(m,c,A,b,
            [(:SOCRotated,1:4)],
            [(:Free,1:2)])
            MathProgBase.optimize!(m)
            @test MathProgBase.status(m) == :Optimal
            @test isapprox(MathProgBase.getobjval(m), -sqrt(2.0), atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[1], 1.0/sqrt(2.0), atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[2], 1.0/sqrt(2.0), atol=tol)
            if duals
                d = MathProgBase.getdual(m)
                @test 2*d[1]*d[2] ≥ d[3]^2 + d[4]^2 - tol
                dualobj = -dot(b,d)
                @test isapprox(dualobj, -sqrt(2.0), atol=tol)

                vardual = MathProgBase.getvardual(m)
                @test isapprox(norm(vardual - (c+ A'd)), 0.0, atol=tol)
                @test isapprox(norm(vardual), 0.0, atol=tol)
            end
        end

        # Problem SOCRotated2 - Infeasible
        # min 0
        # s.t. z ≥ 2
        #      x ≤ 1
        #      y = 1/2
        #      z^2 ≤ 2x*y
        # in conic form:
        # min 0
        # s.t. -2 + z ∈ R₊
        #      -1 + x ∈ R₋
        #     1/2 - y ∈ {0}
        #       (x,y,z) ∈ SOCRoated
        @testset "SOCRotated2 infeasible" begin
            b = [-2, -1, 1/2]
            A = [0 0 -1; -1 0 0; 0 1 0]
            c = [0,0,0]
            constr_cones = [(:NonNeg,1:1),(:NonPos,2:2),(:Zero,3:3)]
            var_cones = [(:SOCRotated,1:3)]

            m = MathProgBase.ConicModel(solver)
            MathProgBase.loadproblem!(m, c, A, b, constr_cones, var_cones)
            MathProgBase.optimize!(m)
            @test MathProgBase.status(m) == :Infeasible

            if duals
                y = MathProgBase.getdual(m)
                @test y[1] > 0
                @test y[2] < 0
                vardual = MathProgBase.getvardual(m)
                @test isapprox(norm(A'y-vardual), 0.0, atol=tol)
                @test 2*vardual[1]*vardual[2] ≥ vardual[3]^2
                @test -dot(b,y) > 0
            end
        end
    end
end

function conicSOCINTtest(solver::MathProgBase.AbstractMathProgSolver;tol=1e-6)
    @testset "Testing SOCINT problems through conic interface with $solver" begin
        # Problem SOCINT1
        # min 0x - 2y - 1z
        #  st  x            == 1
        #      x >= ||(y,z)||
        #      (y,z) binary
        m = MathProgBase.ConicModel(solver)
        MathProgBase.loadproblem!(m,
        [ 0.0, -2.0, -1.0],
        [ 1.0   0.0   0.0],
        [ 1.0],
        [(:Zero,1)],
        [(:SOC,1:3)])
        MathProgBase.setvartype!(m, [:Cont,:Bin,:Bin])
        MathProgBase.optimize!(m)
        @test MathProgBase.status(m) == :Optimal
        @test isapprox(MathProgBase.getobjval(m), -2, atol=tol)
        @test isapprox(MathProgBase.getsolution(m)[1], 1.0, atol=tol)
        @test isapprox(MathProgBase.getsolution(m)[2], 1.0, atol=tol)
        @test isapprox(MathProgBase.getsolution(m)[3], 0.0, atol=tol)
    end
end

function conicEXPtest(solver::MathProgBase.AbstractMathProgSolver;duals=false, tol=1e-6)
    @testset "Testing EXP problems through conic interface with $solver" begin
        # Problem EXP1 - ExpPrimal
        # min x + y + z
        #  st  y e^(x/y) <= z, y > 0 (i.e (x, y, z) are in the exponential primal cone)
        #      x == 1
        #      y == 2
        @testset "EXP1" begin
            m = MathProgBase.ConicModel(solver)
            c = [1.0, 1.0, 1.0]
            A = [0.0 1.0 0.0;
                 1.0 0.0 0.0]
            b = [2.0, 1.0]
            MathProgBase.loadproblem!(m, c, A, b, [(:Zero,1:2)], [(:ExpPrimal, 1:3)])
            MathProgBase.optimize!(m)
            @test MathProgBase.status(m) == :Optimal
            @test isapprox(MathProgBase.getobjval(m), (2*exp(1/2)+3), atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[1], 1.0, atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[2], 2.0, atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[3], 2*exp(1/2), atol=tol)

            if duals
                d = MathProgBase.getdual(m)
                @test isapprox(-dot(b,d), MathProgBase.getobjval(m), atol=tol)
                u,v,w = c+A'd
                # should belong to the ExpDual cone
                @test u < 0
                @test -u*exp(v/w) ≤ exp(1)*w + tol
                s = MathProgBase.getvardual(m)
                @test isapprox(norm(s - (c+A'd)), 0.0, atol=tol)
            end

            # Permute indices - used to cause Mosek failure
            m = MathProgBase.ConicModel(solver)
            c = [1.0, 1.0, 1.0]
            A = [0.0 1.0 0.0;
                 0.0 0.0 1.0]
            b = [2.0, 1.0]
            MathProgBase.loadproblem!(m, c, A, b, [(:Zero,1:2)], [(:ExpPrimal, [3,2,1])])
            MathProgBase.optimize!(m)
            @test MathProgBase.status(m) == :Optimal
            @test isapprox(MathProgBase.getobjval(m), (2*exp(1/2)+3), atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[3], 1.0, atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[2], 2.0, atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[1], 2*exp(1/2), atol=tol)

            if duals
                d = MathProgBase.getdual(m)
                @test isapprox(-dot(b,d), MathProgBase.getobjval(m), atol=tol)
                w,v,u = c+A'd
                # should belong to the ExpDual cone
                @test u < 0
                @test -u*exp(v/w) ≤ exp(1)*w + tol
                s = MathProgBase.getvardual(m)
                @test isapprox(norm(s - (c+A'd)), 0.0, atol=tol)
            end
        end

        # Problem EXP1A - ExpPrimal
        # Same as previous, except we put :ExpPrimal on constr_cones
        @testset "EXP1A" begin
            m = MathProgBase.ConicModel(solver)
            c = [1.0, 1.0, 1.0]
            A = [0.0 1.0 0.0;
                 1.0 0.0 0.0;
                -1.0 0.0 0.0;
                 0.0 -1.0 0.0;
                 0.0 0.0 -1.0]
            b = [2.0, 1.0, 0.0, 0.0, 0.0]
            MathProgBase.loadproblem!(m, c, A, b,
            [(:Zero,1:2), (:ExpPrimal, 3:5)],
            [(:Free, [2,3,1])])
            MathProgBase.optimize!(m)
            @test MathProgBase.status(m) == :Optimal
            @test isapprox(MathProgBase.getobjval(m), (2*exp(1/2)+3), atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[1], 1.0, atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[2], 2.0, atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[3], 2*exp(1/2), atol=tol)

            if duals
                d = MathProgBase.getdual(m)
                @test isapprox(-dot(b,d), MathProgBase.getobjval(m), atol=tol)
                @test isapprox(norm(c+A'd), 0.0, atol=tol)
                u,v,w = d[3:5]
                # should belong to the ExpDual cone
                @test u < 0
                @test -u*exp(v/w) ≤ exp(1)*w + tol
                s = MathProgBase.getvardual(m)
                @test isapprox(norm(s), 0.0, atol=tol)
            end

            # Permute indices - used to cause Mosek failure
            m = MathProgBase.ConicModel(solver)
            c = [1.0, 1.0, 1.0]
            A = [0.0 1.0 0.0;
                -1.0 0.0 0.0;
                 1.0 0.0 0.0;
                 0.0 -1.0 0.0;
                 0.0 0.0 -1.0]
            b = [2.0, 0.0, 1.0, 0.0, 0.0]
            MathProgBase.loadproblem!(m, c, A, b,
            [(:Zero,[1,3]), (:ExpPrimal, [2,4,5])],
            [(:Free, [3,2,1])])
            MathProgBase.optimize!(m)
            @test MathProgBase.status(m) == :Optimal
            @test isapprox(MathProgBase.getobjval(m), (2*exp(1/2)+3), atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[1], 1.0, atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[2], 2.0, atol=tol)
            @test isapprox(MathProgBase.getsolution(m)[3], 2*exp(1/2), atol=tol)

            if duals
                d = MathProgBase.getdual(m)
                @test isapprox(-dot(b,d), MathProgBase.getobjval(m), atol=tol)
                @test isapprox(norm(c+A'd), 0.0, atol=tol)
                u,v,w = d[2],d[4],d[5]
                # should belong to the ExpDual cone
                @test u < 0
                @test -u*exp(v/w) ≤ exp(1)*w + tol
                s = MathProgBase.getvardual(m)
                @test isapprox(norm(s), 0.0, atol=tol)
            end
        end

        # Problem EXP2
        # A problem where ECOS was failing
        if duals
            @testset "EXP2" begin
                c = [0.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0]
                b = [0.0,1.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0]

                I = [8,11,1,4,9,12,1,4,10,13,3,7,6,7,7,8,11,14,15,9,12,14,10,13,14]
                J = [1,1,2,2,2,2,3,3,3,3,4,4,5,5,6,7,7,7,7,8,8,8,9,9,9]
                V = [-1.0,1.0,-1.0,-1.0,-1.0,1.0,-1.0,1.0,-1.0,1.0,-1.0,-0.5,-1.0,-0.5,1.0,-0.3,-0.3,1.0,-1.0,-0.3,-0.3,1.0,-0.3,-0.3,1.0]
                A = sparse(I, J, V, length(b), length(c))

                cone_con = [(:ExpPrimal,1:3),(:ExpPrimal,4:6),(:Zero,7:7),(:NonNeg,8:10),(:NonNeg,11:13),(:NonNeg,14:14),(:Zero,15:15)]
                cone_var = [(:Free,1:9)]
                m = MathProgBase.ConicModel(solver)
                MathProgBase.loadproblem!(m, c, A, b, cone_con, cone_var)
                MathProgBase.optimize!(m)

                primalobj = MathProgBase.getobjval(m)
                primalsol = MathProgBase.getsolution(m)
                # check primal feasibility
                conval = b-A*primalsol
                @test conval[3] + tol ≥ conval[2]*exp(conval[1]/conval[2])
                @test conval[6] + tol ≥ conval[5]*exp(conval[4]/conval[5])
                @test isapprox(conval[7], 0.0, atol=tol)
                for i in 8:14
                    @test conval[i] ≥ -tol
                end

                @test isapprox(conval[15], 0.0, atol=tol)

                y = MathProgBase.getdual(m)
                dualobj = -dot(b,y)
                @test isapprox(primalobj, dualobj, atol=tol)
                dualconval = c + A'*y
                var = MathProgBase.getvardual(m)
                @test isapprox(dualconval, var, atol=tol)
                for i in 1:9
                    @test isapprox(dualconval[i], 0.0, atol=tol)
                end
                if y[1] == 0.0
                    @test y[2] ≥ -tol
                    @test y[3] ≥ -tol
                else
                    @test y[1] < 0.0
                    @test y[3] > 0.0
                    @test -y[1]*log(-y[1]/y[3]) + y[1] - y[2] ≤ tol
                end
                if y[4] == 0.0
                    @test y[5] ≥ -tol
                    @test y[6] ≥ -tol
                else
                    @test y[4] < 0.0
                    @test y[6] > 0.0
                    @test -y[4]*log(-y[4]/y[6]) + y[4] - y[5] ≤ tol
                end
                for i in 8:14
                    @test y[i] ≥ -tol
                end
            end
        end

        # Problem EXP3
        # Another problem where ECOS was failing
        if duals
            @testset "EXP3" begin
                b = [4.0,0.0,1.0,0.0,5.0]
                c = [-1.0,0.0]

                I = [1,2,4,5]
                J = [1,1,2,2]
                V = [2.0,-1.0,-1.0,1.0]
                A = sparse(I,J,V,5,2)
                cone_var = [(:Free,[1,2])]
                cone_con = [(:NonNeg,[1]),(:ExpPrimal,[2,3,4]),(:NonNeg,[5])]
                m = MathProgBase.ConicModel(solver)
                MathProgBase.loadproblem!(m, c, A, b, cone_con, cone_var)
                MathProgBase.optimize!(m)

                primalobj = MathProgBase.getobjval(m)
                primalsol = MathProgBase.getsolution(m)
                # check primal feasibility
                conval = b-A*primalsol
                @test conval[4] + tol ≥ conval[3]*exp(conval[2]/conval[3])
                @test conval[1] ≥ -tol
                @test conval[5] ≥ -tol

                y = MathProgBase.getdual(m)
                dualobj = -dot(b,y)
                @test isapprox(primalobj, dualobj, atol=tol)
                dualconval = c + A'*y
                var = MathProgBase.getvardual(m)
                @test isapprox(norm(dualconval - var), 0.0, atol=tol)
                for i in 1:2
                    @test isapprox(dualconval[1], 0.0, atol=tol)
                end
                @test y[1] ≥ -tol
                @test -y[2]*log(-y[2]/y[4]) + y[2] - y[3] ≤ tol
                @test y[5] ≥ -tol
            end
        end

        # Problem EXPSOC1
        # A problem where MOSEK was failing
        @testset "EXPSOC1" begin
            m = MathProgBase.ConicModel(solver)
            c = [0.0, 0.0, 1.0, 0.0, 0.0]
            A = [-1.0 0.0 0.0 0.0 0.0; 1.0 3.0 1.0 0.0 0.0; 0.0 -1.0 0.0 0.0 0.0; -1.0 0.0 0.0 0.0 0.0; 2.0 3.0 0.0 0.0 0.0; 0.0 0.0 0.0 -1.0 0.0; 0.0 0.0 0.0 1.0 0.0; -2.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 -1.0; 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 -1.0 0.0; 0.0 1.0 0.0 0.0 1.0]
            b = [0.0, 0.0, -1.0, 2.0, 30.0, 1.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 7.0]
            con_cones = Tuple{Symbol,Array{Int64,1}}[(:NonNeg, [1]), (:Zero, [2]), (:NonNeg, [3]), (:NonNeg, [4]), (:NonNeg, [5]), (:SOC, [6, 7, 8]), (:NonNeg, [9]), (:ExpPrimal, [12, 11, 10]), (:NonNeg, [13])]
            var_cones = Tuple{Symbol,Array{Int64,1}}[(:Free, [1, 2, 3, 4, 5])]
            MathProgBase.loadproblem!(m, c, A, b, con_cones, var_cones)
            MathProgBase.optimize!(m)
            @test MathProgBase.status(m) == :Optimal
            @test isapprox(MathProgBase.getobjval(m), -18.082226, atol=tol)

            if duals
                d = MathProgBase.getdual(m)
                @test isapprox(-dot(b,d), MathProgBase.getobjval(m), atol=tol)
                u,v,w = d[12],d[11],d[10]
                # should belong to the ExpDual cone
                @test u < 0
                @test -u*exp(v/w) ≤ exp(1)*w + tol
            end
        end
    end
end

function conicSDPtest(solver::MathProgBase.AbstractMathProgSolver;duals=true, tol=1e-6)
    s2 = sqrt(2)
    @testset "Testing SDP problems through conic interface with $solver" begin
        # Problem SDP1 - sdo1 from MOSEK docs
        # From Mosek.jl/test/mathprogtestextra.jl, under license:
        #   Copyright (c) 2013 Ulf Worsoe, Mosek ApS
        #   Permission is hereby granted, free of charge, to any person obtaining a copy of this
        #   software and associated documentation files (the "Software"), to deal in the Software
        #   without restriction, including without limitation the rights to use, copy, modify, merge,
        #   publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons
        #   to whom the Software is furnished to do so, subject to the following conditions:
        #   The above copyright notice and this permission notice shall be included in all copies or
        #   substantial portions of the Software.
        #
        #     | 2 1 0 |
        # min | 1 2 1 | . X + x1
        #     | 0 1 2 |
        #
        #
        # s.t. | 1 0 0 |
        #      | 0 1 0 | . X + x1 = 1
        #      | 0 0 1 |
        #
        #      | 1 1 1 |
        #      | 1 1 1 | . X + x2 + x3 = 1/2
        #      | 1 1 1 |
        #
        #      (x1,x2,x3) in C^3_q
        #      X in C_sdp
        #
        @testset "SDP1" begin
            m = MathProgBase.ConicModel(solver)
            #     x1   x2   x3    X11  X21  X31  X22  X32  X33
            c = [ 1.0, 0.0, 0.0,  2.0, s2,  0.0, 2.0, s2,  2.0 ]
            A = [ 1.0  0.0  0.0   1.0  0.0  0.0  1.0  0.0  1.0 ;  # A1
                  0.0  1.0  1.0   1.0  s2   s2   1.0  s2   1.0 ]  # A2
            b = [ 1.0, 0.5 ]

            MathProgBase.loadproblem!(m, c, A, b, [(:Zero,1:2)], [(:SOC,1:3),(:SDP,4:9)] )
            MathProgBase.optimize!(m)
            @test MathProgBase.status(m) == :Optimal
            pobj = MathProgBase.getobjval(m)
            @test isapprox(pobj, 7.05710509e-01, atol=tol)

            xx = MathProgBase.getsolution(m)
            x123 = xx[1:3]
            X = xx[4:9]

            if duals
                y = MathProgBase.getdual(m)
                # Check primal objective
                comp_pobj = dot(X,[2.0,s2,0.0, 2.0,s2, 2.0]) + x123[1]
                # Check dual objective
                comp_dobj = -dot(y,b)
                @test isapprox(comp_pobj, comp_dobj, atol=tol)

                var = c + A' * y
                @test (var[2]^2 + var[3]^2 - var[1]^2) < tol # (var[1],var[2],var[3]) in SOC
                s = MathProgBase.getvardual(m)
                @test isapprox(norm(var - s), 0.0, atol=tol)

                M = [s[4]    s[5]/s2 s[6]/s2
                     s[5]/s2 s[7]    s[8]/s2
                     s[6]/s2 s[8]/s2 s[9]]

                @test eigmin(M) > -tol
            end

            # Reverse row and column indices to check nonconsecutive conic indices
            m = MathProgBase.ConicModel(solver)
            row = collect(2:-1:1)
            col = collect(9:-1:1)
            MathProgBase.loadproblem!(m, c[col], A[row,col], b[row], [(:Zero,row)], [(:SOC,col[1:3]),(:SDP,col[4:9])] )
            MathProgBase.optimize!(m)
            @test MathProgBase.status(m) == :Optimal
            pobj = MathProgBase.getobjval(m)
            @test isapprox(pobj, 7.05710509e-01, atol=tol)

            xx = MathProgBase.getsolution(m)
            x123 = xx[col[1:3]]
            X = xx[col[4:9]]

            if duals
                y = MathProgBase.getdual(m)
                # Check primal objective
                comp_pobj = dot(X,[2.0,s2,0.0, 2.0,s2, 2.0]) + x123[1]
                # Check dual objective
                comp_dobj = -dot(y,b[row])
                @test isapprox(comp_pobj, comp_dobj, atol=tol)

                var = c[col] + A[row,col]' * y
                @test (var[col[2]]^2 + var[col[3]]^2 - var[col[1]]^2) < tol # (var[1],var[2],var[3]) in SOC
                s = MathProgBase.getvardual(m)
                @test isapprox(norm(var - s), 0.0, atol=tol)

                M = [s[col[4]]    s[col[5]]/s2 s[col[6]]/s2
                     s[col[5]]/s2 s[col[7]]    s[col[8]]/s2
                     s[col[6]]/s2 s[col[8]]/s2 s[col[9]]]

                @test eigmin(M) > -tol
            end
        end

        # Problem SDP2
        # Caused getdual to fail on SCS and Mosek
        @testset "SDP2" begin
            m = MathProgBase.ConicModel(solver)

            c = [-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-1.0]
            b = [10.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]
            I = [1,2,8,9,10,11,1,3,8,9,10,11,1,4,8,9,10,11,1,5,8,1,6,8,9,10,11,1,7,8,9,10,11,8,10]
            J = [1,1,1,1,1,1,2,2,2,2,2,2,3,3,3,3,3,3,4,4,4,5,5,5,5,5,5,6,6,6,6,6,6,7,7]
            V = [1.0,1.0,-0.44999999999999996,0.45000000000000007,-0.4500000000000001,0.0,1.0,1.0,-0.7681980515339464,0.31819805153394654,-0.13180194846605373,0.0,1.0,1.0,-0.9000000000000001,0.0,0.0,0.0,1.0,1.0,-0.22500000000000003,1.0,1.0,-0.11250000000000003,0.1125,-0.11249999999999999,0.0,1.0,1.0,0.0,0.0,-0.22500000000000003,0.0,1.0,1.0]
            A = sparse(I, J, V, length(b), length(c))
            cone_con = [(:NonNeg,[1]),(:NonPos,[2,3,4,5,6,7]),(:SDP,8:10),(:Zero,11:11)]
            cone_var = [(:NonNeg,[1,2,3,4,5,6]),(:Free,[7])]

            MathProgBase.loadproblem!(m, c, A, b, cone_con, cone_var)
            MathProgBase.optimize!(m)

            @test MathProgBase.status(m) == :Optimal
            pobj = MathProgBase.getobjval(m)

            x = MathProgBase.getsolution(m)
            @test all(x[1:6] .> -tol)
            con = b - A * x
            @test eigmin([con[8] con[9]/s2 ; con[9]/s2 con[10]]) > -tol
            @test con[1] >= -tol
            @test all(con[2:7] .<= tol)
            @test isapprox(con[11], 0.0, atol=tol)

            if duals
                y = MathProgBase.getdual(m)
                @test eigmin([y[8] y[9]/s2 ; y[9]/s2 y[10]]) > -tol
                @test y[1] >= -tol
                @test all(y[2:7] .<= tol)
                var = c + A' * y
                @test all(var[1:6] .>= -tol)
                @test isapprox(var[7], 0.0, atol=tol)
                dobj = -dot(y,b)
                @test isapprox(pobj, dobj, atol=tol)
                s = MathProgBase.getvardual(m)
                @test isapprox(norm(s - var), 0.0, atol=tol)
            end
        end
    end
end


text:
x = 5
y

ext: jl


text:
module OandaTrade

include("Order.jl")
include("Position.jl")
include("Trade.jl")
include("Account.jl")
include("Config.jl")
include("Instrument.jl")
include("Transaction.jl")
include("Pricing.jl")

end

max_stars_repo_name: CalebDepatie/Julianda

max_forks_repo_name: CalebDepatie/Julianda


text:
mutable struct ErgodicPolicyR2 <: FEBOL.Policy
    execution_horizon::Int      # how long should we execute
    horizon_index::Int          #

    em::ErgodicManagerR2
    tm::TrajectoryManager

    cache::Array{Float64,3}
    visualize::Bool

    solver::Symbol

    xd::VVF
    ud::VVF

    verbose::Bool
end

function ErgodicPolicyR2(f;
                       K::Int = 5,
                       planning_horizon::Int = 50,
                       execution_horizon::Int = 5,
                       visualize::Bool=true,
                       solver::Symbol=:PTO,
                       verbose::Bool=false,
                       R::Matrix{Float64}=0.01 * eye(2)
                      )

    # ergodic manager
    N = planning_horizon
    d = Domain(f.n)
    em = ErgodicManagerR2(d, K)

    # trajectory manager
    ci = ConstantInitializer([0.0,0.0])
    h = 1.
    x0 = [0.5,0.5]
    tm = TrajectoryManager(x0,h,N,ci)
    tm.barrier_cost = 1.0

    tm.R = R

    cache = make_cache(f)

    return ErgodicPolicyR2(execution_horizon,
                         execution_horizon,  # index should be >= exec hrzn
                         em,
                         tm,
                         cache,
                         visualize,
                         solver,
                         [[0.0]],
                         [[0.0]],
                         verbose
                        )
end

function FEBOL.reset!(p::ErgodicPolicyR2)
    p.horizon_index = p.execution_horizon
end

function FEBOL.action(m::SearchDomain, x::Vehicle, o, f::DF, p::ErgodicPolicyR2)

    L = m.length
    if p.horizon_index >= p.execution_horizon

        # generate EID from filter
        p.em.phi = mutual_information(f, p.cache)

        # normalize and decompose
        normalize!(p.em)
        decompose!(p.em)

        # right now, let's just make a new trajectory each time
        p.tm.x0 = [x.x/L, x.y/L]
        if p.solver == :SMC
            p.xd, p.ud = smc_trajectory(p.em, p.tm, umax=5/L)
        else
            tic()
            p.xd, p.ud = pto_trajectory(p.em, p.tm, verbose=p.verbose, max_iters=300)
            t = toq()
            if p.verbose
                println("solve time = ", t)
            end
        end

        p.horizon_index = 0
    end

    p.horizon_index += 1

    a = (p.ud[p.horizon_index][1] * L, p.ud[p.horizon_index][2] * L, 0.0)

    if p.verbose
        println("a = ", a)
    end

    if p.visualize
        xd2 = deepcopy(p.xd)
        for xi in xd2
            xi[1] *= L
            xi[2] *= L
        end
        figure("Simulation")
        plot_trajectory(xd2)

        figure("MI")
        ErgodicControlPlots.cla()
        plot(p.em, p.xd)

        # push this into some sort of thing
    end

    return a
end


text:
@testset "Algebra" begin
    @testset "LorentzVector" begin
        zero = Vec4(0., 0., 0., 0.)
        u1 = Vec4(1., 2., 3., 4.)
        u2 = Vec4(-1., 3., 3., 7.)
        @test -zero == zero
        @test +u1 == u1
        @test +u2 == u2
        @test -(-u1) == u1
        @test -(-u2) == u2
        @test zero + u1 == u1
        @test u1 - u1 == zero
        @test u1 + u1 ≈ 2 * u1
        @test u1 + u1 + u1 ≈ 3 * u1
        @test u1 + u2 ≈ Vec4(0., 5., 6., 11.)
        @test u1 - u2 ≈ Vec4(2., -1., 0., -3.)
        @test 3.5 * u1 ≈ u1 * 3.5
        @test -1.2 * u2 ≈ u2 * (-1.2)
        @test u1 / 2 ≈ Vec4(0.5, 1., 1.5, 2.)
        @test (u1 / 11.) * 11. ≈ u1

        @test zero ⋅ u1 == u1 ⋅ zero == 0.
        @test zero ⋅ u2 == u2 ⋅ zero == 0.
        @test u1 ⋅ u1 ≈ -28.
        @test u2 ⋅ u2 ≈ -66.
        @test u1 ⋅ u2 ≈ -44.

        @test Vec4(1, 0, 0, 1/3) - Vec4(1, 0, 0, 1/5) - Vec4(1, 0, 0, 2/15) ≈ Vec4(-1., 0, 0, 0)
        @test Vec4(1/3, 0, 0, 1) - Vec4(1/5, 0, 0, 1) - Vec4(2/15, 0, 0, 1) ≈ Vec4(0, 0, 0, -1.)
    end;

    @testset "SpatialVector" begin
        zero = Vec3(0., 0., 0.)
        u1 = Vec3(1., 2., 4.)
        u2 = Vec3(-1., 3., 7.)
        @test -zero == zero
        @test +u1 == u1
        @test +u2 == u2
        @test -(-u1) == u1
        @test -(-u2) == u2
        @test zero + u1 == u1
        @test u1 - u1 == zero
        @test u1 + u1 ≈ 2 * u1
        @test u1 + u1 + u1 ≈ 3 * u1
        @test u1 + u2 ≈ Vec3(0., 5., 11.)
        @test u1 - u2 ≈ Vec3(2., -1., -3.)
        @test 3.5 * u1 ≈ u1 * 3.5
        @test -1.2 * u2 ≈ u2 * (-1.2)
        @test u1 / 2 ≈ Vec3(0.5, 1., 2.)
        @test (u1 / 11.) * 11. ≈ u1

        @test zero ⋅ u1 == u1 ⋅ zero == 0.
        @test zero ⋅ u2 == u2 ⋅ zero == 0.
        @test u1 ⋅ u1 ≈ 21.
        @test u2 ⋅ u2 ≈ 59.
        @test u1 ⋅ u2 ≈ 33.

        @test Vec3(0, 1, 0) × Vec3(0, 0, 1) == Vec3(1, 0, 0)
        @test Vec3(0, 0, 1) × Vec3(1, 0, 0) == Vec3(0, 1, 0)
        @test Vec3(1, 0, 0) × Vec3(0, 1, 0) == Vec3(0, 0, 1)
        @test u1 × u1 == zero
        @test u2 × u2 == zero
        @test (u1 × u2) ⋅ u1 == 0
        @test (u1 × u2) ⋅ u2 == 0
        @test u1 × u2 == Vec3(2., -11., 5.)

        @test norm(zero) == 0.
        @test norm(u1) ≈ 4.58257569495584
        @test norm(u2) ≈ 7.681145747868608
        @test norm(Vec3(1, 1, 0)) ≈ √2
        @test norm(Vec3(1, 1, 1)) ≈ √3
        @test norm(Vec3(2, 1, 2)) ≈ 3

        @test normalize(Vec3(1, 1, 0)) ≈ Vec3(1/√2, 1/√2, 0)
        @test normalize(Vec3(1, 1, 1)) ≈ Vec3(1/√3, 1/√3, 1/√3)
        @test normalize(Vec3(2, 2, 1)) ≈ Vec3(2/3, 2/3, 1/3)
    end;
end;


text:
"""
    find_maxlength_subarray(array::Array{Int64,1})

Find the maximum length of a subarray with a given continuous depending sequence of numbers.
For this reason, the minimum and maximum element in every subarray will be tracked and
compared with the incremental distance. So if the difference of the value is equal to the
difference of the increment in the array, the length will be updated.

# Arguments
- `array::Array{Int64,1}`: Unsorted array

# Examples
```julia-repl
julia> import ClassicAlgorithmsCollections
julia> arr = [10, 12, 11, 9, 13, 14, 17, 18, 15]
julia> ClassicAlgorithmsCollections.find_maxlength_subarray(arr)
5
```
"""
function find_maxlength_subarray(array::Array{Int64,1})

    # Initialize result
    max_len = 1
    n = length(array)
    for i in 1:n
        min_array = array[i]
        max_array = array[i]

        for j in (i+1):n
            if min_array > array[j]
                min_array = array[j]
            end
            if max_array < array[j]
                max_array = array[j]
            end
            if ((max_array - min_array) == (j - i))
                if max_len < (max_array - min_array)
                    max_len = max_array - min_array + 1
                end
            end
        end
    end
    return max_len
end

max_stars_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_issues_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_forks_repo_name: Anselmoo/ClassicAlgorithmsCollections


text:
function DocQuerySpace(A::Array{Float64,2},k::Int)
    (U1,S1,V1)=svd(A)
    U=U1[:,1:k]
    S=diagm(S1[1:k])
    Vt=V1'
    V=Vt[1:k,:]
    return S*V,U
end

T=PrepTest()

function SVDModel(A::Array{Float64,2},nq::Int,rank::Int)
    DQ=NewDataMatrix()
    DQ.A=A
    DQ.NumQueries=nq
    (ma,na)=size(DQ.A)
    D=DQ.A[:,DQ.NumQueries+1:na]
    Q=DQ.A[:,1:DQ.NumQueries]
    (md,nd)=size(D)
    (mq,nq)=size(Q)
    Costheta_SVD=zeros(nq,nd)
    #Set a tol level
    tol=linspace(0.01,0.99,20)
    lt=length(tol)
    LSI_TDM_Result=TMResults()
    LSI_TDM_Result=InitResults!(LSI_TDM_Result,nq,lt)
    #Normalize D and Q
    D_norm=Normalize(D)
    Q_norm=Normalize(Q)
    DS,QS=DocQuerySpace(D_norm,rank)
    qindex=zeros(nq)
    for z=1:1:nq
        #println("Z=",z)
        qindex[z]=z
        q=QS'*Q_norm[:,z]
        for i=1:nd
            Costheta_SVD[z:z,i:i]=CosTheta(q,DS[:,i])
        end
        for k=1:1:lt
            (LSI_TDM_Result.Dr[z,k],
             LSI_TDM_Result.Dt[z,k],
             LSI_TDM_Result.Nr[z,k],
             LSI_TDM_Result.Rec[z,k],
             LSI_TDM_Result.Prec[z,k])=
             FindResults(k,z,Costheta_SVD[z,:],tol[k])
        end
    end
    (LSI_TDM_Result.RecFinal,LSI_TDM_Result.PrecFinal)=average_RecPrec(LSI_TDM_Result.Rec,LSI_TDM_Result.Prec,qindex)
    return LSI_TDM_Result.RecFinal,LSI_TDM_Result.PrecFinal
    #plotNew_RecPrec(LSI_TDM_Result.RecFinal,LSI_TDM_Result.PrecFinal,"LSI")
end

function SVDModel(Q_C::Corpus,D_C::Corpus,rank::Int)
    #Q_c and D_c are the query and Document corpuses.
    PreProcess!(Q_C)
    PreProcess!(D_C)
    #Compulsory update the lexicon, to identify the keywords
    update_lexicon!(Q_C)
    update_lexicon!(D_C)
    #Get the lexicon of the documents.
    DocLex=lexicon(D_C)
    #Cant use the above as query vectors since, the queries must also be wrt
    #Documents lexicon, and hence the dimensions will match
    D=full(tdm(D_C))
    (md,nd)=size(D)
    nq=length(Q_C.documents)
    Q=zeros(md,nq)
    for j=1:nq
        Q[:,j]=dtv(Q_C.documents[j],DocLex)'
    end
    QD=[Q D]
    r=SVDModel(QD,nq,rank)
    return r
end

function SVDModel(QueryNum::Int,A::Array{Float64,2},NumQueries::Int,rank::Int)
    DQ=NewDataMatrix()
    DQ.A=A
    DQ.NumQueries=NumQueries
    (ma,na)=size(DQ.A)
    D=DQ.A[:,DQ.NumQueries+1:1063]
    Q=DQ.A[:,QueryNum]
    Qt=Q'
    (nq,mq)=size(Qt)
    (md,nd)=size(D)
    Costheta_SVD=zeros(nq,nd)
    #Set a tol level
    tol=linspace(0.01,0.99,20)
    lt=length(tol)
    LSI_1Query_Result=TMResults()
    LSI_1Query_Result=InitResults!(LSI_1Query_Result,nq,lt)
    D_norm=Normalize(D)
    Q_norm=Normalize(Qt')
    DS,QS=DocQuerySpace(D_norm,rank)
    qindex=zeros(nq)
    for z=1:1:nq
        qindex[z]=z
        q=QS'*Q_norm[:,z]
        for i=1:nd
            Costheta_SVD[z:z,i:i]=q'*DS[:,i]/(norm(q,2)*norm(DS[:,i],2))
        end
        for k=1:1:lt
            (LSI_1Query_Result.Dr[z,k],
             LSI_1Query_Result.Dt[z,k],
             LSI_1Query_Result.Nr[z,k],
             LSI_1Query_Result.Rec[z,k],
             LSI_1Query_Result.Prec[z,k])=
            FindResults(k,z,Costheta_SVD[z,:],tol[k])
        end
    end
    (LSI_1Query_Result.RecFinal,LSI_1Query_Result.PrecFinal)=average_RecPrec(LSI_1Query_Result.Rec,LSI_1Query_Result.Prec,qindex)
    return LSI_1Query_Result
    #plotNew_RecPrec(vec(LSI_1Query_Result.Rec),vec(LSI_1Query_Result.Prec),"LSI")
end

function TrySVD(A::Array{Float64,2},k::Int)
    (U,S,V)=svd(A)
    A_k::Array{Float64,2}
    V_t=diagm(S)*V'
    A_k=U[:,1:k]*V_t[1:k,:]
    #Check the Frobenius norm Between A and A_k
    #FN=vec_norm(A-A_k)
    FN=norm(A-A_k)
    return FN
end


text:
@testset "Test Diethelm PECE algorithms" begin
    fun(x, y) = 1-y
    prob = SingleTermFODEProblem(fun, 1.8, [0, 0], (0, 5))
    sol = solve(prob, 0.5, PECE())

    @test isapprox(sol.u, [0.16153482345602124
    0.5289988135096678
    0.848191872231344
    1.0899386876765813
    1.2276684080819034
    1.2684889407947215
    1.2385617809129195
    1.1700696107638846
    1.0918355406175286
    1.0242165055531318
    0.977769035221195]; atol=1e-3)
end

@testset "Test Matrix discrete method" begin
    fun(x, y) = 1-y
    u0 = [0, 0]
    tspan=collect(0.01:0.01:5)
    target = []
    for i in 0.01:0.01:5
        push!(target, i^0.5*mittleff(0.5, 1.5,-i^0.5))
    end

    singleprob = MultiTermsFODEProblem([1, 1], [0.5, 0], 1, u0, (0, 5))
    result = solve(singleprob, 0.01, FODEMatrixDiscrete())

    @test isapprox(result.u, target; atol=1)

    ########### Yet another test case ##########
    yafun(x, y) = 1-y

    #Analytical solution
    yatarget = []

    #MittagLeffler.jl doesn't support array argument
    for i in 0.01:0.01:20
        push!(yatarget, i^1.8*mittleff(1.8,2.8,-i^1.8))
    end

    highsingleprob = MultiTermsFODEProblem([1, 1], [1.8, 0], 1, u0, (0, 20))
    yaresult = solve(highsingleprob, 0.01, FODEMatrixDiscrete())

    @test isapprox(yaresult.u, yatarget; atol=1)

end

@testset "Test Closed Form method" begin
    rightfun(x)=sin(x^2)
    prob = MultiTermsFODEProblem([1 8 26 73 90], [3.5 3.1 2.3 1.2 0.5], rightfun, [30 90 3], [1 0.3 0], [0, 0], (0, 5))
    sol = solve(prob, 0.5, ClosedForm())

    @test isapprox(sol.u, [ 0.0
    0.0854971609134175
    0.38289836071474537
    0.6454389110668579
    0.2783045139283708
    0.043460434414588106
    0.12348574062908478
    0.04623359232392658
   -0.07835912500702993
    0.26027660764161137
    0.2855329523586738]; atol=1e-3)
end

@testset "Test ClosedFormHankelM method" begin
    u0 = [0, 0, 0, 0]
    rightfun(x)=sin(x^2)
    prob = MultiTermsFODEProblem([1 8 26 73 90], [3.5 3.1 2.3 1.2 0.5], rightfun, [30 90], [1 0.3], u0, (0, 1))
    result = solve(prob, 0.5, ClosedFormHankelM())

    @test isapprox(result.u, [0.0; 0.08402140107687359; 0.3754974742112727]; atol=1e-3)
end

@testset "Test GL method for FODESystem" begin
    h=0.5; tf=1
    alpha = [0.99, 0.99, 0.99]
    x0 = [1, 0, 1]
    function testf!(du, u, p, t)
        a, b, c = 10, 28, 8/3
        du[1] = a*(u[2]-u[1])
        du[2] = u[1]*(b-u[3])-u[2]
        du[3] = u[1]*u[2]-c*u[3]
    end
    prob = FODESystem(testf!, alpha, x0, (0, tf))
    result = solve(prob, h, GL())
    @test isapprox(result, [1.0 0.0 1.0
    -4.04478 13.5939 -0.352607
    84.8074 -51.1251 -27.5541]; atol=1e-4)
end

@testset "Test GL method" begin
    fun(x, y) = 1-y
    prob = SingleTermFODEProblem(fun, 0.5, 0, (0, 1))
    sol = solve(prob, 0.1, GL())

    @test isapprox(sol.u, [0.0
    0.31622776601683794
    0.3743416490252569
    0.42454983788325495
    0.4608051796660425
    0.4897286932245971
    0.5136308879844076
    0.5339402943673064
    0.5515444532932976
    0.5670394859746068
    0.5808448788127054]; atol=1e-3)
end

@testset "Test Nonlinear method" begin
    function chua!(du, x, p, t)
        a = 10.725; b = 10.593
        c = 0.268
        m0 = -1.1726
        m1 = -0.7872
        du[1] = a*(x[2]-x[1]-(m1*x[1]+0.5*(m0-m1)*(abs(x[1]+1)-abs(x[1]-1))))
        du[2] = x[1]-x[2]+x[3]
        du[3] = -b*x[2]-c*x[3]
    end
    α = [0.93, 0.99, 0.92];
    x0 = [0.2; -0.1; 0.1];
    h = 0.1; tspan = (0, 0.5);
    prob = FODESystem(chua!, α, x0, tspan)
    result = solve(prob, h, NonLinearAlg())

    @test isapprox(result, [ 0.2 -0.1 0.1
    0.11749    -0.0590683  0.224134
    0.074388   -0.018475   0.282208
    0.0733938   0.0192931  0.286636
    0.117483    0.0534393  0.246248
    0.210073    0.0844175  0.168693]; atol=1e-3)
end

@testset "Test Product Integral Explicit method" begin
    fun(t, y)=1-y
    prob = SingleTermFODEProblem(fun, 0.5, 0, (0, 5))
    sol=solve(prob, 0.5, PIEX())

    @test isapprox(sol.u, [0.0
    0.5840920370824765
    0.6497603208539517
    0.6965874966593231
    0.7309421360673156
    0.7569810570702951
    0.7773192302138408
    0.7936279277715693
    0.8070043316131614
    0.8181895332990293
    0.8276979913681598]; atol=1e-4)
end
#=
@testset "Test Product Integral Implicit method" begin
    fun(t, y)=1-y
    prob = SingleTermFODEProblem(fun, 0.5, 0, 1)
    sol=solve(prob, 0.1, PIIM())

    @test isapprox(sol, [ 0.0
    0.504626504404032
    0.543454077521492
    0.5760954086164125
    0.6028545521580736
    0.6251322531561498
    0.6440148324633577
    0.6602762702077741
    0.6744705915156369
    0.6870025519994545]; atol=1e-4)
end

@testset "Test Product Integral Implicit Trapezoidal method" begin
    fun(t, y)=1-y
    prob = SingleTermFODEProblem(fun, 0.5, 0, 1)
    sol=solve(prob, 0.1, PITrap())

    @test isapprox(sol, [0.0
    0.504626504404032
    0.5185925289496215
    0.5467126570149609
    0.571261200170596
    0.5923262209457972
    0.6105232343055992
    0.6264159719331727
    0.6404461835166287
    0.652951997319657]; atol=1e-4)
end
=#
@testset "Test ChebSpectral method" begin
    n = 10
    rightfun(x) = gamma(7+9/17)/gamma(7+9/17-0.5)*(1+x)^(6+9/17-0.5)
    prob = SingleTermFODEProblem(rightfun, 0.5, 0, 1)
    sol = solve(prob, n, ChebSpectral())

    @test isapprox(sol.u, [-2.103832955612273e-8
    7.611080159690798e-5
    0.01082543428432764
    0.2878301921177669
    2.8447345612767694
   14.117972832745936
   40.9991958543131
   75.6356722636105
   92.37379693143299]; atol=1e-4)
end

@testset "Test Product Integration explicit method for multi-terms FODE" begin
    rightfunex(x, y) = 172/125*cos(4/5*x)
    PIEXMultiprob = MultiTermsFODEProblem([1, 1/16, 4/5, 3/2, 1/25, 6/5], [3, 2.5, 2, 1, 0.5, 0], rightfunex, [0, 0, 0, 0, 0, 0], (0, 10))
    PIEXMultisol = solve(PIEXMultiprob, 0.5, PIEX())

    @test isapprox(PIEXMultisol.u, [  0.0
    0.028666666666666663
    0.20802154481466284
    0.5777512805039694
    1.032617182281499
    1.3542099549416406
    1.2992204949109856
    0.7162691469691946
   -0.35553236137048394
   -1.6531692565209024
   -2.748115186408473
   -3.188228316718636
   -2.6777963530465527
   -1.226608005191384
    0.7994418493409547
    2.76345190855897
    3.9650094859689915
    3.9040343901844707
    2.5017274305718367
    0.18249065149394283
   -2.23584183347657]; atol=1e-4)
end

@testset "Test Product integration with predictor-corrector method for multi-terms FODE" begin
    T = 10; h = 0.5
    rightfun(x, y) = 172/125*cos(4/5*x)
    prob = MultiTermsFODEProblem([1, 1/16, 4/5, 3/2, 1/25, 6/5], [3, 2.5, 2, 1, 0.5, 0], rightfun, [0, 0, 0, 0, 0, 0], (0, T))

    sol = solve(prob, h, PIPECE())

    @test isapprox(sol.u, [  0.0
    0.01942631486406702
    0.13703330688409168
    0.3665463925810521
    0.6301167829803188
    0.7968019190472511
    0.7407968242020087
    0.3985067557399453
   -0.19598223658593597
   -0.9100047185814064
   -1.5470928325016606
   -1.9063815429019095
   -1.8472382231556892
   -1.3378611193669085
   -0.47174020197161115
    0.5537908253525861
    1.49062009988213
    2.108384867938473
    2.2577394376960034
    1.9096795116394878
    1.160023434854465]; atol=1e-4)
end

@testset "Test implicit Product integration rectangular type method for multi-terms FODE" begin
    T = 10; h = 0.5
    rightfun(x, y) = 172/125*cos(4/5*x)
    prob = MultiTermsFODEProblem([1, 1/16, 4/5, 3/2, 1/25, 6/5], [3, 2.5, 2, 1, 0.5, 0], rightfun, [0, 0, 0, 0, 0, 0], (0, T))

    sol = solve(prob, h, PIIMRect())

    @test isapprox(sol.u, [0.0
    0.01586240520261297
    0.1119084041841651
    0.2954711964704977
    0.5008017075335748
    0.6241431840961683
    0.5697362571789858
    0.2914052658364976
   -0.1836104680656457
   -0.7558230882351726
   -1.278848877489753
   -1.6008246449445735
   -1.6091002088212767
   -1.2648474636427782
   -0.6175431717278957
    0.2046270023806023
    1.0255234402793145
    1.664252313922362
    1.9791244106539438
    1.9009297334987802
    1.447593455242784]; atol=1e-4)
end

@testset "Test implicit Product integration trapezoidal type method for multi-terms FODE" begin
    T = 10; h = 0.5
    rightfun(x, y) = 172/125*cos(4/5*x)
    prob = MultiTermsFODEProblem([1, 1/16, 4/5, 3/2, 1/25, 6/5], [3, 2.5, 2, 1, 0.5, 0], rightfun, [0, 0, 0, 0, 0, 0], (0, T))

    sol = solve(prob, h, PIIMTrap())

    @test isapprox(sol.u, [0.0
    0.02157284502548659
    0.153474720938717
    0.41530591952310447
    0.723333042644268
    0.9282329990906286
    0.879936169829285
    0.4964927100233818
   -0.1888267323884115
   -1.0211143586249927
   -1.7629487467873093
   -2.167707731932248
   -2.062971791641991
   -1.415305444047155
   -0.3530008235484339
    0.86303702756878
    1.9127384461430172
    2.5119503312180154
    2.498678841978473
    1.882454666603612
    0.8406427306579508]; atol=1e-4)
end

@testset "Test product integration PECE method for FODESystem" begin
    function test(du, u, p, t)
        du[1] = 1-4*u[1]+u[1]^2*u[2]
        du[2] = 3*u[1]-u[1]^2*u[2]
    end

    alpha = [0.8, 0.7]
    tspan=(0, 0.1); h=0.01
    y0=[1.2; 2.8]
    param=[1 3]
    prob = FODESystem(test, alpha, y0, tspan)
    (t, y) = solve(prob, h, PECE())

    @test isapprox(y, [1.2  1.2061   1.21042  1.21421  1.21767  1.22089  1.22392  1.22678  1.2295   1.2321   1.23459
    2.8  2.78118  2.76953  2.7596   2.75064  2.74235  2.73455  2.72715  2.72008  2.71329  2.70673]; atol=1e-3)
end

@testset "Test explicit product integration method for FODESystem" begin
    function test(du, u, p, t)
        du[1] = 1-4*u[1]+u[1]^2*u[2]
        du[2] = 3*u[1]-u[1]^2*u[2]
    end

    alpha = [0.8, 0.7]
    tspan = (0, 0.1); h=0.01
    y0=[1.2; 2.8]
    param=[1 3]
    prob = FODESystem(test, alpha, y0, tspan)
    (t, y) = solve(prob, h, PIEX())

    @test isapprox(y, [ 1.2  1.20626  1.21061  1.21444  1.21793  1.22118  1.22424  1.22713  1.22989  1.23252  1.23504
    2.8  2.78107  2.76943  2.75949  2.75052  2.74221  2.7344   2.72697  2.71988  2.71305  2.70647]; atol=1e-4)
end

@testset "Test FLMMNewtonGregory method" begin
    a=1; mu=4
    function Brusselator(du, u, p, t)
        du[1] = a-(mu+1)*u[1]+u[1]^2*u[2]
        du[2] = mu*u[1]-u[1]^2*u[2]
    end
    alpha=[0.8; 0.8]
    y0=[0.2; 0.03]
    h=0.1; tspan=(0, 0.5)
    prob = FODESystem(Brusselator, alpha, y0, tspan)
    (t, y) = solve(prob, h, FLMMNewtonGregory())

    @test isapprox(y, [ 0.2   0.200531  0.201161  0.201809  0.202453  0.203089
    0.03  0.165554  0.265678  0.355635  0.439544  0.519211]; atol=1e-4)
end

@testset "Test FLMMBDF" begin
    a=1; mu=4
    function Brusselator(du, u, p, t)
        du[1] = a-(mu+1)*u[1]+u[1]^2*u[2]
        du[2] = mu*u[1]-u[1]^2*u[2]
    end
    alpha=[0.8; 0.8]
    y0=[0.2; 0.03]
    h=0.1; tspan=(0, 0.5)
    prob = FODESystem(Brusselator, alpha, y0, tspan)
    (t, y) = solve(prob, h, FLMMBDF())

    @test isapprox(y, [ 0.2   0.200531  0.201161  0.201809  0.202453  0.203089
    0.03  0.165554  0.265678  0.355635  0.439544  0.519211]; atol=1e-4)
end

@testset "Test FLMMTrap" begin
    a=1; mu=4
    function Brusselator(du, u, p, t)
        du[1] = a-(mu+1)*u[1]+u[1]^2*u[2]
        du[2] = mu*u[1]-u[1]^2*u[2]
    end
    alpha=[0.8; 0.8]
    y0=[0.2; 0.03]
    h=0.1; tspan=(0, 0.5)
    prob = FODESystem(Brusselator, alpha, y0, tspan)
    (t, y) = solve(prob, h, FLMMTrap())

    @test isapprox(y, [0.2   0.200531  0.201161  0.201808  0.202452  0.203088
    0.03  0.165554  0.265678  0.355635  0.439545  0.519211]; atol=1e-4)
end


text:
using JuMP
using Combinatorics

# Inspired from Joey Huchette's test in ConvexHull.jl
function crosspolytopetest(lib::PolyhedraLibrary, n)
    m = Model()

    @variable(m, x[1:n])

    for k in 0:n
        for S in combinations(1:n, k)
            Sᶜ = setdiff(1:n, S)
            @constraint(m, sum(x[i] for i in S) - sum(x[i] for i in Sᶜ) ≤ 1)
        end
    end

    poly = polyhedron(m, lib)

    V = zeros(Int, 2n, n)
    for i in 1:n
        V[2i-1, i] = 1
        V[2i, i] = -1
    end
    generator_fulltest(poly, V)
end


text:
const SCHEMA_DIR = joinpath(ROOT, "schema")
const VERSIONS = Dict([
    "2.0" => "v2",
    "3.0" => "v3",
])

mutable struct OpenAPI
    version::String
    info::Dict{String, Any}
    paths::Dict{String, Any}
    optional_fields::Dict{String, Any}

    function OpenAPI(version::String, info::Dict{String, Any}; optional_fields::Dict{String, Any}=Dict{String, Any}())
        this = new()
        this.version = version
        this.info= info
        this.optional_fields = optional_fields
        this.paths = Dict{String, Any}()
        return this
    end
end
struct InvalidSwaggerSpecificationException <:Exception
    issue::JSONSchema.SingleIssue
end

function Base.showerror(io::IO, e::InvalidSwaggerSpecificationException)
    issue = e.issue
    println(io, """\n
    Validation failed:
    path:         $(isempty(issue.path) ? "top-level" : issue.path)
    instance:     $(issue.x)
    schema key:   $(issue.reason)
    schema value: $(issue.val)""")
end

function validate_spec(spec::Dict{String, Any})

    @assert (haskey(spec, "swagger") || haskey(spec, "openapi")) "Field 'swagger' (v2) or 'openapi' (v3) is missing"

    version_name = haskey(spec, "swagger") ? "swagger" : "openapi"
    version = spec[version_name]
    if version == "3.0"
        spec[version_name] = "3.0.0"
    elseif startswith(version, "3.0")
        version = "3.0"
    end
    @assert haskey(VERSIONS, version) "Version $(version) is not supported"

    openapi_schema = JSONSchema.Schema(JSON.parsefile(joinpath(SCHEMA_DIR , "$(VERSIONS[version]).json")))

    issue = JSONSchema.validate(spec, openapi_schema)
    if issue !== nothing
        throw(InvalidSwaggerSpecificationException(issue))
    end
end

function validate_spec(file::String)
    if endswith(file, ".json")
        return validate_spec(JSON.parsefile(file))
    end
    if endswith(file, ".yml")
        return validate_spec(YAML.load_file(file, dicttype=Dict{String, Any}))
    end
end


text:
# mainly I need the language to perform
# 1. experiment analysis
#   1.1 hypothesis tests in frequentist eyes
#         t.test, chisq.test, poisson.test, prop.test, bino.test
#   1.2 in bayesians' eyes naive simulation of gamma, beta distribution
#        for possion and binomial proportion tset
# 2. modelling and interpretation
#        glm, random forest, xgboost, mcmc, neural net
# 3. nice to have
#        call fortran when needed
#        matmul when having spare time for fun
#        GPU directly or indirectly
#        comfy IDE

using Statistics
using HypothesisTests           # t.test ordinary or welch

using Distributions             # Gamma, Beta
using StatsPlots

#1. experiment

a = rand(Beta(400, 600), 10000)
b = rand(Beta(450, 550), 10000)

var(a)
var(b)
EqualVarianceTTest(a, b)

violin([a b])
histogram([a b])
density([a b], label=["A" "B"])
density([a b], layout=(2,1))

marginalhist(a, b)
plot(a, b)
scatter(a, b)

#2. modelling
~

max_stars_repo_name: scienceofdeduction/julia

max_issues_repo_name: scienceofdeduction/julia

max_forks_repo_name: scienceofdeduction/julia


text:
test_gage_min_distance(
    Array[[1., 0., 0.], [0., 0., 1.]],
    Array[[0., 1., 0.], [0., 0., 1.]],
    0.,
    "two array of arrays"
)


text:
import KernelAbstractions: CPU

using ..Mesh.Grids
using ..Mesh.Elements: interpolationmatrix
using ..MPIStateArrays
using ..DGMethods: SpaceDiscretization
using ..TicToc

"""
    writevtk(prefix, Q::MPIStateArray, dg::SpaceDiscretization [, fieldnames];
             number_sample_points = 0)

Write a vtk file for all the fields in the state array `Q` using geometry and
connectivity information from `dg.grid`. The filename will start with `prefix`
which may also contain a directory path. The names used for each of the fields
in the vtk file can be specified through the collection of strings `fieldnames`;
if not specified the fields names will be `"Q1"` through `"Qk"` where `k` is the
number of states in `Q`, i.e., `k = size(Q,2)`.

If `number_sample_points > 0` then the fields are sampled on an equally spaced,
tensor-product grid of points with 'number_sample_points' in each direction and
the output VTK element type is set to by a VTK lagrange type.

When `number_sample_points == 0` the raw nodal values are saved, and linear VTK
elements are used connecting the degree of freedom boxes.
"""
function writevtk(
    prefix,
    Q::MPIStateArray,
    dg::SpaceDiscretization,
    fieldnames = nothing;
    number_sample_points = 0,
)
    vgeo = dg.grid.vgeo
    h_vgeo = array_device(vgeo) isa CPU ? vgeo : Array(vgeo)
    h_Q = array_device(Q) isa CPU ? Q.data : Array(Q)
    writevtk_helper(
        prefix,
        h_vgeo,
        h_Q,
        dg.grid,
        fieldnames;
        number_sample_points = number_sample_points,
    )
    return nothing
end

"""
    writevtk(prefix, Q::MPIStateArray, dg::SpaceDiscretization, fieldnames,
             state_auxiliary::MPIStateArray, auxfieldnames;
             number_sample_points = 0)

Write a vtk file for all the fields in the state array `Q` and auxiliary state
`state_auxiliary` using geometry and connectivity information from `dg.grid`. The
filename will start with `prefix` which may also contain a directory path. The
names used for each of the fields in the vtk file can be specified through the
collection of strings `fieldnames` and `auxfieldnames`.

If `fieldnames === nothing` then the fields names will be `"Q1"` through `"Qk"`
where `k` is the number of states in `Q`, i.e., `k = size(Q,2)`.

If `auxfieldnames === nothing` then the fields names will be `"aux1"` through
`"auxk"` where `k` is the number of states in `state_auxiliary`, i.e., `k =
size(state_auxiliary,2)`.

If `number_sample_points > 0` then the fields are sampled on an equally spaced,
tensor-product grid of points with 'number_sample_points' in each direction and
the output VTK element type is set to by a VTK lagrange type.

When `number_sample_points == 0` the raw nodal values are saved, and linear VTK
elements are used connecting the degree of freedom boxes.
"""
function writevtk(
    prefix,
    Q::MPIStateArray,
    dg::SpaceDiscretization,
    fieldnames,
    state_auxiliary,
    auxfieldnames;
    number_sample_points = 0,
)
    vgeo = dg.grid.vgeo
    device = array_device(Q)
    (h_vgeo, h_Q, h_aux) =
        device isa CPU ? (vgeo, Q.data, state_auxiliary.data) :
        (Array(vgeo), Array(Q), Array(state_auxiliary))
    writevtk_helper(
        prefix,
        h_vgeo,
        h_Q,
        dg.grid,
        fieldnames,
        h_aux,
        auxfieldnames;
        number_sample_points = number_sample_points,
    )
    return nothing
end

reshaped_view(fields, ind, Np_N1, Nq, nelem) =
    ntuple(length(fields)) do i
        fld = reshape(fields[i], Nq..., nelem)
        fld2 = @view fld[ind..., :]
        reshape(fld2, (Np_N1, nelem))
    end

"""
    writevtk_helper(prefix, vgeo::Array, Q::Array, grid, fieldnames)

Internal helper function for `writevtk`
"""
function writevtk_helper(
    prefix,
    vgeo::Array,
    Q::Array,
    grid,
    fieldnames,
    state_auxiliary = nothing,
    auxfieldnames = nothing;
    number_sample_points,
)
    @assert number_sample_points >= 0

    dim = dimensionality(grid)
    N = polynomialorders(grid)
    Nq = N .+ 1

    nelem = size(Q)[end]

    X = grid.x_vtk[1:dim]
    fields = ntuple(j -> (@view Q[:, j, :]), size(Q, 2))

    auxfields =
        isnothing(state_auxiliary) ? () :
        (
            auxfields = ntuple(
                j -> (@view state_auxiliary[:, j, :]),
                size(state_auxiliary, 2),
            )
        )

    # If any dimension are N = 0 we mirror these out to the boundaries for viz
    # purposed
    if any(N .== 0)
        Nq_N1 = max.(Nq, 2)
        Np_N1 = prod(Nq_N1)
        ind = ntuple(i -> N[i] == 0 ? [1, 1] : Colon(), dim)
        fields = reshaped_view(fields, ind, Np_N1, Nq, nelem)
        auxfields = reshaped_view(auxfields, ind, Np_N1, Nq, nelem)
        Nq = Nq_N1
    end

    # Interpolate to an equally spaced grid if necessary
    if number_sample_points > 0
        FT = eltype(Q)
        # If any dimension are N = 0 we manual set (-1, 1) grids
        ξ = ntuple(
            i -> N[i] == 0 ? FT.([-1, 1]) : referencepoints(grid)[i],
            dim,
        )
        ξdst = range(FT(-1); length = number_sample_points, stop = 1)
        I1d = ntuple(i -> interpolationmatrix(ξ[dim - i + 1], ξdst), dim)
        I = kron(I1d...)
        fields = ntuple(i -> I * fields[i], length(fields))
        auxfields = ntuple(i -> I * auxfields[i], length(auxfields))
        X = ntuple(i -> I * X[i], length(X))
        Nq = ntuple(j -> number_sample_points, dim)
    end

    X = ntuple(i -> reshape(X[i], Nq..., nelem), length(X))

    function get_fields(x, fieldnames, name)
        x = ntuple(i -> reshape(x[i], Nq..., nelem), length(x))
        if fieldnames === nothing
            return ntuple(i -> ("$name$i", x[i]), length(x))
        else
            return ntuple(i -> (fieldnames[i], x[i]), length(x))
        end
    end

    fields = get_fields(fields, fieldnames, "Q")
    auxfields = get_fields(auxfields, auxfieldnames, "aux")

    fields = (fields..., auxfields...)
    if number_sample_points > 0
        return writemesh_highorder(
            prefix,
            X...;
            fields = fields,
            realelems = grid.topology.realelems,
        )
    else
        return writemesh_raw(
            prefix,
            X...;
            fields = fields,
            realelems = grid.topology.realelems,
        )
    end
end


text:
# TableRow

"""
    TableRow{M, N}(catcols, contcols, categorydict) <: Block

`Block` for table rows with M categorical and N continuous columns. `data`
is valid if it satisfies the `AbstractRow` interface in Tables.jl, values
present in indices for categorical and continuous columns are consistent,
and `data` is indexable by the elements of `catcols` and `contcols`.
"""
struct TableRow{M,N,T} <: Block
    catcols::NTuple{M}
    contcols::NTuple{N}
    categorydict::T
end

function TableRow(catcols, contcols, categorydict)
    TableRow{length(catcols),length(contcols)}(catcols, contcols, categorydict)
end

function checkblock(block::TableRow, x)
    columns = Tables.columnnames(x)
    (
        all(col -> col ∈ columns, (block.catcols..., block.contcols...)) &&
        all(
            col ->
                haskey(block.categorydict, col) &&
                    (ismissing(x[col]) || x[col] ∈ block.categorydict[col]),
            block.catcols,
        ) &&
        all(col -> ismissing(x[col]) || x[col] isa Number, block.contcols)
    )
end

function mockblock(block::TableRow)
    cols = (block.catcols..., block.contcols...)
    vals = map(cols) do col
        col in block.catcols ? rand(block.categorydict[col]) : rand()
    end
    return NamedTuple(zip(cols, vals))
end

"""
    setup(TableRow, data[; catcols, contcols])

Create a `TableRow` block from data container `data::TableDataset`. If the
categorical and continuous columns are not specified manually, try to
guess them from the dataset's column types.
"""
function setup(::Type{TableRow}, data; catcols = nothing, contcols = nothing)
    catcols_, contcols_ = getcoltypes(data)
    catcols = isnothing(catcols) ? catcols_ : catcols
    contcols = isnothing(contcols) ? contcols_ : contcols

    return TableRow(
        catcols,
        contcols,
        gettransformdict(data, DataAugmentation.Categorify, catcols),
    )
end

function Base.show(io::IO, block::TableRow)
    print(io, ShowCase(block, (:catcols, :contcols), show_params = false, new_lines = true))
end

# ## Interpretation

function showblock!(io, ::ShowText, block::TableRow, obs)
    rowdata = vcat(
        [obs[col] for col in block.catcols],
        [obs[col] for col in block.contcols],
    )
    rownames = [block.catcols..., block.contcols...]
    tabledata = hcat(rownames, rowdata)
    PrettyTables.pretty_table(
        io, tabledata;
        alignment=[:r, :l],
        highlighters=PrettyTables.Highlighter((obs, i, j) -> (j == 2), bold=true),
        noheader=true, tf=PrettyTables.tf_borderless,)
end


text:
struct BasicBlockSimulateState
    trace::Symbol
    score::Symbol
    stmts::Vector{Expr}
end

function process!(ir::BasicBlockIR, state::BasicBlockSimulateState, node::JuliaNode)
    trace = state.trace
    (typ, trace_field) = get_value_info(node)
    push!(state.stmts, quote
        $trace.$trace_field = $(expr_read_from_trace(node, trace))
    end)
end

function process!(ir::BasicBlockIR, state::BasicBlockSimulateState,
                  node::Union{ArgsChangeNode,AddrChangeNode})
    trace = state.trace
    (typ, trace_field) = get_value_info(node)
    push!(state.stmts, quote
        $trace.$trace_field = nothing
    end)
end

function process!(ir::BasicBlockIR, state::BasicBlockSimulateState, node::AddrDistNode)
    trace, score = state.trace, state.score
    addr = node.address
    dist = QuoteNode(node.dist)
    args = get_args(trace, node)
    push!(state.stmts, quote
        $trace.$addr = random($dist, $(args...))
        $score += logpdf($dist, $trace.$addr, $(args...))
        $trace.$is_empty_field = false
    end)
    if has_output(node)
        (_, trace_field) = get_value_info(node)
        push!(state.stmts, quote
            $trace.$trace_field = $trace.$addr
        end)
    end
end

function process!(ir::BasicBlockIR, state::BasicBlockSimulateState, node::AddrGeneratorNode)
    trace, score = state.trace, state.score
    addr = node.address
    gen = QuoteNode(node.gen)
    args = get_args(trace, node)
    call_record = gensym("call_record")
    push!(state.stmts, quote
        $trace.$addr = simulate($gen, $(Expr(:tuple, args...)))
        $call_record = get_call_record($trace.$addr)
        $score += $call_record.score
        $trace.$is_empty_field = $trace.$is_empty_field && !has_choices($trace.$addr)
    end)
    if has_output(node)
        (_, trace_field) = get_value_info(node)
        push!(state.stmts, quote
            $trace.$trace_field = $call_record.retval
        end)
    end
end

function codegen_simulate(gen::Type{T}, args) where {T <: BasicGenFunction}
    trace_type = get_trace_type(gen)
    ir = get_ir(gen)
    stmts = Expr[]

    # initialize trace and score
    trace = gensym("trace")
    score = gensym("score")
    push!(stmts, quote
        $trace = $trace_type()
        $score = 0.
        $trace.$is_empty_field = true
    end)

    # unpack arguments
    arg_names = Symbol[arg_node.name for arg_node in ir.arg_nodes]
    push!(stmts, Expr(:(=), Expr(:tuple, arg_names...), :args))

    # record arguments in trace
    for arg_node in ir.arg_nodes
        push!(stmts, quote $trace.$(value_field(arg_node)) = $(arg_node.name) end)
    end

    # record parameters in trace
    for param in ir.params
        value_node = ir.value_nodes[param.name]
        push!(stmts, quote $trace.$(value_field(value_node)) = gen.params[$(QuoteNode(param.name))] end)
    end

    # process expression nodes in topological order
    state = BasicBlockSimulateState(trace, score, stmts)
    for node in ir.expr_nodes_sorted
        process!(ir, state, node)
    end

    if ir.output_node === nothing
        retval = :nothing
    else
        retval = quote $trace.$(value_field(something(ir.output_node))) end
    end

    push!(stmts, quote
        $trace.$call_record_field = CallRecord($score, $retval, args)
        return $trace
    end)
    Expr(:block, stmts...)
end

push!(Gen.generated_functions, quote
@generated function Gen.simulate(gen::Gen.BasicGenFunction, args)
    Gen.codegen_simulate(gen, args)
end
end)

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
# https://diffeqflux.sciml.ai/dev/examples/LV-Univ/

using DiffEqFlux, Flux, Optim, OrdinaryDiffEq, Plots
using DataInterpolations

using JLD2
using Dates

t100 = DateTime("2018-07-25T22:23:51.39")
t101 = DateTime("2018-07-25T22:23:51.721")
dt = (t101-t100).value*1e-3 |> Float32
JLD2.@load joinpath(@__DIR__, "oneOdo.jld2") jvd X100 X101

vx0, vy0 = jvd[1,3], jvd[1,4] # 0, 0
x0 = [0.0;0;vx0;vy0]
tspan = (0.0f0, dt)
tsteps = range(tspan[1],tspan[2],length=size(jvd,1))  # 0.0f0:1.0:25.0f0

#
throttle_itp = LinearInterpolation( jvd[:,1], tsteps )
sterring_itp = LinearInterpolation( jvd[:,2], tsteps )

# take in x and u
model_g = FastChain(FastDense( 6, 16, tanh),
                    FastDense(16, 16, tanh),
                    FastDense(16, 4) )
#

# The model weights are destructured into a vector of parameters
p_model = initial_params(model_g)
n_weights = length(p_model)

# # Parameters of the second equation (linear dynamics)
# p_system = Float32[0.0; 0.0, vx0, vy0]

p_all = p_model # [p_model; p_system]
θ = Float32[p_all; x0]

# dx/dt = Fx + g(x,u)
function dxdt_odo!(dx, x, p, t)
  # Destructure the parameters
  @assert length(p) == n_weights "dxdt_odo! model weight lengths do not match"
  # model_weights = p[1:n_weights]
  # state = p[(end-3):end]

  # Getting the throttle and steering inputs (from interpolators)
  thr = throttle_itp(t)
  str = sterring_itp(t)

  # Update in place
  # Velocity model dynamics and TBD odometry model
  # The neural network outputs odometry from joystick values
  dx[1:4] .= model_g([x;thr;str], p)
  dx[1:2] .+= x[3:4]
  dx
end

prob_odo = ODEProblem(dxdt_odo!, x0, tspan, p_all)
sol_odo = concrete_solve(prob_odo, Tsit5(), x0, p_all,
                          abstol = 1e-8, reltol = 1e-6)

function predict_univ(θ)
  # solving for both the ODE solution and network parameters
  x0_ = θ[(end-3):end]
  p_ = θ[1:(end-4)]
  return Array(concrete_solve(prob_odo, Tsit5(), x0_, p_,
                              saveat = tsteps))
end

function loss_univ(θ)
  pred = predict_univ(θ)
  return sum(abs2, (X101[1:2,1] - X100[1:2,1]) - pred[1:2,end]) +
         sum(abs2, pred[3:4,1] .- pred[3:4,:])
end
l = loss_univ(θ)

list_plots = []
iter = 0
callback = function (θ, l)
  global list_plots, iter

  if iter == 0
    list_plots = []
  end
  iter += 1

  println(l)

  plt = Plots.plot(predict_univ(θ), ylim = (0, 6))
  push!(list_plots, plt)
  display(plt)
  return false
end

result_univ = DiffEqFlux.sciml_train(loss_univ, θ,
                                     ADAM(0.01),
                                     cb = callback,
                                     maxiters=1000)
                                     # BFGS(initial_stepnorm = 0.01),
#

#


text:
# Principal Component Analysis

#### PCA type

type PCA
    mean::Vector{Float64}       # sample mean: of length d (mean can be empty, which indicates zero mean)
    proj::Matrix{Float64}       # projection matrix: of size d x p
    prinvars::Vector{Float64}   # principal variances: of length p
    tprinvar::Float64           # total principal variance, i.e. sum(prinvars)
    tvar::Float64               # total input variance
end

## constructor

function PCA(mean::Vector{Float64}, proj::Matrix{Float64}, pvars::Vector{Float64}, tvar::Float64)
    d, p = size(proj)
    (isempty(mean) || length(mean) == d) ||
        throw(DimensionMismatch("Dimensions of mean and proj are inconsistent."))
    length(pvars) == p ||
        throw(DimensionMismatch("Dimensions of proj and pvars are inconsistent."))
    tpvar = sum(pvars)
    tpvar <= tvar || error("principal variance cannot exceed total variance.")
    PCA(mean, proj, pvars, tpvar, tvar)
end

## properties

indim(M::PCA) = size(M.proj, 1)
outdim(M::PCA) = size(M.proj, 2)

Base.mean(M::PCA) = fullmean(indim(M), M.mean)

projection(M::PCA) = M.proj

principalvar(M::PCA, i::Integer) = M.prinvars[i]
principalvars(M::PCA) = M.prinvars

tprincipalvar(M::PCA) = M.tprinvar
tresidualvar(M::PCA) = M.tvar - M.tprinvar
tvar(M::PCA) = M.tvar

principalratio(M::PCA) = M.tprinvar / M.tvar

## use

transform{T<:Real}(M::PCA, x::AbstractVecOrMat{T}) = At_mul_B(M.proj, centralize(x, M.mean))
reconstruct{T<:Real}(M::PCA, y::AbstractVecOrMat{T}) = decentralize(M.proj * y, M.mean)

## show & dump

function show(io::IO, M::PCA)
    pr = @sprintf("%.5f", principalratio(M))
    print(io, "PCA(indim = $(indim(M)), outdim = $(outdim(M)), principalratio = $pr)")
end

function dump(io::IO, M::PCA)
    show(io, M)
    println(io)
    print(io, "principal vars: ")
    printvecln(io, M.prinvars)
    println(io, "total var = $(tvar(M))")
    println(io, "total principal var = $(tprincipalvar(M))")
    println(io, "total residual var  = $(tresidualvar(M))")
    println(io, "mean:")
    printvecln(io, mean(M))
    println(io, "projection:")
    printarrln(io, projection(M))
end

#### PCA Training

## auxiliary

const default_pca_pratio = 0.99

function check_pcaparams(d::Int, mean::Vector, md::Int, pr::Float64)
    isempty(mean) || length(mean) == d ||
        throw(DimensionMismatch("Incorrect length of mean."))
    md >= 1 || error("maxoutdim must be a positive integer.")
    0.0 < pr <= 1.0 || error("pratio must be a positive real value with pratio <= 1.0.")
end

function choose_pcadim(v::AbstractVector, ord::Vector{Int}, vsum::Float64, md::Int, pr::Float64)
    md = min(length(v), md)
    k = 1
    a = v[ord[1]]
    thres = vsum * pr
    while k < md && a < thres
        a += v[ord[k += 1]]
    end
    return k
end

## core algorithms

function pcacov(C::DenseMatrix{Float64}, mean::Vector{Float64};
                maxoutdim::Int=size(C,1),
                pratio::Float64=default_pca_pratio)

    check_pcaparams(size(C,1), mean, maxoutdim, pratio)
    Eg = eigfact!(Symmetric(copy(C)))
    ev = Eg.values
    ord = sortperm(ev; rev=true)
    vsum = sum(ev)
    k = choose_pcadim(ev, ord, vsum, maxoutdim, pratio)
    v, P = extract_kv(Eg, ord, k)
    PCA(mean, P, v, vsum)
end

function pcasvd(Z::DenseMatrix{Float64}, mean::Vector{Float64}, tw::Real;
                maxoutdim::Int=min(size(Z)...),
                pratio::Float64=default_pca_pratio)

    check_pcaparams(size(Z,1), mean, maxoutdim, pratio)
    Svd = svdfact(Z)
    v = Svd.S::Vector{Float64}
    U = Svd.U::Matrix{Float64}
    for i = 1:length(v)
        @inbounds v[i] = abs2(v[i]) / tw
    end
    ord = sortperm(v; rev=true)
    vsum = sum(v)
    k = choose_pcadim(v, ord, vsum, maxoutdim, pratio)
    si = ord[1:k]
    PCA(mean, U[:,si], v[si], vsum)
end

## interface functions

function fit(::Type{PCA}, X::DenseMatrix{Float64};
             method::Symbol=:auto,
             maxoutdim::Int=size(X,1),
             pratio::Float64=default_pca_pratio,
             mean=nothing)

    d, n = size(X)

    # choose method
    if method == :auto
        method = d < n ? :cov : :svd
    end

    # process mean
    mv = preprocess_mean(X, mean)

    # delegate to core
    if method == :cov
        C = cov(X; vardim=2, mean=isempty(mv) ? 0 : mv)::Matrix{Float64}
        M = pcacov(C, mv; maxoutdim=maxoutdim, pratio=pratio)
    elseif method == :svd
        Z = centralize(X, mv)
        M = pcasvd(Z, mv, float64(n); maxoutdim=maxoutdim, pratio=pratio)
    else
        error("Invalid method name $(method)")
    end

    return M::PCA
end


text:
import BLPDemand
import Dialysis

#try
#  include(joinpath(pkgdir(BLPDemand), "test", "runtests.jl"))
#catch e
#  println("Error in BLPDemand tests, but continuing anyway")
#end
try
  include(joinpath(pkgdir(Dialysis), "test", "runtests.jl"))
catch e
  println("Error in Dialysis tests, but continuing anyway")
end


text:
include("b_p.jl")
include("q_e.jl")
include("tau_e.jl")

include("q_95.jl")
include("beta_n.jl")

include("f_bs.jl")
include("f_cd.jl")
include("f_in.jl")

include("p_bar.jl")
include("q_dv.jl")

include("w_m.jl")

include("r_p.jl")
include("l_p.jl")

include("f_rp.jl")

include("beta_p.jl")
include("b_v.jl")

include("n_hat.jl")

include("v_loop.jl")

include("r_le.jl")
include("v_te.jl")

include("omega_ce.jl")
include("omega_pe.jl")

include("plasma_parameter.jl")

include("gray_cost.jl")


text:
abstract type MLJFluxProbabilistic <: MLJModelInterface.Probabilistic end
abstract type MLJFluxDeterministic <: MLJModelInterface.Deterministic end

const MLJFluxModel = Union{MLJFluxProbabilistic,MLJFluxDeterministic}

const doc_regressor(model_name) = """

    $model_name(; hyparameters...)

Instantiate an MLJFlux model. Available hyperparameters:

- `builder`: Default = `MLJFlux.Linear(σ=Flux.relu)` (regressors) or
   `MLJFlux.Short(n_hidden=0, dropout=0.5, σ=Flux.σ)` (classifiers)

-  `optimiser`: The optimiser to use for training. Default =
   `Flux.ADAM()`

- `loss`: The loss function used for training. Default = `Flux.mse`
   (regressors) and `Flux.crossentropy` (classifiers)

- `n_epochs`: Number of epochs to train for. Default = `10`

-  `batch_size`: The batch_size for the data. Default = 1

-  `lambda`: The regularization strength. Default = 0. Range = [0, ∞)

- `alpha`: The L2/L1 mix of regularization. Default = 0. Range = [0, 1]

-  `rng`: The random number generator (RNG) passed to builders, for
   weight intitialization, for example. Can be any `AbstractRNG` or
   the seed (integer) for a `MersenneTwister` that is reset on every
   cold restart of model (machine) training. Default =
   `GLOBAL_RNG`.

-  `acceleration`: Use `CUDALibs()` for training on GPU; default is `CPU1()`.

- `optimiser_changes_trigger_retraining`: True if fitting an
   associated machine should trigger retraining from scratch whenever
   the optimiser changes. Default = `false`

"""

doc_classifier(model_name) = doc_regressor(model_name)*"""
- `finaliser`: Operation applied to the unnormalized output of the
  final layer to obtain probabilities (outputs summing to
  one). The shape of the inputs and outputs
  of this operator must match.  Default = `Flux.softmax`.

"""

for Model in [:NeuralNetworkClassifier, :ImageClassifier]

    ex = quote
        mutable struct $Model{B,F,O,L} <: MLJFluxProbabilistic
            builder::B
            finaliser::F
            optimiser::O   # mutable struct from Flux/src/optimise/optimisers.jl
            loss::L        # can be called as in `loss(yhat, y)`
            epochs::Int    # number of epochs
            batch_size::Int  # size of a batch
            lambda::Float64  # regularization strength
            alpha::Float64   # regularizaton mix (0 for all l2, 1 for all l1)
            rng::Union{AbstractRNG,Int64}
            optimiser_changes_trigger_retraining::Bool
            acceleration::AbstractResource  # eg, `CPU1()` or `CUDALibs()`
        end

        function $Model(; builder::B   = Short()
                        , finaliser::F = Flux.softmax
                        , optimiser::O = Flux.Optimise.ADAM()
                        , loss::L      = Flux.crossentropy
                        , epochs       = 10
                        , batch_size   = 1
                        , lambda       = 0
                        , alpha        = 0
                        , rng          = Random.GLOBAL_RNG
                        , optimiser_changes_trigger_retraining = false
                        , acceleration = CPU1()
                        ) where {B,F,O,L}

            model = $Model{B,F,O,L}(builder
                                    , finaliser
                                    , optimiser
                                    , loss
                                    , epochs
                                    , batch_size
                                    , lambda
                                    , alpha
                                    , rng
                                    , optimiser_changes_trigger_retraining
                                    , acceleration
                                    )

            message = clean!(model)
            isempty(message) || @warn message

            return model
        end

        @doc doc_classifier($Model) $Model

    end
    eval(ex)

end

for Model in [:NeuralNetworkRegressor, :MultitargetNeuralNetworkRegressor]

    ex = quote
        mutable struct $Model{B,O,L} <: MLJFluxDeterministic
            builder::B
            optimiser::O  # mutable struct from Flux/src/optimise/optimisers.jl
            loss::L       # can be called as in `loss(yhat, y)`
            epochs::Int   # number of epochs
            batch_size::Int # size of a batch
            lambda::Float64 # regularization strength
            alpha::Float64  # regularizaton mix (0 for all l2, 1 for all l1)
            rng::Union{AbstractRNG,Integer}
            optimiser_changes_trigger_retraining::Bool
            acceleration::AbstractResource  # eg, `CPU1()` or `CUDALibs()`
        end

        function $Model(; builder::B   = Linear()
                        , optimiser::O = Flux.Optimise.ADAM()
                        , loss::L      = Flux.mse
                        , epochs       = 10
                        , batch_size   = 1
                        , lambda       = 0
                        , alpha        = 0
                        , rng          = Random.GLOBAL_RNG
                        , optimiser_changes_trigger_retraining=false
                        , acceleration  = CPU1()
                        ) where {B,O,L}

            model = $Model{B,O,L}(builder
                                  , optimiser
                                  , loss
                                  , epochs
                                  , batch_size
                                  , lambda
                                  , alpha
                                  , rng
                                  , optimiser_changes_trigger_retraining
                                  , acceleration)

            message = clean!(model)
            isempty(message) || @warn message

            return model
        end

        @doc $doc_regressor($Model) $Model

    end
    eval(ex)

end

const Regressor =
    Union{NeuralNetworkRegressor, MultitargetNeuralNetworkRegressor}


text:
## Automatically generated using the following and then hand-edited to
## remove duplicates and coerce constants to the appropriate type.

## using Clang.wrap_c
## context = wrap_c.init(; output_file="superlu_h.jl",
##                       header_library=x->"/usr/include/superlu",
##                       header_wrapped=(x,y)->contains(y, "superlu"),
##                       clang_diagnostics=true,
##                       clang_args=["-v"])
## context.options.wrap_structs = true
## wrap_c.wrap_c_headers(context, ["/usr/include/superlu/slu_ddefs.h"])

const EMPTY = int32(-1)
const FALSE = int32(0)
const TRUE = int32(1)
const NO_MEMTYPE = int32(4)

const NODROP = (0x0000)
const DROP_BASIC = (0x0001)
const DROP_PROWS = (0x0002)
const DROP_COLUMN = (0x0004)
const DROP_AREA = (0x0008)
const DROP_SECONDARY = (0x000E)
const DROP_DYNAMIC = (0x0010)
const DROP_INTERP = (0x0100)
const MILU_ALPHA = (1.0e-2)

typealias int_t Cint

# begin enum Stype_t
typealias Stype_t Cint   # technically, Uint32 but printing a Cint is less verbose
const SLU_NC = int32(0)
const SLU_NCP = int32(1)
const SLU_NR = int32(2)
const SLU_SC = int32(3)
const SLU_SCP = int32(4)
const SLU_SR = int32(5)
const SLU_DN = int32(6)
const SLU_NR_loc = int32(7)
# end enum Stype_t

# begin enum Dtype_t
typealias Dtype_t Cint
const SLU_S = int32(0)
const SLU_D = int32(1)
const SLU_C = int32(2)
const SLU_Z = int32(3)
# end enum Dtype_t

# begin enum Mtype_t
typealias Mtype_t Cint
const SLU_GE = int32(0)
const SLU_TRLU = int32(1)
const SLU_TRUU = int32(2)
const SLU_TRL = int32(3)
const SLU_TRU = int32(4)
const SLU_SYL = int32(5)
const SLU_SYU = int32(6)
const SLU_HEL = int32(7)
const SLU_HEU = int32(8)
# end enum Mtype_t

type SuperMatrix
    Stype::Stype_t
    Dtype::Dtype_t
    Mtype::Mtype_t
    nrow::int_t
    ncol::int_t
    Store::Ptr{None}
end

type NCformat
    nnz::int_t
    nzval::Ptr{None}
    rowind::Ptr{int_t}
    colptr::Ptr{int_t}
end

type NRformat
    nnz::int_t
    nzval::Ptr{None}
    colind::Ptr{int_t}
    rowptr::Ptr{int_t}
end

type SCformat
    nnz::int_t
    nsuper::int_t
    nzval::Ptr{None}
    nzval_colptr::Ptr{int_t}
    rowind::Ptr{int_t}
    rowind_colptr::Ptr{int_t}
    col_to_sup::Ptr{int_t}
    sup_to_col::Ptr{int_t}
end
type SCPformat
    nnz::int_t
    nsuper::int_t
    nzval::Ptr{None}
    nzval_colbeg::Ptr{int_t}
    nzval_colend::Ptr{int_t}
    rowind::Ptr{int_t}
    rowind_colbeg::Ptr{int_t}
    rowind_colend::Ptr{int_t}
    col_to_sup::Ptr{int_t}
    sup_to_colbeg::Ptr{int_t}
    sup_to_colend::Ptr{int_t}
end
type NCPformat
    nnz::int_t
    nzval::Ptr{None}
    rowind::Ptr{int_t}
    colbeg::Ptr{int_t}
    colend::Ptr{int_t}
end
type DNformat
    lda::int_t
    nzval::Ptr{None}
end
type NRformat_loc
    nnz_loc::int_t
    m_loc::int_t
    fst_row::int_t
    nzval::Ptr{None}
    rowptr::Ptr{int_t}
    colind::Ptr{int_t}
end

# begin enum yes_no_t
typealias yes_no_t Cint
const NO = int32(0)
const YES = int32(1)
# end enum yes_no_t

# begin enum fact_t
typealias fact_t Cint
const DOFACT = int32(0)
const SamePattern = int32(1)
const SamePattern_SameRowPerm = int32(2)
const FACTORED = int32(3)
# end enum fact_t

# begin enum rowperm_t
typealias rowperm_t Cint
const NOROWPERM = int32(0)
const LargeDiag = int32(1)
const MY_PERMR = int32(2)
# end enum rowperm_t

# begin enum colperm_t
typealias colperm_t Cint
const NATURAL = int32(0)
const MMD_ATA = int32(1)
const MMD_AT_PLUS_A = int32(2)
const COLAMD = int32(3)
const METIS_AT_PLUS_A = int32(4)
const PARMETIS = int32(5)
const ZOLTAN = int32(6)
const MY_PERMC = int32(7)
# end enum colperm_t

# begin enum trans_t
typealias trans_t Cint
const NOTRANS = int32(0)
const TRANS = int32(1)
const CONJ = int32(2)
# end enum trans_t

# begin enum DiagScale_t
typealias DiagScale_t Cint
const NOEQUIL = int32(0)
const ROW = int32(1)
const COL = int32(2)
const BOTH = int32(3)
# end enum DiagScale_t

# begin enum IterRefine_t
typealias IterRefine_t Cint
const NOREFINE = int32(0)
const SLU_SINGLE = int32(1)
const SLU_DOUBLE = int32(2)
const SLU_EXTRA = int32(3)
# end enum IterRefine_t

# begin enum MemType
typealias MemType Cint
const LUSUP = int32(0)
const UCOL = int32(1)
const LSUB = int32(2)
const USUB = int32(3)
const LLVL = int32(4)
const ULVL = int32(5)
# end enum MemType

# begin enum stack_end_t
typealias stack_end_t Cint
const HEAD = int32(0)
const TAIL = int32(1)
# end enum stack_end_t

# begin enum LU_space_t
typealias LU_space_t Cint
const SYSTEM = int32(0)
const USER = int32(1)
# end enum LU_space_t

# begin enum norm_t
typealias norm_t Cint
const ONE_NORM = int32(0)
const TWO_NORM = int32(1)
const INF_NORM = int32(2)
# end enum norm_t

# begin enum milu_t
typealias milu_t Cint
const SILU = int32(0)
const SMILU_1 = int32(1)
const SMILU_2 = int32(2)
const SMILU_3 = int32(3)
# end enum milu_t

# begin enum PhaseType - these are names of 1-based positions in SuperLUStat_t vectors.
typealias PhaseType Int
const COLPERM = 1
const ROWPERM = 2
const RELAX = 3
const ETREE = 4
const EQUIL = 5
const SYMBFAC = 6
const DIST = 7
const FACT = 8
const COMM = 9
const SOL_COMM = 10
const RCOND = 11
const SOLVE = 12
const REFINE = 13
const TRSV = 14
const GEMV = 15
const FERR = 16
const NPHASES = 16 # number of phases in SuperLUStat_t vectors.
# end enum PhaseType

typealias flops_t Cfloat
typealias Logical Cuchar
type superlu_options_t
    Fact::fact_t
    Equil::yes_no_t
    ColPerm::colperm_t
    Trans::trans_t
    IterRefine::IterRefine_t
    DiagPivotThresh::Cdouble
    SymmetricMode::yes_no_t
    PivotGrowth::yes_no_t
    ConditionNumber::yes_no_t
    RowPerm::rowperm_t
    ILU_DropRule::Cint
    ILU_DropTol::Cdouble
    ILU_FillFactor::Cdouble
    ILU_Norm::norm_t
    ILU_FillTol::Cdouble
    ILU_MILU::milu_t
    ILU_MILU_Dim::Cdouble
    ParSymbFact::yes_no_t
    ReplaceTinyPivot::yes_no_t
    SolveInitialized::yes_no_t
    RefineInitialized::yes_no_t
    PrintStat::yes_no_t
    nnzL::Cint
    nnzU::Cint
    num_lookaheads::Cint
    lookahead_etree::yes_no_t
    SymPattern::yes_no_t
end
type e_node
    size::Cint
    mem::Ptr{None}
end
type ExpHeader
    size::Cint
    mem::Ptr{None}
end
type LU_stack_t
    size::Cint
    used::Cint
    top1::Cint
    top2::Cint
    array::Ptr{None}
end
type SuperLUStat_t
    panel_histo::Ptr{Cint}
    utime::Ptr{Cdouble}
    ops::Ptr{flops_t}
    TinyPivots::Cint
    RefineSteps::Cint
    expansions::Cint
end
type mem_usage_t
    for_lu::Cfloat
    total_needed::Cfloat
end
type GlobalLU_t
    xsup::Ptr{Cint}
    supno::Ptr{Cint}
    lsub::Ptr{Cint}
    xlsub::Ptr{Cint}
    lusup::Ptr{Cdouble}
    xlusup::Ptr{Cint}
    ucol::Ptr{Cdouble}
    usub::Ptr{Cint}
    xusub::Ptr{Cint}
    nzlmax::Cint
    nzumax::Cint
    nzlumax::Cint
    n::Cint
    MemModel::LU_space_t
    num_expansions::Cint
    expanders::Ptr{ExpHeader}
    stack::LU_stack_t
end

const i1 = one(Cint)