text:
@testset "Basic tests" begin
    optimizer = POI.ParametricOptimizer(GLPK.Optimizer())

    MOI.set(optimizer, MOI.Silent(), true)

    x = MOI.add_variables(optimizer, 2)
    y, cy = MOI.add_constrained_variable(optimizer, POI.Parameter(0))
    z = MOI.VariableIndex(4)
    cz = MOI.ConstraintIndex{MOI.SingleVariable,POI.Parameter}(4)

    for x_i in x
        MOI.add_constraint(
            optimizer,
            MOI.SingleVariable(x_i),
            MOI.GreaterThan(0.0),
        )
    end

    @test_throws ErrorException("Cannot constrain a parameter") MOI.add_constraint(
        optimizer,
        MOI.SingleVariable(y),
        MOI.EqualTo(0.0),
    )

    @test_throws ErrorException("Variable not in the model") MOI.add_constraint(
        optimizer,
        MOI.SingleVariable(z),
        MOI.GreaterThan(0.0),
    )

    cons1 = MOI.ScalarAffineFunction(
        MOI.ScalarAffineTerm.([1.0, 1.0], [x[1], y]),
        0.0,
    )

    c1 = MOI.add_constraint(optimizer, cons1, MOI.EqualTo(2.0))

    obj_func = MOI.ScalarAffineFunction(
        MOI.ScalarAffineTerm.([1.0, 1.0], [x[1], y]),
        0.0,
    )
    MOI.set(
        optimizer,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        obj_func,
    )
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    MOI.optimize!(optimizer)

    @test MOI.get(optimizer, MOI.ObjectiveValue()) == 2

    @test MOI.get(optimizer, MOI.VariablePrimal(), x[1]) == 2

    @test_throws ErrorException("Variable not in the model") MOI.get(
        optimizer,
        MOI.VariablePrimal(),
        z,
    )

    MOI.set(optimizer, MOI.ConstraintSet(), cy, POI.Parameter(1.0))

    @test_throws ErrorException("Parameter not in the model") MOI.set(
        optimizer,
        MOI.ConstraintSet(),
        cz,
        POI.Parameter(1.0),
    )

    MOI.optimize!(optimizer)

    @test MOI.get(optimizer, MOI.ObjectiveValue()) == 2

    @test MOI.get(optimizer, MOI.VariablePrimal(), x[1]) == 1

    new_obj_func = MOI.ScalarAffineFunction(
        MOI.ScalarAffineTerm.([1.0, 1.0], [x[1], x[2]]),
        0.0,
    )
    MOI.set(
        optimizer,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        new_obj_func,
    )
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    MOI.optimize!(optimizer)

    @test MOI.get(optimizer, MOI.ObjectiveValue()) == 1

    @test MOI.get(optimizer, MOI.SolverName()) ==
          "ParametricOptimizer with GLPK attached"

    @test MOI.supports(optimizer, MOI.VariableName(), MOI.VariableIndex)
    @test MOI.supports(optimizer, MOI.ConstraintName(), MOI.ConstraintIndex)
    @test MOI.get(optimizer, MOI.ObjectiveSense()) == MOI.MIN_SENSE
    @test MOI.get(optimizer, MOI.VariableName(), x[1]) == ""
    @test MOI.get(optimizer, MOI.ConstraintName(), c1) == ""
end

@testset "Quadratic objective parameter x parameter" begin
    ipopt = Ipopt.Optimizer()
    MOI.set(ipopt, MOI.RawParameter("print_level"), 0)
    opt_in = MOIU.CachingOptimizer(MOIU.Model{Float64}(), ipopt)
    optimizer = POI.ParametricOptimizer(opt_in)

    A = [0.0 1.0; 1.0 0.0]
    a = [1.0, 1.0]

    x = MOI.add_variables(optimizer, 2)

    for x_i in x
        MOI.add_constraint(
            optimizer,
            MOI.SingleVariable(x_i),
            MOI.GreaterThan(0.0),
        )
    end

    y, cy = MOI.add_constrained_variable(optimizer, POI.Parameter(1))
    z, cz = MOI.add_constrained_variable(optimizer, POI.Parameter(1))

    quad_terms = MOI.ScalarQuadraticTerm{Float64}[]
    push!(quad_terms, MOI.ScalarQuadraticTerm(A[1, 2], y, z))

    objective_function = MOI.ScalarQuadraticFunction(
        MOI.ScalarAffineTerm.(a, x),
        quad_terms,
        0.0,
    )

    MOI.set(
        optimizer,
        MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(),
        objective_function,
    )
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    MOI.optimize!(optimizer)

    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 1.0, atol = ATOL)
    @test MOI.get(optimizer, MOI.VariablePrimal(), x[1]) == 0

    @test_throws ErrorException(
        "Cannot calculate the dual of a multiplicative parameter",
    ) MOI.get(optimizer, MOI.ConstraintDual(), cy)
    @test_throws ErrorException(
        "Cannot calculate the dual of a multiplicative parameter",
    ) MOI.get(optimizer, MOI.ConstraintDual(), cz)

    MOI.set(optimizer, MOI.ConstraintSet(), cy, POI.Parameter(2.0))
    MOI.optimize!(optimizer)
    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 2.0, atol = ATOL)

    MOI.set(optimizer, MOI.ConstraintSet(), cz, POI.Parameter(3.0))
    MOI.optimize!(optimizer)
    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 6.0, atol = ATOL)

    MOI.set(optimizer, POI.ParameterValue(), y, 5)
    MOI.set(optimizer, POI.ParameterValue(), z, 5.0)
    @test_throws ErrorException MOI.set(
        optimizer,
        POI.ParameterValue(),
        MOI.VariableIndex(10872368175),
        5.0,
    )
    MOI.optimize!(optimizer)
    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 25.0, atol = ATOL)
end

@testset "Quadratic objective parameter in affine part" begin
    ipopt = Ipopt.Optimizer()
    MOI.set(ipopt, MOI.RawParameter("print_level"), 0)
    opt_in = MOIU.CachingOptimizer(MOIU.Model{Float64}(), ipopt)
    optimizer = POI.ParametricOptimizer(opt_in)

    A = [0.0 1.0; 1.0 0.0]
    a = [2.0, 1.0]

    x = MOI.add_variables(optimizer, 2)

    for x_i in x
        MOI.add_constraint(
            optimizer,
            MOI.SingleVariable(x_i),
            MOI.GreaterThan(0.0),
        )
    end

    y, cy = MOI.add_constrained_variable(optimizer, POI.Parameter(1))
    z, cz = MOI.add_constrained_variable(optimizer, POI.Parameter(1))

    quad_terms = MOI.ScalarQuadraticTerm{Float64}[]
    push!(quad_terms, MOI.ScalarQuadraticTerm(A[1, 1], x[1], x[1]))
    push!(quad_terms, MOI.ScalarQuadraticTerm(A[1, 2], x[1], x[2]))
    push!(quad_terms, MOI.ScalarQuadraticTerm(A[2, 2], x[2], x[2]))

    objective_function = MOI.ScalarQuadraticFunction(
        MOI.ScalarAffineTerm.(a, [y, z]),
        quad_terms,
        0.0,
    )

    MOI.set(
        optimizer,
        MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(),
        objective_function,
    )
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    MOI.optimize!(optimizer)

    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 3.0, atol = ATOL)
    @test isapprox(
        MOI.get(optimizer, MOI.VariablePrimal(), x[1]),
        0,
        atol = ATOL,
    )

    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cy),
        2.0,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cz),
        1.0,
        atol = ATOL,
    )

    MOI.set(optimizer, MOI.ConstraintSet(), cy, POI.Parameter(2.0))
    MOI.optimize!(optimizer)
    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 5.0, atol = ATOL)

    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cy),
        2.0,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cz),
        1.0,
        atol = ATOL,
    )

    MOI.set(optimizer, MOI.ConstraintSet(), cz, POI.Parameter(3.0))
    MOI.optimize!(optimizer)
    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 7.0, atol = ATOL)

    MOI.set(optimizer, MOI.ConstraintSet(), cy, POI.Parameter(5.0))
    MOI.set(optimizer, MOI.ConstraintSet(), cz, POI.Parameter(5.0))
    MOI.optimize!(optimizer)
    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 15.0, atol = ATOL)
end

@testset "Quadratic constraint parameter x parameter" begin
    ipopt = Ipopt.Optimizer()
    MOI.set(ipopt, MOI.RawParameter("print_level"), 0)
    opt_in = MOIU.CachingOptimizer(MOIU.Model{Float64}(), ipopt)
    optimizer = POI.ParametricOptimizer(opt_in)

    A = [2.0 1.0; 1.0 2.0]
    a = [1.0, 1.0]

    c = [2.0, 1.0]

    x = MOI.add_variables(optimizer, 2)

    for x_i in x
        MOI.add_constraint(
            optimizer,
            MOI.SingleVariable(x_i),
            MOI.GreaterThan(0.0),
        )
    end

    MOI.add_constraint(optimizer, MOI.SingleVariable(x[1]), MOI.LessThan(20.0))

    y, cy = MOI.add_constrained_variable(optimizer, POI.Parameter(0))
    z, cz = MOI.add_constrained_variable(optimizer, POI.Parameter(0))

    quad_terms = MOI.ScalarQuadraticTerm{Float64}[]

    push!(quad_terms, MOI.ScalarQuadraticTerm(A[1, 1], y, y))
    push!(quad_terms, MOI.ScalarQuadraticTerm(A[1, 2], y, z))
    push!(quad_terms, MOI.ScalarQuadraticTerm(A[2, 2], z, z))

    constraint_function = MOI.ScalarQuadraticFunction(
        MOI.ScalarAffineTerm.(a, x),
        quad_terms,
        0.0,
    )

    MOI.add_constraint(optimizer, constraint_function, MOI.LessThan(30.0))

    obj_func =
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(c, [x[1], x[2]]), 0.0)
    MOI.set(
        optimizer,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        obj_func,
    )
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    MOI.optimize!(optimizer)

    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 50.0, atol = ATOL)
    @test isapprox(
        MOI.get(optimizer, MOI.VariablePrimal(), x[1]),
        20.0,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.VariablePrimal(), x[2]),
        10.0,
        atol = ATOL,
    )

    MOI.set(optimizer, MOI.ConstraintSet(), cy, POI.Parameter(2.0))
    MOI.optimize!(optimizer)
    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 42.0, atol = ATOL)

    MOI.set(optimizer, MOI.ConstraintSet(), cz, POI.Parameter(1.0))
    MOI.optimize!(optimizer)
    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 36.0, atol = ATOL)

    MOI.set(optimizer, MOI.ConstraintSet(), cy, POI.Parameter(-1.0))
    MOI.set(optimizer, MOI.ConstraintSet(), cz, POI.Parameter(-1.0))
    MOI.optimize!(optimizer)
    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 45.0, atol = ATOL)
end

@testset "Quadratic constraint no parameters" begin
    ipopt = Ipopt.Optimizer()
    MOI.set(ipopt, MOI.RawParameter("print_level"), 0)
    opt_in = MOIU.CachingOptimizer(MOIU.Model{Float64}(), ipopt)
    optimizer = POI.ParametricOptimizer(opt_in)

    A = [0.0 1.0; 1.0 0.0]
    a = [0.0, 0.0]
    c = [2.0, 1.0]

    x = MOI.add_variables(optimizer, 2)

    for x_i in x
        MOI.add_constraint(
            optimizer,
            MOI.SingleVariable(x_i),
            MOI.GreaterThan(1.0),
        )
        MOI.add_constraint(
            optimizer,
            MOI.SingleVariable(x_i),
            MOI.LessThan(5.0),
        )
    end

    quad_terms = MOI.ScalarQuadraticTerm{Float64}[]

    push!(quad_terms, MOI.ScalarQuadraticTerm(A[1, 2], x[1], x[2]))

    constraint_function = MOI.ScalarQuadraticFunction(
        MOI.ScalarAffineTerm.(a, x),
        [MOI.ScalarQuadraticTerm(A[1, 2], x[1], x[2])],
        0.0,
    )

    MOI.add_constraint(optimizer, constraint_function, MOI.LessThan(9.0))

    obj_func =
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(c, [x[1], x[2]]), 0.0)
    MOI.set(
        optimizer,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        obj_func,
    )
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    MOI.optimize!(optimizer)

    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 11.8, atol = ATOL)
    @test isapprox(
        MOI.get(optimizer, MOI.VariablePrimal(), x[1]),
        5.0,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.VariablePrimal(), x[2]),
        1.8,
        atol = ATOL,
    )
end

@testset "Quadratic constraint parameter x variable" begin
    ipopt = Ipopt.Optimizer()
    MOI.set(ipopt, MOI.RawParameter("print_level"), 0)
    opt_in = MOIU.CachingOptimizer(MOIU.Model{Float64}(), ipopt)
    optimizer = POI.ParametricOptimizer(opt_in)

    A = [2.0 1.0; 1.0 2.0]
    a = [1.0, 1.0]

    c = [2.0, 1.0]

    x = MOI.add_variables(optimizer, 2)

    for x_i in x
        MOI.add_constraint(
            optimizer,
            MOI.SingleVariable(x_i),
            MOI.GreaterThan(0.0),
        )
    end

    MOI.add_constraint(optimizer, MOI.SingleVariable(x[1]), MOI.LessThan(20.0))

    y, cy = MOI.add_constrained_variable(optimizer, POI.Parameter(0))

    quad_terms = MOI.ScalarQuadraticTerm{Float64}[]

    push!(quad_terms, MOI.ScalarQuadraticTerm(A[1, 1], x[1], x[1]))
    push!(quad_terms, MOI.ScalarQuadraticTerm(A[1, 2], x[1], y))
    push!(quad_terms, MOI.ScalarQuadraticTerm(A[2, 2], y, y))

    constraint_function = MOI.ScalarQuadraticFunction(
        MOI.ScalarAffineTerm.(a, x),
        quad_terms,
        0.0,
    )

    MOI.add_constraint(optimizer, constraint_function, MOI.LessThan(30.0))

    obj_func =
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(c, [x[1], x[2]]), 0.0)
    MOI.set(
        optimizer,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        obj_func,
    )
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    MOI.optimize!(optimizer)

    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 30.25, atol = ATOL)
    @test isapprox(
        MOI.get(optimizer, MOI.VariablePrimal(), x[1]),
        0.5,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.VariablePrimal(), x[2]),
        29.25,
        atol = ATOL,
    )

    MOI.set(optimizer, MOI.ConstraintSet(), cy, POI.Parameter(2.0))
    MOI.optimize!(optimizer)
    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 22.0, atol = ATOL)
end

@testset "Quadratic constraint variable x parameter" begin
    ipopt = Ipopt.Optimizer()
    MOI.set(ipopt, MOI.RawParameter("print_level"), 0)
    opt_in = MOIU.CachingOptimizer(MOIU.Model{Float64}(), ipopt)
    optimizer = POI.ParametricOptimizer(opt_in)

    A = [2.0 1.0; 1.0 2.0]
    a = [1.0, 1.0]

    c = [2.0, 1.0]

    x = MOI.add_variables(optimizer, 2)

    for x_i in x
        MOI.add_constraint(
            optimizer,
            MOI.SingleVariable(x_i),
            MOI.GreaterThan(0.0),
        )
    end

    MOI.add_constraint(optimizer, MOI.SingleVariable(x[1]), MOI.LessThan(20.0))

    y, cy = MOI.add_constrained_variable(optimizer, POI.Parameter(0))

    quad_terms = MOI.ScalarQuadraticTerm{Float64}[]

    push!(quad_terms, MOI.ScalarQuadraticTerm(A[2, 2], y, y))
    push!(quad_terms, MOI.ScalarQuadraticTerm(A[1, 2], y, x[1]))
    push!(quad_terms, MOI.ScalarQuadraticTerm(A[2, 2], x[1], x[1]))

    constraint_function = MOI.ScalarQuadraticFunction(
        MOI.ScalarAffineTerm.(a, x),
        quad_terms,
        0.0,
    )

    MOI.add_constraint(optimizer, constraint_function, MOI.LessThan(30.0))

    obj_func =
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(c, [x[1], x[2]]), 0.0)
    MOI.set(
        optimizer,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        obj_func,
    )
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    MOI.optimize!(optimizer)

    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 30.25, atol = ATOL)
    @test isapprox(
        MOI.get(optimizer, MOI.VariablePrimal(), x[1]),
        0.5,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.VariablePrimal(), x[2]),
        29.25,
        atol = ATOL,
    )

    MOI.set(optimizer, MOI.ConstraintSet(), cy, POI.Parameter(2.0))
    MOI.optimize!(optimizer)
    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 22.0, atol = ATOL)
end

@testset "Quadratic constraint variable x variable + parameter in affine part" begin
    ipopt = Ipopt.Optimizer()
    MOI.set(ipopt, MOI.RawParameter("print_level"), 0)
    opt_in = MOIU.CachingOptimizer(MOIU.Model{Float64}(), ipopt)
    optimizer = POI.ParametricOptimizer(opt_in)

    A = [2.0 1.0; 1.0 2.0]
    a = [1.0, 1.0]
    c = [2.0, 1.0]

    x = MOI.add_variables(optimizer, 2)

    for x_i in x
        MOI.add_constraint(
            optimizer,
            MOI.SingleVariable(x_i),
            MOI.GreaterThan(0.0),
        )
    end

    y, cy = MOI.add_constrained_variable(optimizer, POI.Parameter(0))

    quad_terms = MOI.ScalarQuadraticTerm{Float64}[]

    push!(quad_terms, MOI.ScalarQuadraticTerm(A[1, 1], x[1], x[1]))
    push!(quad_terms, MOI.ScalarQuadraticTerm(A[1, 2], x[1], x[2]))
    push!(quad_terms, MOI.ScalarQuadraticTerm(A[2, 2], x[2], x[2]))

    constraint_function = MOI.ScalarQuadraticFunction(
        MOI.ScalarAffineTerm.(a, [x[1], y]),
        quad_terms,
        0.0,
    )

    cons_index =
        MOI.add_constraint(optimizer, constraint_function, MOI.LessThan(25.0))

    obj_func =
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(c, [x[1], x[2]]), 0.0)
    MOI.set(
        optimizer,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        obj_func,
    )
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    MOI.optimize!(optimizer)

    @test isapprox(
        MOI.get(optimizer, MOI.ObjectiveValue()),
        9.0664,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.VariablePrimal(), x[1]),
        4.3665,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.VariablePrimal(), x[2]),
        1 / 3,
        atol = ATOL,
    )

    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cy),
        MOI.get(optimizer, MOI.ConstraintDual(), cons_index),
        atol = ATOL,
    )

    MOI.set(optimizer, MOI.ConstraintSet(), cy, POI.Parameter(2.0))
    MOI.optimize!(optimizer)
    @test isapprox(
        MOI.get(optimizer, MOI.ObjectiveValue()),
        8.6609,
        atol = ATOL,
    )
end

@testset "Quadratic constraint variable x variable + parameter in affine part - variation to assess duals" begin
    ipopt = Ipopt.Optimizer()
    MOI.set(ipopt, MOI.RawParameter("print_level"), 0)
    opt_in = MOIU.CachingOptimizer(MOIU.Model{Float64}(), ipopt)
    optimizer = POI.ParametricOptimizer(opt_in)

    A = [2.0 1.0; 1.0 2.0]
    a = [1.0, 2.0]
    c = [2.0, 1.0]

    x = MOI.add_variables(optimizer, 2)

    for x_i in x
        MOI.add_constraint(
            optimizer,
            MOI.SingleVariable(x_i),
            MOI.GreaterThan(0.0),
        )
    end

    y, cy = MOI.add_constrained_variable(optimizer, POI.Parameter(0))

    quad_terms = MOI.ScalarQuadraticTerm{Float64}[]

    push!(quad_terms, MOI.ScalarQuadraticTerm(A[1, 1], x[1], x[1]))
    push!(quad_terms, MOI.ScalarQuadraticTerm(A[1, 2], x[1], x[2]))
    push!(quad_terms, MOI.ScalarQuadraticTerm(A[2, 2], x[2], x[2]))

    constraint_function = MOI.ScalarQuadraticFunction(
        MOI.ScalarAffineTerm.(a, [x[1], y]),
        quad_terms,
        0.0,
    )

    cons_index =
        MOI.add_constraint(optimizer, constraint_function, MOI.LessThan(25.0))

    obj_func =
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(c, [x[1], x[2]]), 0.0)
    MOI.set(
        optimizer,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        obj_func,
    )
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    MOI.optimize!(optimizer)

    @test isapprox(
        MOI.get(optimizer, MOI.ObjectiveValue()),
        9.0664,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.VariablePrimal(), x[1]),
        4.3665,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.VariablePrimal(), x[2]),
        1 / 3,
        atol = ATOL,
    )

    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cy),
        2 * MOI.get(optimizer, MOI.ConstraintDual(), cons_index),
        atol = ATOL,
    )

    MOI.set(optimizer, MOI.ConstraintSet(), cy, POI.Parameter(2.0))
    MOI.optimize!(optimizer)
    @test isapprox(
        MOI.get(optimizer, MOI.ObjectiveValue()),
        8.2376,
        atol = ATOL,
    )
end

@testset "Special cases of getters" begin
    ipopt = Ipopt.Optimizer()
    MOI.set(ipopt, MOI.RawParameter("print_level"), 0)
    opt_in = MOIU.CachingOptimizer(MOIU.Model{Float64}(), ipopt)
    optimizer = POI.ParametricOptimizer(opt_in)

    A = [2.0 1.0; 1.0 2.0]
    a = [1.0, 1.0]
    c = [2.0, 1.0]

    x = MOI.add_variables(optimizer, 2)

    for x_i in x
        MOI.add_constraint(
            optimizer,
            MOI.SingleVariable(x_i),
            MOI.GreaterThan(0.0),
        )
    end

    y, cy = MOI.add_constrained_variable(optimizer, POI.Parameter(0))

    quad_terms = MOI.ScalarQuadraticTerm{Float64}[]

    push!(quad_terms, MOI.ScalarQuadraticTerm(A[1, 1], x[1], y))
    push!(quad_terms, MOI.ScalarQuadraticTerm(A[1, 2], x[1], y))
    push!(quad_terms, MOI.ScalarQuadraticTerm(A[2, 2], x[2], y))

    constraint_function = MOI.ScalarQuadraticFunction(
        MOI.ScalarAffineTerm.(a, [x[1], y]),
        quad_terms,
        0.0,
    )

    cons_index =
        MOI.add_constraint(optimizer, constraint_function, MOI.LessThan(25.0))

    obj_func = MOI.ScalarQuadraticFunction(
        MOI.ScalarAffineTerm.(c, [x[1], x[2]]),
        [MOI.ScalarQuadraticTerm(A[2, 2], x[2], y)],
        0.0,
    )
    MOI.set(
        optimizer,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        obj_func,
    )
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    @test MOI.get(optimizer, MOI.ObjectiveFunctionType()) ==
          MOI.ScalarQuadraticFunction{Float64}
    @test MOI.get(optimizer, MOI.NumberOfVariables()) == 3
end


text:
using HDF5

function word2vec(path::String)
    run(`python word2vec.py $path`)
    lines = open(readlines, "word2vec.model")
    rm("word2vec.model")
    w = String[]
    v = Float32[]
    for i = 2:length(lines)
        line = lines[i]
        items = split(line, " ")
        push!(w, items[1])
        for k = 2:length(items)
            push!(v, parse(Float32,items[k]))
        end
    end
    v = reshape(v, length(v)÷length(w), length(w))
    h5write("$(path).100d.h5", "words", w)
    h5write("$(path).100d.h5", "vectors", v)
end

word2vec(ARGS[1])

ext: jl


text:
# julia> include("/Users/ian/jd/10-19 Code/12 Julia/12.03 PE/pe4.jl")

function ispal(x::Int64)
    if sum(abs.([parse(Int64,i) for i=string(x)].-reverse([parse(Int64,i) for i=string(x)]))) == 0
        return true
    else
        return false
    end
end

pals = [k for k=[i*j for i=1:999 for j=1:i] if ispal(k)]
println("ans = $(sort(pals)[end])")

max_stars_repo_path: learning exercises/pe/pe4.jl

max_stars_repo_name: ian-double-u/NNfromScratch

max_issues_repo_path: learning exercises/pe/pe4.jl

max_issues_repo_name: ian-double-u/NNfromScratch

max_forks_repo_path: learning exercises/pe/pe4.jl

max_forks_repo_name: ian-double-u/NNfromScratch


text:
function BicubicInterpolatingFunction(arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::Vector{Vector{jdouble}}, arg3::Vector{Vector{jdouble}}, arg4::Vector{Vector{jdouble}}, arg5::Vector{Vector{jdouble}})
    return BicubicInterpolatingFunction((Vector{jdouble}, Vector{jdouble}, Vector{Vector{jdouble}}, Vector{Vector{jdouble}}, Vector{Vector{jdouble}}, Vector{Vector{jdouble}}), arg0, arg1, arg2, arg3, arg4, arg5)
end

function is_valid_point(obj::BicubicInterpolatingFunction, arg0::jdouble, arg1::jdouble)
    return jcall(obj, "isValidPoint", jboolean, (jdouble, jdouble), arg0, arg1)
end

function value(obj::BicubicInterpolatingFunction, arg0::jdouble, arg1::jdouble)
    return jcall(obj, "value", jdouble, (jdouble, jdouble), arg0, arg1)
end


text:
module TaaSSim
# package code goes here

#optimisation
using JuMP
using Clp
using GLPK
using Cbc
using LightGraphs
using LightXML

#tester functions
export build_network

#file functions
export
	readTablesFromFile, readTablesFromData, readDlmFile, readVehicleFile, readRequestFile, readNodeFile

#network functions
export
	getDistance, getTravelTime, getTravelTimes, getTravelDistance, getTravelCost, getTripPoint, getRequestedTimes, getRequestedLoads,
	getMaxTravelTimes, getCapacities,
	generateCordeauCostMatrix, generateCordeauLoads, generateCordeauTravelTimes, generateTimeWindows, generateCordeauInputs, formulateCordeauModel, solveCordeauModel

#types
export
	Vehicle, Request, Node, Location, Test

include("defs.jl")

include("file/file_io.jl")
include("file/read_darp_files.jl")

include("types/types.jl")

include("network/darp_network.jl")

end # module

max_stars_repo_name: jprendo/TaaSSim


text:
"""
    loopback(; showversion=false, idx=-1, kwargs...)

Return all "binary-encoded" arrays without change. For a positive `idx` also return
a (JSON-encoded) element of one of these arrays. Can be used for testing.
"""
function loopback(; idx=-1, kwargs...)

    bigarrs = Dict(pairs(kwargs))

    if idx < 0
        return (; bigarrs)
    end

    # otherwise return an element

    if :testarr in keys(bigarrs)
        a=kwargs[:testarr]
    else
        first_arrkey = sort([(keys(bigarrs))...])[1] # first by alphabet
        a=kwargs[first_arrkey]
    end

    elem = -1
    if ! isimage(a)
        try
            elem = a[idx...]
            if typeof(elem) in (ComplexF32, ComplexF64)
                elem = (re = real.(elem), im = imag.(elem))
            elseif elem isa Bool
                elem = UInt8(elem)
            end
        catch
        end
    end

    return (; elem, bigarrs)
end

testfunctions = (; loopback)


text:
module ConstrainedLasso

using Convex, GLMNet, ECOS, Random, LinearAlgebra, Statistics

# package code goes here
include("constrsparsereg.jl")
include("classopath.jl")
include("constrsparsereg_admm.jl")
include("genlasso.jl")

export lsq_classopath, lsq_constrsparsereg, lsq_constrsparsereg_admm, genlasso

end # module


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


text:
#! /usr/bin/env julia

path = joinpath(dirname(@__FILE__), "..", "src")

if is_windows()
    if success(`where cl.exe`)
        info("Compiling with cl.exe")
        run(`cl.exe /D_USRDLL /D_WINDLL /EHsc /Fo$(path) $(path)cclipper.cpp $(path)clipper.cpp /MT /link /DLL /OUT:$(path)cclipper.dll`)
    else
        import WinRPM
        info("Compiling with WinRPM gcc-c++")
        WinRPM.install("gcc-c++"; yes = true)
        WinRPM.install("gcc"; yes = true)
        WinRPM.install("headers"; yes = true)
        WinRPM.install("winpthreads-devel"; yes = true)

        gpp = Pkg.dir("WinRPM","deps","usr","x86_64-w64-mingw32","sys-root","mingw","bin","g++")
        RPMbindir = Pkg.dir("WinRPM","deps","usr","x86_64-w64-mingw32","sys-root","mingw","bin")
        incdir = Pkg.dir("WinRPM","deps","usr","x86_64-w64-mingw32","sys-root","mingw","include")

        push!(Base.Libdl.DL_LOAD_PATH, RPMbindir)
        ENV["PATH"] = ENV["PATH"] * ";" * RPMbindir
        if success(`$gpp --version`)
            cd(path) do
                run(`$gpp -c -fPIC -std=c++11 clipper.cpp cclipper.cpp -I $incdir`)
                run(`$gpp -shared -o cclipper.dll clipper.o cclipper.o`)
            end
        else
            error("no windows c++ compiler (cl.exe) found, and WinRPM with g++ is failing as well.")
        end
    end
end

if is_unix()
    cd(path) do
        # Note: on Mac OS X, g++ is aliased to clang++.
        run(`g++ -c -fPIC -std=c++11 clipper.cpp cclipper.cpp`)
        run(`g++ -shared -o cclipper.so clipper.o cclipper.o`)
    end
end


text:
using LinearAlgebra
using PyPlot

"""
Surrogate model. Defaults to Gaussian process.

Available radial basis function for regression include:
	:gaussian, :linear, :squared, :multiquadratic
    (:none indicates to use Gaussian process)
"""
function surrogate_model(𝒟; basis=:none, kwargs...)
    # Data
    X = map(first, 𝒟)
	y = map(last, 𝒟)

    if basis == :none # Use Gaussian processes
        kern = SE(0., 0.) # squared exponential covairance kernel
		X = cmat(X)
        gp = GP(X, y, MeanZero(), kern)
        optimize!(gp, NelderMead())
        f = x->predict_f(gp, reshape(x', (2,1)))[1][1]
    else
        # Using regression with radial basis
        if basis == :gaussian
            ψ = r->gaussian_basis(r)
        elseif basis == :linear
            ψ = r->r
        elseif basis == :squared
            ψ = r->r^2
        elseif basis == :multiquadratic
            ψ = (r;σ=1)->(r^2 + σ^2)^(1/2)
        else
            error("No regression matching $basis")
        end

        bases = radial_bases(ψ, X) # use x-values as centers

        # Predictor
        f = regression(X, y, bases)
    end

    return f::Function
end

# Kochenderfer, Mykel J.. Algorithms for Optimization (The MIT Press)

function design_matrix(X)
	n, m = length(X[1]), length(X)
	return [j==0 ? 1.0 : X[i][j] for i in 1:m, j in 0:n]
end

function linear_regression(X, y)
	θ = pinv(design_matrix(X))*y
	return x -> θ⋅[1; x]
end

function example_linear_regression()
	# Data
	D = [(1,1), (2,3), (3,3), (4,4)]
	X = map(first, D)
	y = map(last, D)

	# Predictor
	f = linear_regression(X,y)

	scatter(X, y)
	xn = 0:0.1:4
	PyPlot.plot(xn, f.(xn))
	ylim([0,5])
	xlim([0,5])
end

function regression(X, y, bases)
	B = [b(x) for x in X, b in bases]
	θ = pinv(B)*y
	return x -> sum(θ[i] * bases[i](x) for i in 1:length(θ))
end

function example_basis_regression(k=4)
	# Data
	D = [(1,1), (2,3), (3,3), (4,4)]
	X = map(first, D)
	y = map(last, D)

	bases = polynomial_bases(1, k) # 1D polynomial of degree k

	# Predictor
	f = regression(X, y, bases)

	clf()
	scatter(X, y)
	xn = 0:0.1:4
	PyPlot.plot(xn, f.(xn))
	ylim([0,5])
	xlim([0,5])
end

function example_radial_basis_regression()
	# Data
	D = [(1,1), (2,3), (3,3), (4,4)]
	X = map(first, D)
	y = map(last, D)

	ψ = r->gaussian_basis(r)
	bases = radial_bases(ψ, X) # use x-values as centers

	# Predictor
	f = regression(X, y, bases)

	clf()
	scatter(X, y)
	xn = 0:0.1:4
	PyPlot.plot(xn, f.(xn))
	ylim([0,5])
	xlim([0,5])
end

function example_radial_basis_regression_2d()
	# Data
	D = [([1,1],1), ([2,2],3), ([3,3],3), ([2,2],4)]
	X = map(first, D)
	y = map(last, D)

	ψ = r->gaussian_basis(r)
	bases = radial_bases(ψ, X) # use x-values as centers

	# Predictor
	f = regression(X, y, bases)

	clf()
	scatter(map(first, X), map(last, X), s=7, c="black", edgecolors="white", linewidths=1/2)

    r = range(0, stop=4, length=100)
    # note reverse and x, y switch (imshow has a weird mapping)
    imshow([f([x,y]) for y in reverse(r), x in r], extent=[r[1], r[end], r[1], r[end]], cmap="viridis_r")
end

function plot_regression()
	X = 0:0.1:10

	for i in 1:7
		Y = [i*x + randn() for x in X]
		PyPlot.plot(X, Y, marker=".", linestyle="none")

		f = linear_regression(X, Y)

		PyPlot.plot(X, f.(X))

		f4 = regression(X, Y, [x->x+0.8cos(5x)])
		PyPlot.plot(X, f4.(X))
	end
end

# Kochenderfer, Mykel J.. Algorithms for Optimization (The MIT Press)
function sinusoidal_bases_1d(j, k, a, b)
	T = b[j] - a[j]
	bases = Function[x->1/2]
	for i in 1:k
		push!(bases, x->sin(2π*i*x[j]/T))
		push!(bases, x->cos(2π*i*x[j]/T))
	end
	return bases
end

function sinusoidal_bases(k, a, b)
	n = length(a)
	bases = [sinusoidal_bases_1d(i, k, a, b) for i in 1:n]
	terms = Function[]
	for ks in Iterators.product([0:2k for i in 1:n]...)
		powers = [div(k+1,2) for k in ks]
		if sum(powers) ≤ k
			push!(terms,  x->prod(b[j+1](x) for (j,b) in zip(ks,bases)))
		end
	end
	return terms
end

# Kochenderfer, Mykel J.. Algorithms for Optimization (The MIT Press)
radial_bases(ψ, C, p=2) = [x->ψ(norm(x - c, p)) for c in C]

gaussian_basis(r, σ=1) = exp(-r^2 / (2σ^2))

polynomial_bases_1d(i, k) = [x->x[i]^p for p in 0:k]
function polynomial_bases(n, k)
	bases = [polynomial_bases_1d(i, k) for i in 1:n]
	terms = Function[]
	for ks in Iterators.product([0:k for i in 1:n]...)
		if sum(ks) ≤ k
			push!(terms,  x->prod(b[j+1](x) for (j,b) in zip(ks,bases)))
		end
	end
	return terms
end

bases = [(x)->x^2 * log(x)]

function regression(X, y, bases, λ)
	B = [b(x) for x in X, b in bases]
	θ = (B'B + λ*I)\B'y
	return x -> sum(θ[i] * bases[i](x) for i in 1: length(θ))
end


text:
using Flux

info("Hooking train loop")

function Flux.train!(loss, data, opt; cb = () -> ())
  loss(first(data)...)
  opt()
  cb()
end

file(x) = joinpath(@__DIR__, "..", x)

models = [
  ("MNIST MLP","mnist/mlp.jl"),
  ("MNIST Conv","mnist/conv.jl"),
  ("MNIST Autoencoder","mnist/autoencoder.jl")]

info("Testing CPU models")
for (name, p) in models
  info(name)
  include(file(p))
end
info("MNIST VAE")
info("mnist/vae.jl")

if Base.find_in_path("CuArrays") != nothing
  using CuArrays
  info("Testing GPU models")
  for (name, p) in models
    info(name)
    include(file(p))
  end
end

max_stars_repo_name: GantZA/model-zoo

max_issues_repo_name: GantZA/model-zoo

max_forks_repo_name: GantZA/model-zoo


text:
using Oceananigans.BoundaryConditions: CBC
using Oceananigans.Architectures: arch_array
using Oceananigans.Grids: halo_size

struct FieldBoundaryBuffers{W, E, S, N}
    west :: W
    east :: E
   south :: S
   north :: N
end

FieldBoundaryBuffers() = FieldBoundaryBuffers(nothing, nothing, nothing, nothing)
FieldBoundaryBuffers(grid, data, ::Missing) = nothing
FieldBoundaryBuffers(grid, data, ::Nothing) = nothing

function FieldBoundaryBuffers(grid, data, boundary_conditions)

    Hx, Hy, Hz = halo_size(grid)

    west  = create_buffer_x(architecture(grid), data, Hx, boundary_conditions.west)
    east  = create_buffer_x(architecture(grid), data, Hx, boundary_conditions.east)
    south = create_buffer_y(architecture(grid), data, Hy, boundary_conditions.south)
    north = create_buffer_y(architecture(grid), data, Hy, boundary_conditions.north)

    return FieldBoundaryBuffers(west, east, south, north)
end

create_buffer_x(arch, data, H, bc)    = nothing
create_buffer_y(arch, data, H, bc)    = nothing

create_buffer_x(arch, data, H, ::CBC) = (send = arch_array(arch, zeros(H, size(parent(data), 2), size(parent(data), 3))),
                                         recv = arch_array(arch, zeros(H, size(parent(data), 2), size(parent(data), 3))))
create_buffer_y(arch, data, H, ::CBC) = (send = arch_array(arch, zeros(size(parent(data), 1), H, size(parent(data), 3))),
                                         recv = arch_array(arch, zeros(size(parent(data), 1), H, size(parent(data), 3))))

Adapt.adapt_structure(to, buff::FieldBoundaryBuffers) =
    FieldBoundaryBuffers(Adapt.adapt(to, buff.west),
                         Adapt.adapt(to, buff.east),
                         Adapt.adapt(to, buff.north),
                         Adapt.adapt(to, buff.south))


text:
################################################# FILE DESCRIPTION #########################################################

# The Graph datatype is the core datastructure used in Graft.jl. The Graph datatype has the following fields:
# 1. nv     : The number of vertices in the graph.
# 2. ne     : The number of edges int he graph.
# 3. indxs  : The adjacency matrix for the graph. The SparseMatrixCSC type is used here, both
#             as an adjacency matrix and as an index table, that maps edges onto their entries in the
#             edge dataframe.
# 4. vdata  : A dataframe used to store vertex data. This dataframe is indexed by the internally used
#             vertex identifiers.
# 5. edata  : An edge dataframe used to store edge data. This dataframe is indexed by indxs datastructure.
# 6. lmap   : A label map that maps externally used labels onto the internally used vertex identifiers and vice versa.

################################################# IMPORT/EXPORT ############################################################
export
# Types
Graph,
# Subgraph
subgraph,
# Methods
indxs, vdata, edata, lmap

type Graph
   nv::Int
   ne::Int
   indxs::SparseMatrixCSC{Int,Int}
   vdata::AbstractDataFrame
   edata::AbstractDataFrame
   lmap::LabelMap
end

################################################# GENERATION ################################################################

include("generator.jl")

################################################# ACCESSORS #################################################################

""" The number of vertices in the graph """
nv(g::Graph) = g.nv

""" The number of edges in the graph """
ne(g::Graph) = g.ne

""" Retrieve the adjacency matrix / edge index table """
indxs(g::Graph) = g.indxs

""" Retrieve the vertex dataframe """
vdata(g::Graph) = g.vdata

""" Retrieve the edge dataframe """
edata(g::Graph) = g.edata

""" Retrieve the label map """
lmap(g::Graph)  = g.lmap

#################################################  BASICS ###################################################################

(==)(g1::Graph, g2::Graph) = nv(g1) == nv(g2) && edges(g1) == edges(g2) && lmap(g1) == lmap(g2)

Base.isequal(g1::Graph, g2::Graph) = g1 == g2 && vdata(g1) == vdata(g2) && edata(g1) == edata(g2)

function Base.copy(g::Graph)
   Graph(nv(g), ne(g), copy(indxs(g)), copy(vdata(g)), copy(edata(g)), copy(lmap(g)))
end

function Base.deepcopy(g::Graph)
   Graph(nv(g), ne(g), deepcopy(indxs(g)), deepcopy(vdata(g)), deepcopy(edata(g)), deepcopy(lmap(g)))
end

################################################# COMBINATORIAL ############################################################

include("combinatorial.jl")

################################################# MUTATION #################################################################

include("mutation.jl")

################################################# VDATA ####################################################################

include("vdata.jl")

################################################# VDATA ####################################################################

include("edata.jl")

################################################# LABELLING ################################################################

""" Set labels for all vertices in the graph """
function setlabel!(g::Graph, ls::Vector)
   if length(ls) == nv(g)
      g.lmap = setlabel!(lmap(g), ls)
      return
   else
      error("Trying to assign $(length(ls)) labels to $(nv(g)) vertices")
   end
end

""" Use a vertex property as the vertex label """
function setlabel!(g::Graph, propname::Symbol)
   ls = getvprop(g, :, propname)
   g.lmap = setlabel!(lmap(g), ls)
   return
end

""" Remove all vertex labels """
function setlabel!(g::Graph)
   g.lmap = setlabel!(lmap(g))
   return
end

""" Relabel a single vertex in the graph """
function relabel!(g::Graph, v::VertexID, l)
   g.lmap = relabel!(lmap(g), v, l)
   return
end

""" Relabel a list of vertices in the graph """
function relabel!(g::Graph, vs::VertexList, ls::AbstractVector)
   g.lmap = relabel!(lmap(g), vs, ls)
   return
end

""" Check if the input label is valid """
haslabel(g::Graph, x) = haslabel(lmap(g), x)

""" Translate the input label into the internally used vertex identifier """
decode(g::Graph, x) = decode(lmap(g), x)

""" Translate the input vertex identifier into its externally used label """
encode(g::Graph) = encode(lmap(g))
encode(g::Graph, x) = encode(lmap(g), x)

################################################# DISPLAY ##################################################################

include("display.jl")


text:
abstract type RgbdReader end

function get_timestamp_index(reader::T, timestamp::Int64) where T<:RgbdReader
    index = findfirst(x->x==timestamp, reader.timestamps)
    if isa(index, Nothing)
        return -1
    else
        return index
    end
end

function get_id_index(reader::T, id::Int64) where T<:RgbdReader
    index = findfirst(x->x==id, reader.ids)
    if isa(index, Nothing)
        return -1
    else
        return index
    end
end


text:
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()
Pkg.gc()

max_stars_repo_name: ndgnuh/RiceBPH

max_issues_repo_name: ndgnuh/RiceBPH

max_forks_repo_name: ndgnuh/RiceBPH


text:
using MarkupGenerator, Test, Random

const ASSETS_DIR = joinpath(dirname(@__FILE__), "..", "assets")

# 1. Create element without name, attributes and children
# @testset "element with no name, attributes and children" begin
#     # Test value: an element name as string
#     el = "xml"
#     # Test expected output
#     @test element() == string("<",el,"/>")
#     println(element())
#
#
#     # Test with integer value (only strings accepted)
#     el = 1
#     @test_throws MethodError element(el)
# end
#
# # 2. Create element without attributes and children
# @testset "element with custom name but no attributes,children" begin
#     # Test value: an element name as string
#     el = "svg"
#     # Test expected output
#     @test element(el) == string("<",el,"/>")
#     println(element(el))
#
#     # Test with integer value (only strings accepted)
#     el = 1
#     @test_throws MethodError element(el)
# end
#
# # 3. Create element without attributes but with children
# @testset "element with custom name, no attributes but with children" begin
#     # Test values: an element name as string and a string as a nested element
#     el = "svg"
#     children = "content"
#     # Test expected output
#     @test element(el,children) == string("<",el,">",children,"</",el,">")
#     println(element(el,children))
#
#     # Reverse test with integer values (only strings accepted)
#     el = 1
#     children = 2
#     @test_throws MethodError element(el,children)
# end
#
# # 4. Create element with attributes but without children
# @testset "element with attributes but without children" begin
#     # Test values: an element name as string and a string as a nested element
#     el = "svg"
#     key = "id"
#     value = "root"
#     attributes = Dict(key=>value)
#     # Test expected output
#     @test element(el,attributes) == string("<",el," ",key,"=\"",value,"\"/>")
#     println(element(el,attributes))
#
#     # Reverse test with integer value for attribute (Dict expected)
#     attributes = 1
#     @test_throws MethodError element(el,attributes)
#
# end
#
# # 5. Create element with attributes and children
# @testset "element with custom name, attributes and children" begin
#     # Test values: an element name as string and a string as a nested element
#     el = "svg"
#     key = "id"
#     value = "root"
#     children = "content"
#     attributes = Dict(key=>value)
#     # Test expected output
#     @test element(el,attributes,children) == string("<",el," ",key,"=\"",value,"\">",children,"</",el,">")
#     println(element(el,attributes,children))
#
#     # Reverse test with integer value for attributes and children (Dict and string expected)
#     attributes = 1
#     children = Dict()
#     @test_throws MethodError element(el,attributes,children)
#
# end
#
# # 6. Read attributes from JSON file
# @testset "Read attributes values from JSON file" begin
#     svg_json = get_json(string(PKG_ROOT_DIR,"/assets/svg/","svg.json"))
#     a = Dict("href" => "#")
#
#     @test svg_json["a"] == a
#     println(svg_json["a"])
#
#     svg_json = get_json(string(PKG_ROOT_DIR,"/assets/svg/","xyz.json"))
#     @test_throws MethodError get_json(svg_json)["a"]
#
#
# end
#
# # 7. Read attributes from JSON file and use it as values for an element
# @testset "Read attributes values from JSON file" begin
#     svg_json = get_json(string(PKG_ROOT_DIR,"/assets/svg/","svg.json"))
#     el = "defs"
#     attributes = svg_json[el]
#     result = "<defs/>"
#
#     @test element(el,attributes) == result
#     println(element(el,attributes))
#
#     attributes = false
#
#     @test_throws MethodError element(el,attributes)
#
#
# end
#
# # 8. Test generated svg elements
# @testset "Test svg-specific elements" begin
#
#     at = Dict("id" => "my rect")
#
#     println(svg_rect())
#     println(svg_rect(at))
#     println(svg_rect(Dict(),"content"))
#     println(svg_rect(at,"content"))
#
#
#
# end
#
# # 9. Test svg document function
# @testset "Test svg document function" begin
#
#     mydoc = SvgContent()
#     # mydoc.style = svg_style(Dict("id"=>"mystyles"),svg_css(CSS))
#     # mydoc.main = svg_g(Dict("id"=>"main"),svg_rect())
#
#     # println(svg_document(Dict("id"=>"root"),mydoc))
#
# end
#
# # 10. Test gradient element
# @testset "Test gradient element" begin
#
#     myrad = GradientContent()
#     myrad.stop1 = svg_stop(Dict("stop-color"=>"#ffffff"))
#
#     # println(radialGradient_template(Dict("id"=>"myrad"),myrad))
#     # println(linearGradient_template(Dict("id"=>"myrad"),myrad))
#
# end
#
# # 11. Test viewBox attribute
# @testset "Test viewBox attribute" begin
#
#     a = [0.0,0.0,100.34,200.55]
#     r = string("0.0 0.0 100.34 200.55")
#
#     @test join_str(a) == r
#     # println(join_str(a))
#
# end
#
# # 12. Test id assignments for styles
# @testset "Test create complete sample svg document" begin
#
#     #  Common parameters
#     x = 0
#     y = 0
#     width = 800
#     height = 600
#
#     my_rad_id = set_random_id()
#     myrad = GradientContent()
#     myrad.stop1 = svg_stop(Dict("stop-color"=>"#ffcc00", "offset"=>"10%"))
#     myrad.stop2 = svg_stop(Dict("stop-color"=>"#0000ff", "offset"=>"100%"))
#     rad = radialGradient_template(Dict("id"=>my_rad_id),myrad)
#
#
#     my_rect_style_id = set_random_id()
#
#     my_rect_style = CSS["rect"]
#     my_rect_style["fill"] = set_href(my_rad_id)
#     my_rect_style_el = svg_class("rect",my_rect_style,my_rect_style_id)
#     println(my_rect_style_el)
#
#     my_rect_attributes = SVG["rect"]
#     # # Delete style attribute because we use class instead
#     delete!(my_rect_attributes, "style")
#
#     my_rect_attributes["class"] = my_rect_style_id
#     my_rect_attributes["x"] = string(x)
#     my_rect_attributes["y"] = string(y)
#     my_rect_attributes["width"] = string(width)
#     my_rect_attributes["height"] = string(height)
#     my_rect = svg_rect(my_rect_attributes)
#
#
#
#     mydoc = SvgContent()
#
#     mydoc.style = svg_style(Dict(),svg_css([my_rect_style_el]))
#
#     mydoc.defs = svg_defs(Dict("id"=>"defs"),rad)
#     mydoc.main = svg_g(Dict("id"=>"main"), my_rect)
#
#     my_svg_attributes = SVG["svg"]
#     my_svg_attributes["viewBox"] = join_str([x,y,width,height])
#
#     println(svg_document(Dict("id"=>"root"),mydoc))
#
#
# end

# 13. Test doc recipe
@testset "Test doc recipe" begin

    params = get_attributes(RCP,"svg_doc_recipe")

    params["w"] = "800"
    params["h"] = "800"
    println(svg_doc_recipe(params))

end

# 14. Test canvas recipe
@testset "Test canvas recipe" begin
    params = get_attributes(RCP,"svg_canvas_recipe")

    # println(params)

    params["w"] = "800"
    params["h"] = "800"
    params["factor"] = "8"
    params["line_class"] = "no"
    params["circle_class"] = "no"
    params["rect_class"] = "yes"
    println(svg_canvas_recipe(params))
end


text:
import Documenter

Documenter.deploydocs(
    repo = "github.com/bramtayl/Keys.jl.git",
    target = "build",
    deps = nothing,
    make = nothing
)


text:
const DOC_CONSTRUCTOR =
"""
    UnivariateFinite(support,
                     probs;
                     pool=nothing,
                     augmented=false,
                     ordered=false)

Construct a discrete univariate distribution whose finite support is
the elements of the vector `support`, and whose corresponding
probabilities are elements of the vector `probs`. Alternatively,
construct an abstract *array* of `UnivariateFinite` distributions by
choosing `probs` to be an array of one higher dimension than the array
generated.

Here the word "probabilities" is an abuse of terminology as there is
no requirement that the that probabilities actually sum to one. The
only requirement is that the probabilities have a common type `T` for
which `zero(T)` is defined. In particular, `UnivariateFinite` objects
implement arbitrary non-negative, signed, or complex measures over
finite sets of labelled points. A `UnivariateDistribution` will be a
bona fide probability measure when constructed using the
`augment=true` option (see below) or when `fit` to data. And the
probabilities of a `UnivariateFinite` object `d` must be non-negative,
with a non-zero sum, for `rand(d)` to be defined and interpretable.

Unless `pool` is specified, `support` should have type
 `AbstractVector{<:CategoricalValue}` and all elements are assumed to
 share the same categorical pool, which may be larger than `support`.

*Important.* All levels of the common pool have associated
probabilities, not just those in the specified `support`. However,
these probabilities are always zero (see example below).

If `probs` is a matrix, it should have a column for each class in
`support` (or one less, if `augment=true`). More generally, `probs`
will be an array whose size is of the form `(n1, n2, ..., nk, c)`,
where `c = length(support)` (or one less, if `augment=true`) and the
constructor then returns an array of `UnivariateFinite` distributions
of size `(n1, n2, ..., nk)`.

```
using CategoricalDistributions, CategoricalArrays, Distributions
samples = categorical(['x', 'x', 'y', 'x', 'z'])
julia> Distributions.fit(UnivariateFinite, samples)
           UnivariateFinite{Multiclass{3}}
     ┌                                        ┐
   x ┤■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.6
   y ┤■■■■■■■■■■■■ 0.2
   z ┤■■■■■■■■■■■■ 0.2
     └                                        ┘

julia> d = UnivariateFinite([samples[1], samples[end]], [0.1, 0.9])
UnivariateFinite{Multiclass{3}(x=>0.1, z=>0.9)
           UnivariateFinite{Multiclass{3}}
     ┌                                        ┐
   x ┤■■■■ 0.1
   z ┤■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.9
     └                                        ┘

julia> rand(d, 3)
3-element Array{Any,1}:
 CategoricalValue{Symbol,UInt32} 'z'
 CategoricalValue{Symbol,UInt32} 'z'
 CategoricalValue{Symbol,UInt32} 'z'

julia> levels(samples)
3-element Array{Symbol,1}:
 'x'
 'y'
 'z'

julia> pdf(d, 'y')
0.0
```

### Specifying a pool

Alternatively, `support` may be a list of raw (non-categorical)
elements if `pool` is:

- some `CategoricalArray`, `CategoricalValue` or `CategoricalPool`,
  such that `support` is a subset of `levels(pool)`

- `missing`, in which case a new categorical pool is created which has
  `support` as its only levels.

In the last case, specify `ordered=true` if the pool is to be
considered ordered.

```
julia> UnivariateFinite(['x', 'z'], [0.1, 0.9], pool=missing, ordered=true)
         UnivariateFinite{OrderedFactor{2}}
     ┌                                        ┐
   x ┤■■■■ 0.1
   z ┤■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.9
     └                                        ┘

samples = categorical(['x', 'x', 'y', 'x', 'z'])
julia> d = UnivariateFinite(['x', 'z'], [0.1, 0.9], pool=samples)
     ┌                                        ┐
   x ┤■■■■ 0.1
   z ┤■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.9
     └                                        ┘

julia> pdf(d, 'y') # allowed as `'y' in levels(samples)`
0.0

v = categorical(['x', 'x', 'y', 'x', 'z', 'w'])
probs = rand(100, 3)
probs = probs ./ sum(probs, dims=2)
julia> d1 = UnivariateFinite(['x', 'y', 'z'], probs, pool=v)
100-element UnivariateFiniteVector{Multiclass{4},Symbol,UInt32,Float64}:
 UnivariateFinite{Multiclass{4}}(x=>0.194, y=>0.3, z=>0.505)
 UnivariateFinite{Multiclass{4}}(x=>0.727, y=>0.234, z=>0.0391)
 UnivariateFinite{Multiclass{4}}(x=>0.674, y=>0.00535, z=>0.321)
   ⋮
 UnivariateFinite{Multiclass{4}}(x=>0.292, y=>0.339, z=>0.369)
```

### Probability augmentation

If `augment=true` the provided array is augmented by inserting
appropriate elements *ahead* of those provided, along the last
dimension of the array. This means the user only provides probabilities
for the classes `c2, c3, ..., cn`. The class `c1` probabilities are
chosen so that each `UnivariateFinite` distribution in the returned
array is a bona fide probability distribution.

```julia
julia> UnivariateFinite([0.1, 0.2, 0.3], augment=true, pool=missing)
3-element UnivariateFiniteArray{Multiclass{2}, String, UInt8, Float64, 1}:
 UnivariateFinite{Multiclass{2}}(class_1=>0.9, class_2=>0.1)
 UnivariateFinite{Multiclass{2}}(class_1=>0.8, class_2=>0.2)
 UnivariateFinite{Multiclass{2}}(class_1=>0.7, class_2=>0.3)

d2 = UnivariateFinite(['x', 'y', 'z'], probs[:, 2:end], augment=true, pool=v)
julia> pdf(d1, levels(v)) ≈ pdf(d2, levels(v))
true
```

---

    UnivariateFinite(prob_given_class; pool=nothing, ordered=false)

Construct a discrete univariate distribution whose finite support is
the set of keys of the provided dictionary, `prob_given_class`, and
whose values specify the corresponding probabilities.

The type requirements on the keys of the dictionary are the same as
the elements of `support` given above with this exception: if
non-categorical elements (raw labels) are used as keys, then
`pool=...` must be specified and cannot be `missing`.

If the values (probabilities) are arrays instead of scalars, then an
abstract array of `UnivariateFinite` elements is created, with the
same size as the array.

"""

# # TYPES - PLAIN AND ARRAY

# R - reference type <: Unsigned
# V - type of class labels (eg, Char in `categorical(['a', 'b'])`)
# P - raw probability type
# S - scitype of samples

# Note that the keys of `prob_given_ref` need not exhaust all the
# refs of all classes but will be ordered (LittleDicts preserve order)
struct UnivariateFinite{S,V,R,P}
    scitype::Type{S}
    decoder::CategoricalDecoder{V,R}
    prob_given_ref::LittleDict{R,P,Vector{R}, Vector{P}}
end

@doc DOC_CONSTRUCTOR UnivariateFinite

"""
    UnivariateFiniteArray

Array type whose elements are `UnivariateFinite` distributions sharing
a common sample space (`CategoricalArrays` pool).

See [`UnivariateFinite`](@ref) for constructor.

"""
struct UnivariateFiniteArray{S,V,R,P,N} <:
    AbstractArray{UnivariateFinite{S,V,R,P},N}
    scitype::Type{S}
    decoder::CategoricalDecoder{V,R}
    prob_given_ref::LittleDict{R,Array{P,N},Vector{R}, Vector{Array{P,N}}}
end

const UnivariateFiniteVector{S,V,R,P} = UnivariateFiniteArray{S,V,R,P,1}

# private:
const SingletonOrArray{S,V,R,P} = Union{UnivariateFinite{S,V,R,P},
                                        UnivariateFiniteArray{S,V,R,P}}

# # CHECKS AND ERROR MESSAGES

# checks that scalar probabilities lie in [0, 1] and checks that
# vector probabilities sum to one have now been dropped, except where
# `augment=true` is specified.

# not exported:
const Prob{P} = Union{P, AbstractArray{P}} where P

# TODO: have error functions return exceptions, not throw them
# TODO: are some of these now obsolete?

const ERR_01 = DomainError("Probabilities must be in [0,1].")
err_dim(support, probs) = DimensionMismatch(
    "Probability array is incompatible "*
    "with the number of classes, $(length(support)), which should "*
    "be equal to `$(size(probs)[end])`, the last dimension "*
    "of the probability array. Perhaps you meant to set `augment=true`? ")
err_dim_augmented(support, probs) = DimensionMismatch(
    "Probability array to be augmented is incompatible "*
    "with the number of classes, $(length(support)), which should "*
    "be one more than `$(size(probs)[end])`, the last dimension "*
    "of the probability array. ")
const ERR_AUG = ArgumentError(
    "Array cannot be augmented. There are "*
    "sums along the last axis exceeding one. ")

function _check_pool(pool)
    ismissing(pool) || pool == nothing ||
        @warn "Specified pool ignored, as class labels being "*
    "generated automatically. "
    return nothing
end
_check_probs_01(probs) =
    all(0 .<= probs .<= 1) || throw(ERR_01)
_check_probs(probs) = (_check_probs_01(probs); _check_probs_sum(probs))
_check_augmentable(support, probs) = _check_probs_01(probs) &&
    size(probs)[end] + 1 == length(support) ||
    throw(err_dim_augmented(support, probs))

## AUGMENTING ARRAYS TO MAKE THEM PROBABILITY ARRAYS

_unwrap(A::Array) = A
_unwrap(A::Vector) = first(A)

isbinary(support) = length(support) == 2

# augmentation inserts the sum-subarray *before* the array:
_augment_probs(support, probs) =
    _augment_probs(Val(isbinary(support)), support, probs,)
function _augment_probs(::Val{false},
                        support,
                        probs::AbstractArray{P,N}) where {P,N}
    _check_augmentable(support, probs)
    aug_size = size(probs) |> collect
    aug_size[end] += 1
    augmentation = _unwrap(one(P) .- sum(probs, dims=N))
    all(0 .<= augmentation .<= 1) || throw(ERR_AUG)
    aug_probs = Array{P}(undef, aug_size...)
    aug_probs[fill(:, N - 1)..., 2:end] = probs
    aug_probs[fill(:, N - 1)..., 1] = augmentation
    return aug_probs
end
function _augment_probs(::Val{true},
                        support,
                        probs::AbstractArray{P,N}) where {P,N}
    _check_probs_01(probs)
    aug_size = [size(probs)..., 2]
    augmentation = one(P) .- probs
    all(0 .<= augmentation .<= 1) || throw(ERR_AUG)
    aug_probs = Array{P}(undef, aug_size...)
    aug_probs[fill(:, N)..., 2] = probs
    aug_probs[fill(:, N)..., 1] = augmentation
    return aug_probs
end

## CONSTRUCTORS - FROM DICTIONARY

# The following constructor will get called by all the others. It
# returns a UnivariateFinite object *or* a UnivariateFiniteArray,
# depending on the values of the dictionary - scalar or array - which
# represent the probabilities, one for each class in the support.
function UnivariateFinite(
    prob_given_class::AbstractDict{<:CategoricalValue, <:Prob};
    kwargs...)

    # this constructor ignores kwargs

    probs = values(prob_given_class) |> collect
#    _check_probs_01.(probs)
#    _check_probs_sum(probs)

    # `LittleDict`s preserve order of keys, which we need for rand():
    _support  = keys(prob_given_class) |> collect |> sort

    # retrieve decoder and classes from first key:
    class1         = first(_support)
    parent_decoder = decoder(class1)
    parent_classes = classes(class1)

    # TODO: throw pre-defined exception below

    issubset(_support, parent_classes) ||
        error("Categorical elements are not from the same pool. ")

    pairs = [int(c) => prob_given_class[c]
                for c in _support]

    probs1 = first(values(prob_given_class))
    S = scitype(class1) # this `scitype` defined in utilities.jl
    if  probs1 isa AbstractArray
        return UnivariateFiniteArray(S, parent_decoder, LittleDict(pairs...))
    else
        return UnivariateFinite(S, parent_decoder, LittleDict(pairs...))
    end
end

# TODO: throw and test pre-defined exceptions in this method:
function UnivariateFinite(d::AbstractDict{V,<:Prob};
                          pool=nothing,
                          ordered=false) where V

    ismissing(pool) &&
        throw(ArgumentError(
            "You cannot specify `pool=missing` "*
            "if passing `UnivariateFinite` a dictionary"))

    pool === nothing && throw(ArgumentError(
        "You must specify `pool=c` "*
        "where `c` is a "*
        "`CategoricalArray`, `CategoricalArray` or "*
        "CategoricalPool`"))

    ordered && @warn "Ignoring `ordered` key-word argument as using "*
    "specified pool to order. "

    raw_support = keys(d) |> collect
    _classes = classes(pool)
    issubset(raw_support, _classes) ||
        error("Specified support, $raw_support, not contained in "*
              "specified pool, $(levels(classes)). ")
    support = filter(_classes) do c
        c in raw_support
    end

    prob_given_class =
        LittleDict([c=>d[CategoricalArrays.DataAPI.unwrap(c)] for c in support])

    return UnivariateFinite(prob_given_class)
end

## CONSTRUCTORS - FROM ARRAYS

# example: _get_on_last(A, 4) = A[:, :, 4] if A has 3 dims:
_get_on_last(probs::AbstractArray{<:Any,N}, i) where N =
    probs[fill(:,N-1)..., i]

# 1. Univariate Finite from a vector of classes or raw labels and
# array of probs; first, a dispatcher:
function UnivariateFinite(
    support::AbstractVector,
    probs;
    kwargs...)

    if support isa AbstractArray{<:CategoricalValue}
        if :pool in keys(kwargs)
            @warn "Ignoring value of `pool` as the specified "*
            "support defines one already. "
        end
        if :ordered in keys(kwargs)
            @warn "Ignoring value of `ordered` as the "*
            "specified support defines an order already. "
        end
    end

    return _UnivariateFinite(Val(isbinary(support)),
                             support,
                             probs;
                             kwargs...)
end

# The core method, ultimately called by 1.0, 1.1, 1.2, 1.3 below.
function _UnivariateFinite(support::AbstractVector{CategoricalValue{V,R}},
                           probs::AbstractArray{P},
                           N;
                           augment=false,
                           kwargs...) where {V,R,P}

    unique(support) == support ||
        error("Non-unique vector of classes specified")

    _probs = augment ? _augment_probs(support, probs) : probs

    augment || length(support) == size(_probs) |> last ||
        throw(err_dim(support, _probs))

    # it's necessary to force the typing of the LittleDict otherwise it
    # flips to Any type (unlike regular Dict):

    if N == 0
        prob_given_class = LittleDict{CategoricalValue{V,R},P}()
    else
        prob_given_class =
            LittleDict{CategoricalValue{V,R}, AbstractArray{P,N}}()
    end
    for i in eachindex(support)
        prob_given_class[support[i]] = _get_on_last(_probs, i)
    end

    # calls dictionary constructor above:
    return UnivariateFinite(prob_given_class; kwargs...)
end

# 1.0 support does not consist of categorical elements:
function _UnivariateFinite(support,
                           probs::AbstractArray,
                           N;
                           augment=false,
                           pool=nothing,
                           ordered=false)

    # If we got here, then the vector `support` is not
    # `AbstractVector{<:CategoricalValue}`

    if pool === nothing || ismissing(pool)
        if pool === nothing
            @warn "No `CategoricalValue` found from which to extract a "*
            "complete pool of classes. "*
            "Creating a new pool (ordered=$ordered). "*
            "You can:\n"*
            " (i) specify `pool=missing` to suppress this warning; or\n"*
            " (ii) use an existing pool by specifying `pool=c` "*
            "where `c` is a "*
            "`CategoricalArray`, `CategoricalValue` or "*
            "CategoricalPool`.\n"*
            "In case (i) "*
            "specify `ordered=true` if samples are to be `OrderedFactor`. "
        end
        v = categorical(support, ordered=ordered, compress=true)
        levels!(v, support)
        _support = classes(v)
    else
        _classes = classes(pool)
        issubset(support, _classes) ||
            error("Specified support, $support, not contained in "*
                  "specified pool, $(levels(classes)). ")
        _support = filter(_classes) do c
            c in support
        end
    end

    # calls core method:
    return _UnivariateFinite(_support, probs, N;
                             augment=augment, pool=pool, ordered=ordered)
end

# 1.1 generic (non-binary) case:
_UnivariateFinite(::Val{false},
                  support::AbstractVector,
                  probs::AbstractArray{<:Any,M};
                  augment=false,
                  kwargs...) where M =
                      _UnivariateFinite(support,
                                        probs,
                                        M - 1;
                                        augment=augment,
                                        kwargs...)

# 1.2 degenerate (binary) case:
_UnivariateFinite(::Val{true},
                  support::AbstractVector,
                  probs::AbstractArray{<:Any,M};
                  augment=false,
                  kwargs...) where M =
                      _UnivariateFinite(support,
                                        probs,
                                        augment ? M : M - 1;
                                        augment=augment,
                                        kwargs...)

# 1.3 corner case, probs a scalar:
_UnivariateFinite(::Val{true},
                  support::AbstractVector,
                  probs;
                  kwargs...) =
                      UnivariateFinite(support, [probs,]; kwargs...)[1]

# 2. probablity only; unspecified support:
function UnivariateFinite(probs::AbstractArray{<:Any,N};
                          pool=nothing,
                          augment=false,
                          kwargs...) where N
    _check_pool(pool)

    # try to infer number of classes:
    if N == 1
        if augment
            c = 2
        else
            c = length(probs)
        end
    elseif N == 2
        if augment
            c = size(probs, 2) + 1
        else
            c = size(probs, 2)
        end
    else
        throw(ArgumentError(
            "You need to explicitly specify a support for "*
            "probablility arrays of three "*
            "or more dimensions. "))
    end

    support = ["class_$i" for i in 1:c]

    UnivariateFinite(support,
                     probs;
                     pool=pool,
                     augment=augment,
                     kwargs...)
end

UnivariateFiniteArray(args...; kwargs...) = UnivariateFinite(args...; kwargs...)


text:
#=
Alternating Projections methods, first try
- Julia version: 1.1.0
- Author: Oleg Soloviev
- Date: 2019-09-01 =#

"""
AlternatingProjections

Module contains types and functions for formulation and solving feasibility problem
    using projections-based methods.
"""
module AlternatingProjections

using LinearAlgebra
using FFTW
import Base.size # to add size method for FeasibleSet's subtypes

include("AbstractProblems.jl")
export Problem, Algorithm, IterativeAlgorithm, solve

export ProjectionsMethod, FeasibleSet, FeasibilityProblem, project, project!, reflect, reflect!, ConvexSet, apsolve, TwoSetsFP, TransformedSet, LinearTransformedSet,
FourierTransformedSet, forward!, backward!, generatingset

# Common type of a problem is feasibility problem, for which we need to introduce concept of
# a feasuible set
"""
    FeasibleSet

Abstract type representing a set for feasibility problem.
    For any set we should be able to take it's representative element.

    Fro any concrete subtype, we define a project operator (maybe set-valued) to be defined.

"""
abstract type FeasibleSet end

getelement(s::FeasibleSet) = error("Don't know how to take an element of $typeof(s)")

"""
Big class of feasibility problems./
"""
abstract type FeasibilityProblem <: Problem end

"""
    TwoSetsFP(A,B)

Two sets feasibility problem: for given sets A and B find x∈A∩B, if A∩B≠∅, or x∈A closest to B in some sense.

"""
struct TwoSetsFP <: FeasibilityProblem
    A::FeasibleSet
    B::FeasibleSet
end

"""
    ConvexSet

General type, no projection method is specified. For a convex set a projection is unique.

"""
abstract type ConvexSet <: FeasibleSet end

# Projections

"""
    project!(xp, x, A)

Project `x` on set `A` using preallocated `xp`.
"""
function project!(xp, x, feasset::FeasibleSet)
    error("Don't know how to project on ", typeof(feasset))
end

"""
    project!(x, A)

Project `x` on set `A` and storing result in `x`.
"""
function project!(x, feasset::FeasibleSet)
    error("Don't know how to project on ", typeof(feasset))
end

"""
    project(x, A)

Project `x` on set `A` and return the result of projection.
"""
project(x, feasset::FeasibleSet) = project!(copy(x), x, feasset)

# Reflections
"""
    reflect!(xr, x, A)

Reflect `x` on set `A` using preallocated `xp`, R(x) = 2P(x) -x, where P is projection on the set.
"""
function reflect!(xr, x, feasset::FeasibleSet)
    project!(xr,x,feasset)
    #  xr .= 2 .* xr .- x
    @. xr = 2 * xr - x
    return xr
end

reflect(x, feasset::FeasibleSet) = 2 * project(x, feasset) - x

"""
ProjectionsMethod is class of iterative algorithms for solving feasibility problems based on projections on the feasible sets.
    Examples are Alternating Projections (AP), also known in the literature as Projections on the Convex Sets (POCS),
    Douglas-Rachford algorithm (DR), their combination DRAP and many others.

"""
abstract type  ProjectionsMethod <: IterativeAlgorithm end

include("AP.jl")
include("DR.jl")
include("DRAP.jl")

include("TransformedSet.jl")

# Constraints
include("SupportConstraint.jl")
include("AmplitudeConstraint.jl")

# Iterators
include("iterators.jl")

end # module


text:
module PascalsTriangle

# Exceptions
export OutOfBoundsError, NonAdjacentError

# Types
export Entry, Row, Column, LazyColumn
export Centre, Center, LazyCentre, LazyCenter

# Accessors
export rownumber, rowposition, colnumber, value

# Checks
export isfirst, isatleft, isatright, isadjacent, areadjacent
export isinterior, issubtractable, aresubtractable, isvalid

# Conversion
export toarray

# Modifiers
export up, down, left, right, prev, next
export up!, down!, left!, right!, prev!, next!

include("docstrings.jl")
include("exceptions.jl")
include("entries.jl")
include("rows.jl")
include("columns.jl")
include("centres.jl")
include("movement.jl")

end # of module


text:
"""
    Contact{T,N}

    Abstract type containing contact information associated with Body objects.
"""
abstract type Contact{T,N} end

# constraint Jacobians
function constraint_jacobian_configuration(relative::Symbol, model::Contact,
    xp::AbstractVector, vp::AbstractVector, qp::Quaternion, ωp::AbstractVector,
    xc::AbstractVector, vc::AbstractVector, qc::Quaternion, ωc::AbstractVector,
    timestep)

    # distance Jacobian
    ∂d∂x = ∂distance∂x(relative, model.collision, xp, qp, xc, qc)
    ∂d∂q = ∂distance∂q(relative, model.collision, xp, qp, xc, qc)

    # contact point velocity Jacobian
    ∂vt∂x = ∂relative_tangential_velocity∂x(relative, model, xp, qp, vp, ωp, xc, qc, vc, ωc)
    ∂vt∂q = ∂relative_tangential_velocity∂q(relative, model, xp, qp, vp, ωp, xc, qc, vc, ωc)

    X = [
            ∂d∂x;
            szeros(1, 3);
            model.friction_parameterization * ∂vt∂x;
        ]

    Q = [
            ∂d∂q;
            szeros(1, 4);
            model.friction_parameterization * ∂vt∂q
        ]

    return [X Q]
end

function constraint_jacobian_velocity(relative::Symbol, model::Contact,
    xp::AbstractVector, vp::AbstractVector, qp::Quaternion, ωp::AbstractVector,
    xc::AbstractVector, vc::AbstractVector, qc::Quaternion, ωc::AbstractVector,
    timestep)

    # distance Jacobian
    ∂d∂x = ∂distance∂x(relative, model.collision, xp, qp, xc, qc)
    ∂d∂q = ∂distance∂q(relative, model.collision, xp, qp, xc, qc)

    # relative tangential velocity Jacobian
    ∂vt∂q = ∂relative_tangential_velocity∂q(relative, model, xp, qp, vp, ωp, xc, qc, vc, ωc)
    ∂vt∂v = ∂relative_tangential_velocity∂v(relative, model, xp, qp, vp, ωp, xc, qc, vc, ωc)
    ∂vt∂ω = ∂relative_tangential_velocity∂ω(relative, model, xp, qp, vp, ωp, xc, qc, vc, ωc)

    # recover current orientation
    if relative == :parent
        x = next_position(xp, -vp, timestep)
        ∂x∂v = linear_integrator_jacobian_velocity(x, vp, timestep)
        q = next_orientation(qp, -ωp, timestep)
        ∂q∂ω = rotational_integrator_jacobian_velocity(q, ωp, timestep);
    elseif relative == :child
        x = next_position(xc, -vc, timestep)
        ∂x∂v = linear_integrator_jacobian_velocity(x, vc, timestep)
        q = next_orientation(qc, -ωc, timestep)
        ∂q∂ω = rotational_integrator_jacobian_velocity(q, ωc, timestep);
    end

    V = [
            ∂d∂x * ∂x∂v;
            szeros(1, 3);
            model.friction_parameterization * ∂vt∂v
        ]

    Ω = [
            ∂d∂q * ∂q∂ω;
            szeros(1, 3);
            model.friction_parameterization * (∂vt∂ω + ∂vt∂q * ∂q∂ω);
        ]

    return [V Ω]
end

function impulse_map(relative::Symbol, model::Contact, pbody::Node, cbody::Node, timestep)
    # configurations
    xp, qp = next_configuration(pbody.state, timestep)
    xc, qc = next_configuration(cbody.state, timestep)

    # mapping
    X = force_mapping(relative, model, xp, qp, xc, qc)
    c = contact_point(relative, model.collision, xp, qp, xc, qc)

    if relative == :parent
        r = c - xp
        Q = rotation_matrix(inv(qp)) * skew(r) * X
    elseif relative == :child
        r = c - xc
        Q = rotation_matrix(inv(qc)) * skew(r) * X
    end

    return [
                X;
                Q;
           ]
end

function impulse_map_jacobian(relative::Symbol, jacobian::Symbol, model::Contact, pbody::Node, cbody::Node, λ, timestep)

    # configurations
    xp, qp = next_configuration(pbody.state, timestep)
    xc, qc = next_configuration(cbody.state, timestep)

    # mapping
    X = force_mapping(relative, model, xp, qp, xc, qc)

    # force Jacobian
    Xx = ∂force_mapping_jvp∂x(relative, jacobian, model, xp, qp, xc, qc, λ)
    Xq = ∂force_mapping_jvp∂q(relative, jacobian, model, xp, qp, xc, qc, λ)

    # contact point
    c = contact_point(relative, model.collision, xp, qp, xc, qc)

    # torque Jacobian
    if relative == :parent
        r = c - xp
        q = qp
    elseif relative == :child
        r = c - xc
        q = qc
    end

    Qx = rotation_matrix(inv(q)) * skew(r) * Xx
    Qx -= rotation_matrix(inv(q)) * skew(X * λ) * (∂contact_point∂x(relative, jacobian, model.collision, xp, qp, xc, qc) - (relative == jacobian ? 1.0 : 0.0) * I(3))

    Qq = rotation_matrix(inv(q)) * skew(r) * Xq
    Qq -= rotation_matrix(inv(q)) * skew(X * λ) * ∂contact_point∂q(relative, jacobian, model.collision, xp, qp, xc, qc)
    Qq += ∂rotation_matrix∂q(inv(q), skew(r) * X * λ) * Tmat()

    return [
                Xx Xq;
                Qx Qq;
           ]
end

# force mapping
function force_mapping(relative::Symbol, model::Contact,
    xp::AbstractVector, qp::Quaternion,
    xc::AbstractVector, qc::Quaternion)

    X = [
        contact_normal(model.collision, xp, qp, xc, qc)' szeros(3, 1) contact_tangent(model.collision, xp, qp, xc, qc)' * model.friction_parameterization'
    ]

    if relative == :parent
        return X
    elseif relative == :child
        return -1.0 * X
    end
end

# mapping * λ
function ∂force_mapping_jvp∂x(relative::Symbol, jacobian::Symbol,
    model::Contact,
    xp::AbstractVector, qp::Quaternion,
    xc::AbstractVector, qc::Quaternion,
    λ::AbstractVector)

    # normal
    X = ∂contact_normal_transpose∂x(jacobian, model.collision, xp, qp, xc, qc) * λ[1]

    # tangent
    λ_tangent = model.friction_parameterization' * λ[SUnitRange(3, length(λ))]
    X += ∂contact_tangent_one_transpose∂x(jacobian, model.collision, xp, qp, xc, qc) * λ_tangent[1]
    X += ∂contact_tangent_two_transpose∂x(jacobian, model.collision, xp, qp, xc, qc) * λ_tangent[2]

    if relative == :parent
        return X
    elseif relative == :child
        return -1.0 * X
    end
end

# mapping * λ
function ∂force_mapping_jvp∂q(relative::Symbol, jacobian::Symbol,
    model::Contact,
    xp::AbstractVector, qp::Quaternion,
    xc::AbstractVector, qc::Quaternion,
    λ::AbstractVector)

    # normal
    X = ∂contact_normal_transpose∂q(jacobian, model.collision, xp, qp, xc, qc) * λ[1]

    # tangent
    λ_tangent = model.friction_parameterization' * λ[SUnitRange(3, length(λ))]
    X += ∂contact_tangent_one_transpose∂q(jacobian, model.collision, xp, qp, xc, qc) * λ_tangent[1]
    X += ∂contact_tangent_two_transpose∂q(jacobian, model.collision, xp, qp, xc, qc) * λ_tangent[2]

    if relative == :parent
        return X
    elseif relative == :child
        return -1.0 * X
    end
end

# utilities
Base.length(model::Contact{T,N}) where {T,N} = N
neutral_vector(model::Contact{T,N}) where {T,N} = sones(T, Int(N / 2))
cone_degree(model::Contact{T,N}) where {T,N} = Int(N / 2)


text:
using MathProgBase
using Compat
using Compat.Test

# Here the type represents the complete instance, but it
# could also store instance data.
mutable struct HS071 <: MathProgBase.AbstractNLPEvaluator
end

# hs071
# min x1 * x4 * (x1 + x2 + x3) + x3
# st  x1 * x2 * x3 * x4 >= 25
#     x1^2 + x2^2 + x3^2 + x4^2 = 40
#     1 <= x1, x2, x3, x4 <= 5
# Start at (1,5,5,1)
# End at (1.000..., 4.743..., 3.821..., 1.379...)

function MathProgBase.initialize(d::HS071, requested_features::Vector{Symbol})
    for feat in requested_features
        if !(feat in [:Grad, :Jac, :Hess])
            error("Unsupported feature $feat")
            # TODO: implement Jac-vec and Hess-vec products
            # for solvers that need them
        end
    end
end

MathProgBase.features_available(d::HS071) = [:Grad, :Jac, :Hess]

MathProgBase.eval_f(d::HS071, x) = x[1] * x[4] * (x[1] + x[2] + x[3]) + x[3]

function MathProgBase.eval_g(d::HS071, g, x)
    g[1] = x[1]   * x[2]   * x[3]   * x[4]
    g[2] = x[1]^2 + x[2]^2 + x[3]^2 + x[4]^2
end

function MathProgBase.eval_grad_f(d::HS071, grad_f, x)
    grad_f[1] = x[1] * x[4] + x[4] * (x[1] + x[2] + x[3])
    grad_f[2] = x[1] * x[4]
    grad_f[3] = x[1] * x[4] + 1
    grad_f[4] = x[1] * (x[1] + x[2] + x[3])
end

MathProgBase.jac_structure(d::HS071) = [1,1,1,1,2,2,2,2],[1,2,3,4,1,2,3,4]
# lower triangle only
MathProgBase.hesslag_structure(d::HS071) = [1,2,2,3,3,3,4,4,4,4],[1,1,2,1,2,3,1,2,3,4]

function MathProgBase.eval_jac_g(d::HS071, J, x)
    # Constraint (row) 1
    J[1] = x[2]*x[3]*x[4]  # 1,1
    J[2] = x[1]*x[3]*x[4]  # 1,2
    J[3] = x[1]*x[2]*x[4]  # 1,3
    J[4] = x[1]*x[2]*x[3]  # 1,4
    # Constraint (row) 2
    J[5] = 2*x[1]  # 2,1
    J[6] = 2*x[2]  # 2,2
    J[7] = 2*x[3]  # 2,3
    J[8] = 2*x[4]  # 2,4
end

function MathProgBase.eval_hesslag(d::HS071, H, x, σ, μ)
    # Again, only lower left triangle
    # Objective
    H[1] = σ * (2*x[4])               # 1,1
    H[2] = σ * (  x[4])               # 2,1
    H[3] = 0                          # 2,2
    H[4] = σ * (  x[4])               # 3,1
    H[5] = 0                          # 3,2
    H[6] = 0                          # 3,3
    H[7] = σ* (2*x[1] + x[2] + x[3])  # 4,1
    H[8] = σ * (  x[1])               # 4,2
    H[9] = σ * (  x[1])               # 4,3
    H[10] = 0                         # 4,4

    # First constraint
    H[2] += μ[1] * (x[3] * x[4])  # 2,1
    H[4] += μ[1] * (x[2] * x[4])  # 3,1
    H[5] += μ[1] * (x[1] * x[4])  # 3,2
    H[7] += μ[1] * (x[2] * x[3])  # 4,1
    H[8] += μ[1] * (x[1] * x[3])  # 4,2
    H[9] += μ[1] * (x[1] * x[2])  # 4,3

    # Second constraint
    H[1]  += μ[2] * 2  # 1,1
    H[3]  += μ[2] * 2  # 2,2
    H[6]  += μ[2] * 2  # 3,3
    H[10] += μ[2] * 2  # 4,4

end

function nlptest(solver)
    @testset "Testing NLP with $solver" begin
        m = MathProgBase.NonlinearModel(solver)
        l = [1,1,1,1]
        u = [5,5,5,5]
        lb = [25, 40]
        ub = [Inf, 40]
        MathProgBase.loadproblem!(m, 4, 2, l, u, lb, ub, :Min, HS071())
        MathProgBase.setwarmstart!(m,[1,5,5,1])

        MathProgBase.optimize!(m)
        stat = MathProgBase.status(m)

        @test stat == :Optimal
        x = MathProgBase.getsolution(m)
        @test isapprox(x[1], 1.0000000000000000, atol=1e-5)
        @test isapprox(x[2], 4.7429996418092970, atol=1e-5)
        @test isapprox(x[3], 3.8211499817883077, atol=1e-5)
        @test isapprox(x[4], 1.3794082897556983, atol=1e-5)
        @test isapprox(MathProgBase.getobjval(m), 17.014017145179164, atol=1e-5)

        # Test that a second call to optimize! works
        MathProgBase.setwarmstart!(m,[1,5,5,1])
        MathProgBase.optimize!(m)
        stat = MathProgBase.status(m)
        @test stat == :Optimal
        if hasmethod(MathProgBase.freemodel!, Tuple{typeof(m)})
            MathProgBase.freemodel!(m)
        end
    end
end

# Same as above but no hessian callback
mutable struct HS071_2 <: MathProgBase.AbstractNLPEvaluator
end

# hs071
# min x1 * x4 * (x1 + x2 + x3) + x3
# st  x1 * x2 * x3 * x4 >= 25
#     x1^2 + x2^2 + x3^2 + x4^2 = 40
#     1 <= x1, x2, x3, x4 <= 5
# Start at (1,5,5,1)
# End at (1.000..., 4.743..., 3.821..., 1.379...)

function MathProgBase.initialize(d::HS071_2, requested_features::Vector{Symbol})
    for feat in requested_features
        if !(feat in [:Grad, :Jac])
            error("Unsupported feature $feat")
            # TODO: implement Jac-vec and Hess-vec products
            # for solvers that need them
        end
    end
end

MathProgBase.features_available(d::HS071_2) = [:Grad, :Jac]

MathProgBase.eval_f(d::HS071_2, x) = MathProgBase.eval_f(HS071(), x)

MathProgBase.eval_g(d::HS071_2, g, x) = MathProgBase.eval_g(HS071(), g, x)

MathProgBase.eval_grad_f(d::HS071_2, grad_f, x) = MathProgBase.eval_grad_f(HS071(), grad_f, x)

MathProgBase.jac_structure(d::HS071_2) = MathProgBase.jac_structure(HS071())

MathProgBase.eval_jac_g(d::HS071_2, J, x) = MathProgBase.eval_jac_g(HS071(), J, x)

function nlptest_nohessian(solver)
    @testset "Testing NLP without Hessian with $solver" begin
        m = MathProgBase.NonlinearModel(solver)
        l = [1,1,1,1]
        u = [5,5,5,5]
        lb = [25, 40]
        ub = [Inf, 40]
        MathProgBase.loadproblem!(m, 4, 2, l, u, lb, ub, :Min, HS071_2())
        MathProgBase.setwarmstart!(m,[1,5,5,1])

        MathProgBase.optimize!(m)
        stat = MathProgBase.status(m)

        @test stat == :Optimal
        x = MathProgBase.getsolution(m)
        @test isapprox(x[1], 1.0000000000000000, atol=1e-5)
        @test isapprox(x[2], 4.7429996418092970, atol=1e-5)
        @test isapprox(x[3], 3.8211499817883077, atol=1e-5)
        @test isapprox(x[4], 1.3794082897556983, atol=1e-5)
        @test isapprox(MathProgBase.getobjval(m), 17.014017145179164, atol=1e-5)

        # Test that a second call to optimize! works
        MathProgBase.setwarmstart!(m,[1,5,5,1])
        MathProgBase.optimize!(m)
        stat = MathProgBase.status(m)
        @test stat == :Optimal
        if hasmethod(MathProgBase.freemodel!, Tuple{typeof(m)})
            MathProgBase.freemodel!(m)
        end
    end
end

# a test for convex nonlinear solvers

mutable struct QCQP <: MathProgBase.AbstractNLPEvaluator
end

# min x - y
# st  x + x^2 + x*y + y^2 <= 1
#     -2 <= x, y <= 2
# solution: x+y = -1/3
# optimal objective -1-4/sqrt(3)

function MathProgBase.initialize(d::QCQP, requested_features::Vector{Symbol})
    for feat in requested_features
        if !(feat in [:Grad, :Jac, :Hess])
            error("Unsupported feature $feat")
            # TODO: implement Jac-vec and Hess-vec products
            # for solvers that need them
        end
    end
end

MathProgBase.features_available(d::QCQP) = [:Grad, :Jac, :Hess]

MathProgBase.eval_f(d::QCQP, x) = x[1]-x[2]

function MathProgBase.eval_g(d::QCQP, g, x)
    g[1] = x[1] + x[1]^2 + x[1]*x[2] + x[2]^2
end

function MathProgBase.eval_grad_f(d::QCQP, grad_f, x)
    grad_f[1] = 1
    grad_f[2] = -1
end

MathProgBase.jac_structure(d::QCQP) = [1,1],[1,2]
# lower triangle only
MathProgBase.hesslag_structure(d::QCQP) = [1,2,2],[1,1,2]

function MathProgBase.eval_jac_g(d::QCQP, J, x)
    J[1] = 1 + 2x[1] + x[2]
    J[2] = x[1]+2x[2]
end

function MathProgBase.eval_hesslag(d::QCQP, H, x, σ, μ)
    # Again, only lower left triangle
    # Objective, linear

    # First constraint
    H[1] = 2μ[1] # d/dx^2
    H[2] = μ[1]  # d/dxdy
    H[3] = 2μ[1] # d/dy^2

end

function convexnlptest(solver)
    @testset "Testing convex NLP with $solver" begin
        m = MathProgBase.NonlinearModel(solver)
        l = [-2,-2]
        u = [2,2]
        lb = [-Inf]
        ub = [1.0]
        MathProgBase.loadproblem!(m, 2, 1, l, u, lb, ub, :Min, QCQP())

        MathProgBase.optimize!(m)
        stat = MathProgBase.status(m)

        @test stat == :Optimal
        x = MathProgBase.getsolution(m)
        @test isapprox(x[1]+x[2], -1/3, atol=1e-3)
        @test isapprox(MathProgBase.getobjval(m), -1-4/sqrt(3), atol=1e-5)

        # Test that a second call to optimize! works
        MathProgBase.optimize!(m)
        stat = MathProgBase.status(m)
        @test stat == :Optimal
        if hasmethod(MathProgBase.freemodel!, Tuple{typeof(m)})
            MathProgBase.freemodel!(m)
        end
    end
end

# a test for unconstrained nonlinear solvers

mutable struct Rosenbrock <: MathProgBase.AbstractNLPEvaluator
end

# min (1-x)^2 + 100*(y-x^2)^2
# solution: (x,y) = (1,1)
# optimal objective 0

function MathProgBase.initialize(d::Rosenbrock, requested_features::Vector{Symbol})
    for feat in requested_features
        if !(feat in [:Grad, :Hess])
            error("Unsupported feature $feat")
            # TODO: implement Jac-vec and Hess-vec products
            # for solvers that need them
        end
    end
end

MathProgBase.features_available(d::Rosenbrock) = [:Grad, :Hess]

MathProgBase.eval_f(d::Rosenbrock, x) = (1-x[1])^2 + 100*(x[2]-x[1]^2)^2

MathProgBase.eval_g(d::Rosenbrock, g, x) = nothing

function MathProgBase.eval_grad_f(d::Rosenbrock, grad_f, x)
    grad_f[1] = -2*(1-x[1]) - 400*x[1]*(x[2]-x[1]^2)
    grad_f[2] = 200*(x[2]-x[1]^2)
end

MathProgBase.hesslag_structure(d::Rosenbrock) = [1,2,2],[1,1,2]

function MathProgBase.eval_hesslag(d::Rosenbrock, H, x, σ, μ)
    # Again, only lower left triangle

    H[1] = σ*(2+100*(8x[1]^2-4*(x[2]-x[1]^2))) # d/dx^2
    H[2] = -σ*400*x[1]  # d/dxdy
    H[3] = σ*200 # d/dy^2

end

function rosenbrocktest(solver)
    @testset "Testing NLP on the Rosenbrock function with $solver" begin
        m = MathProgBase.NonlinearModel(solver)
        l = [-Inf,-Inf]
        u = [Inf,Inf]
        lb = Float64[]
        ub = Float64[]
        MathProgBase.loadproblem!(m, 2, 0, l, u, lb, ub, :Min, Rosenbrock())

        MathProgBase.setwarmstart!(m,[10.0,10.0])
        MathProgBase.optimize!(m)
        stat = MathProgBase.status(m)

        @test stat == :Optimal
        x = MathProgBase.getsolution(m)
        @test isapprox(x[1], 1.0, atol=1e-5)
        @test isapprox(x[2], 1.0, atol=1e-5)
        @test isapprox(MathProgBase.getobjval(m), 0.0, atol=1e-5)

        # Test that a second call to optimize! works
        MathProgBase.optimize!(m)
        stat = MathProgBase.status(m)
        @test stat == :Optimal
        if hasmethod(MathProgBase.freemodel!, Tuple{typeof(m)})
            MathProgBase.freemodel!(m)
        end
    end
end


text:
module IO

export outMCMC
export summaryMCMC

using DelimitedFiles
using CSV
using StatsBase
using MCMCChains
using StatsPlots

macro name(arg)
    x = string(arg)
    quote
        $x
    end
end

function outMCMC(folder::String,thisVar,output)
        out0 = open(folder*"/$(thisVar)Out", "a")
        writedlm(out0, output)
        close(out0)
end

function summaryMCMC(param;summary=false,outFolder=pwd()*"/outMCMC")
        param = CSV.read("$outFolder/$(param)Out",CSV.Tables.matrix,header=false)
                if summary==true
                        chn = Chains(param)
                        display(chn)
                        display(plot(chn))
                        param = mean(Matrix(param),dims=1)
                else param = mean(Matrix(param),dims=1)
                end
        return param
end

end


text:
const integral_offset = 1e-8

const integral_zero = 0.0 + integral_offset
const integral_one = 1.0 - integral_offset

const min_I_P = 0
const max_I_P = 40
const no_pts_I_P = 301

const min_eta_CD = 0.01
const max_eta_CD = 1.10
const no_pts_eta_CD = 31

const min_rho = 0.0 + 1e-8
const max_rho = 1.0 - 1e-8
const no_pts_rho = 21

const min_T_bar = 3.0
const max_T_bar = 36.0
const no_pts_T_bar = 51

const max_roots = 10

const __is_logging = true

const H_mode_scaling = Dict(
  :constant => 0.1445,
  :I => 0.93,
  :R => 1.39,
  :a => 0.58,
  :kappa => 0.78,
  :n => 0.41,
  :B => 0.15,
  :A => 0.19,
  :P => 0.69
)

const L_mode_scaling = Dict(
  :constant => 0.0578,
  :I => 0.96,
  :R => 1.89,
  :a => -0.06,
  :kappa => 0.64,
  :n => 0.40,
  :B => 0.03,
  :A => 0.20,
  :P => 0.73
)

const scaling_keys = keys(H_mode_scaling)

const E_F = 17.6
const E_Li = 4.8

const reduced_mass = 1.124656e6

const gamov_const = 34.3827

const bosch_hale_coeffs = [
  +1.17302e-9,
  +1.51361e-2,
  +7.51886e-2,
  +4.60643e-3,
  +1.35000e-2,
  -1.06750e-4,
  +1.36600e-5
]

const secondary_limits = OrderedDict(
  :beta => "TB",
  :kink => "KF",
  :pcap => "PC",
  :wall => "WL",
  :heat => "DV"
)

const secondary_params = OrderedDict(
  :beta => "beta_N",
  :kink => "q_95",
  :pcap => "P_E",
  :wall => "P_W",
  :heat => "q_DV"
)

const wave_gamma = 8.562

const wave_M = 1.0

const wave_n = 0.77

const wave_m = 1.38

const wave_x_r = 2.47

const wave_c = 0.778


text:
function R_bench(N=10_000_000, K = 100; libpath = "")
  R"""
  if ($libpath != "") {
    .libPaths($libpath)
  }

  # set higher memory limit on Windows
  memory.limit(2^31-1)

  library(data.table)
  library(pipeR)
  library(magrittr)

  source("src/r/1_generate_synthetic.r")
  timings <- rbench(N = $N, K = $K)
  gc()
  """
  @rget timings
  return timings
end


text:
module AdditionalUnits

# using pkgs
using Unitful

# pkg constants
const _UF = Unitful

# dimensions
@dimension          𝐂           "𝐂"         Currency

# units
## currency
@refunit            €           "€"         Euro        𝐂           true
## electricity
@unit               VA          "VA"        VoltAmpere  1.0u"W"     true
@unit               Wh          "Wh"        WattHour    3600.0u"J"  true
## mass
@unit               ton         "ton"       Ton         1000.0u"kg" true

# init function to make this work at run-time
const localunits = _UF.basefactors
const localpromotion = _UF.promotion
function __init__()
    merge!(_UF.basefactors, localunits)
    merge!(_UF.promotion, localpromotion)
    _UF.register(AdditionalUnits)
end

end


text:
# ? activate Splines
import Pkg
Pkg.activate("Splines")

using PlotlyJS
using Revise

import Splines: BSpline, ⊗
import Splines.Utils: x, y, z

"""
    Fitting splines in R³ to create interpolations in RGB color space.
"""

Revise.revise()

t = 1:6*π
nodes = Array{Float64}(undef, 3, length(t))
nodes[1, :] = (sin.(t) .+ 1) .* 255/2
nodes[2, :] = (cos.(t) .+ 1) .* 255/2
nodes[3, :] = Array(range(0, stop=255, length=size(t, 1)))

r_nodes = [[0, 0, 0] [150, 0, 0] [255, 0, 0] [255, 0, 0]]
b_nodes = [[0, 0, 0] [0, 150, 0] [0, 255, 0] [0, 255, 0]]
g_nodes = [[0, 0, 0] [0, 0, 150] [0, 0, 255] [0, 0, 255]]

lines = []

# plot RGB lines
for points in [r_nodes, b_nodes, g_nodes]

    spline = BSpline(points, 1; δt=0.005)

    plotted_line = scatter(
        x=spline.coordinates[1, :],
        y=spline.coordinates[2, :],
        z=spline.coordinates[3, :],
        mode="lines",
        type="scatter3d",
        line = attr(color=spline.coordinates, width=5), name=nothing)
    push!(lines, plotted_line)
end

# plot splines of different degrees

for degree in (1, 3, 5, 7,)
    spline = BSpline(nodes, degree; δt=0.005)
    plotted_line = scatter(
        x=spline.coordinates[1, :],
        y=spline.coordinates[2, :],
        z=spline.coordinates[3, :],
        mode="lines",
        type="scatter3d",
        line = attr(color=spline.coordinates, width=20), name=nothing)
    push!(lines, plotted_line)
end

# plot
layout = Layout(
    showlegend=false,
    scene=attr(
            xaxis_title="R",
            yaxis_title="G",
            zaxis_title="B",
        xaxis=attr(
            nticks=4,
            range=[0, 255],
            axis=([], false),
            name="R",
        ),
        yaxis=attr(
            nticks=4,
            range=[0, 255],
            axis=([], false),
        ),
        zaxis=attr(
            nticks=4,
            range=[0, 255],
            axis=([], false),
        ),
    ),
)

display(plot([lines...], layout))
@info "Done!"

max_stars_repo_name: FedeClaudi/Splines

max_issues_repo_name: FedeClaudi/Splines

max_forks_repo_name: FedeClaudi/Splines


text:
function addSummaryStat!(df::DataFrame, x::Symbol, y::Symbol; order=:none, func = mean,
    colors=[:black], kwaargs...)

	if length(colors) == 1
       colors = repeat(colors, length(unique(df[!,x])))
    end

    if func == mean
        temp = sort!(combine(groupby(df, x), y => func => :mean, y=> std => :sd))

		if order != :none
			customSort!(temp, [(:construct,order)])
		end

        for i in 1:size(temp)[1]
            plot!([i-0.1, i-0.1],
                [temp.mean[i] - temp.sd[i], temp.mean[i] + temp.sd[i]],
                linewidth=3, color=colors[i]; kwaargs...)
            scatter!([i-0.1, i-0.1],
                [temp.mean[i], temp.mean[i]],
                markersize=2, marker=:square,
                color=:white,
                markerstrokewidth=2,
                markerstrokecolor=:white)
        end

    elseif func == median
        temp = sort!(combine(groupby(df, x),
            y => func => :median,
            y => (x -> percentile(x,25)) => :ybot,
            y => (x -> percentile(x,75)) => :ytop))

        for i in 1:size(temp)[1]
            plot!([i-0.1, i-0.1],
                [temp.ybot[i], temp.ytop[i]],
                linewidth=3;kwaargs...)
            scatter!([i-0.1, i-0.1],
                [temp.median[i], temp.median[i]],
                markersize=2, marker=:square,
                color=:white,
                markerstrokewidth=2,
                markerstrokecolor=:white)
        end
    end
end


text:
module PrecompileFuncs

using OptimDash
using OptimDash.Optim
using OptimDash.PlotlyJS

function precompile_optim()
    chn = Base.Channel{Vector{Float64}}(Inf)
    user_funcs = "f(x) = x[1]^2 + 3*x[2]^2"
    funcs_init = "[3.0, 3.0]"
    lb_str = "[-10.0, -10.0]"
    ub_str = "[10.0, 10.0]"
    eval_f, init_f, lb_f, ub_f = OptimDash.eval_definition(
        user_funcs, funcs_init, lb_str, ub_str
    )
    function g(f)
        function _g(x)
            put!(chn, x)
            return Base.invokelatest(f, x)
        end
        return _g
    end
    alg = BFGS()
    results = optimize(g(eval_f), init_f, alg, Optim.Options(x_tol = 1e-5))

    bound_range = (ub_f - lb_f)./1000
    range_value = [collect(1:10) for _ in 1:2]
    z =  rand(1:10,10,10)
    heat_plt = heatmap(
        x = range_value[1], y = range_value[2],
        z = z
    )
    # scatter_plt = scatter(
    #     ;x=x_graph, y=y_graph,
    #     mode="lines", line_color="black"
    # )
    fig = OptimDash.PlotlyJS.plot(heat_plt)
    return "success"
end
end


text:
using Documenter
using Ghost

makedocs(
    sitename = "Ghost",
    format = Documenter.HTML(),
    modules = [Ghost],
    pages = [
        "Main" => "index.md",
        "Linearized traces" => "trace.md",
        "Tape anatomy" => "tape.md",
        "Loops" => "loops.md",
        "Reference" => "reference.md",
    ],
)

deploydocs(
    repo = "github.com/dfdx/Ghost.jl.git",
    devbranch = "main",
)


text:
# f.jl
# F-tests
#
# Copyright (C) 2019   Raphael Saavedra
#
# 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.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

export VarianceFTest

struct VarianceFTest <: HypothesisTest
    n_x::Int   # size of sample 1
    n_y::Int   # size of sample 2
    df_x::Int  # degrees of freedom 1
    df_y::Int  # degrees of freedom 2
    F::Real    # test statistic
end

"""
    VarianceFTest(x::AbstractVector{<:Real}, y::AbstractVector{<:Real})

Perform an F-test of the null hypothesis that two real-valued vectors `x` and `y` have equal variances.

Implements: [`pvalue`](@ref)

# References

  * George E. P. Box, "Non-Normality and Tests on Variances", Biometrika 40 (3/4): 318–335, 1953.

# External links

  * [F-test of equality of variances on Wikipedia](https://en.wikipedia.org/wiki/F-test_of_equality_of_variances)
"""
function VarianceFTest(x::AbstractVector{<:Real}, y::AbstractVector{<:Real})
    n1, n2 = length(x), length(y)
    F = var(x) / var(y)
    return VarianceFTest(n1, n2, n1-1, n2-1, F)
end

testname(::VarianceFTest) = "Variance F-test"
population_param_of_interest(x::VarianceFTest) = ("variance ratio", 1.0, x.F)
default_tail(test::VarianceFTest) = :both

function show_params(io::IO, x::VarianceFTest, ident)
    println(io, ident, "number of observations: [$(x.n_x), $(x.n_y)]")
    println(io, ident, "F statistic:            $(x.F)")
    println(io, ident, "degrees of freedom:     [$(x.df_x), $(x.df_y)]")
end

function pvalue(x::VarianceFTest; tail=:both)
    dist = FDist(x.df_x, x.df_y)
    if tail == :both
        return 1 - 2*abs(cdf(dist, x.F) - 0.5)
    elseif tail == :right
        return 1 - cdf(dist, x.F)
    elseif tail == :left
        return cdf(dist, x.F)
    else
        throw(ArgumentError("tail=$(tail) is invalid"))
    end
end


text:
module MainComp

export main

include("Portran.jl")
include("Styles.jl")

using JSON

function main()
    global file_to_parse = ARGS[1]
    file_name = file_to_parse[1:end-4]
    io = read(file_to_parse, String)

    try
        tks = tokenise(io)
        # println(join(tks.tokens,"\n"))
        env = syntactic_parse(tks)
        write("$file_name.json", string(json(env, 4)))
        transpile(tks, env, file_name*".c")
        run(`gcc $file_name.c -o $file_name.exe`)
        run(`./$file_name.exe`)
    catch
    end
end # function main

end # Module

using .MainComp
main()

max_stars_repo_name: GabrielTrettel/Portran

max_issues_repo_name: GabrielTrettel/Portran

max_forks_repo_name: GabrielTrettel/Portran