text:
include("basic/gamsworld.jl")

@testset "basic tests" begin

@testset "no objective and start value" begin
    println("==================================")
    println("No objective and start value")
    println("==================================")
    juniper = DefaultTestSolver(log_levels=[:Table])

    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        juniper...)
    )

    @variable(m, x, Int, start=3)

    @constraint(m, x >= 0)
    @constraint(m, x <= 5)
    @NLconstraint(m, x^2 >= 17)

    status = solve(m)
    rand_num = rand()
    status = solve(m)
    @test rand() != rand_num

    inner = internalmodel(m)
    @test JuMP.termination_status(m) == MOI.LOCALLY_SOLVED
    @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
    @test JuMP.dual_status(m) == MOI.FEASIBLE_POINT
    @test status == MOI.LOCALLY_SOLVED
    @test isapprox(JuMP.value(x), 5, atol=sol_atol)
    @test Juniper.getnsolutions(inner) == 1
    @test inner.primal_start[1] == 3
end

@testset "bruteforce" begin
    println("==================================")
    println("Bruteforce")
    println("==================================")

    m = Model()

    @variable(m, 1 <= x[1:4] <= 5, Int)

    special_minimizer_fct(x) = x
    grad(x) = 1.0
    grad2(x) = 0.0

    register_args = [:special_minimizer_fct, 1, special_minimizer_fct, grad, grad2]
    JuMP.register(m, register_args...)

    @NLobjective(m,Min,special_minimizer_fct(x[1]))

    juniper_all_solutions = DefaultTestSolver(
        branch_strategy=:StrongPseudoCost,
        list_of_solutions = true,
        strong_restart = true,
        debug = true,
    )

    JuMP.set_optimizer(m, optimizer_with_attributes(
        Juniper.Optimizer,
        juniper_all_solutions...
        )
    )

    @constraint(m, x[1] >= 0.9)
    @constraint(m, x[1] <= 1.1)
    @NLconstraint(m, (x[1]-x[2])^2 >= 0.1)
    @NLconstraint(m, (x[2]-x[3])^2 >= 0.1)
    @NLconstraint(m, (x[1]-x[3])^2 >= 0.1)
    @NLconstraint(m, (x[1]-x[4])^2 >= 0.1)
    @NLconstraint(m, (x[2]-x[4])^2 >= 0.1)
    @NLconstraint(m, (x[3]-x[4])^2 >= 0.1)

    JuMP.set_optimizer_attribute(m, "all_solutions", true)
    JuMP.optimize!(m)
    bm = JuMP.backend(m)
    status = MOI.get(bm, MOI.TerminationStatus())

    innermodel = bm.optimizer.model.inner
    debugDict = innermodel.debugDict
    @test getnstate(debugDict,:Integral) == 24
    @test different_hashes(debugDict) == true
    counter_test(debugDict,Juniper.getnbranches(innermodel))

    list_of_solutions = Juniper.getsolutions(innermodel)
    @test length(unique(list_of_solutions)) == Juniper.getnsolutions(innermodel)

    @test status == MOI.LOCALLY_SOLVED
    @test Juniper.getnsolutions(innermodel) == 24
end

@testset "bruteforce optimizer without attributes" begin
    m = Model(Juniper.Optimizer)
    set_optimizer_attribute(m, "nl_solver", Ipopt.Optimizer)
    set_optimizer_attribute(m, "mip_solver", Cbc.Optimizer)

    @variable(m, 1 <= x[1:4] <= 5, Int)

    @objective(m, Min, x[1])

    @constraint(m, x[1] >= 0.9)
    @constraint(m, x[1] <= 1.1)
    @NLconstraint(m, (x[1]-x[2])^2 >= 0.1)
    @NLconstraint(m, (x[2]-x[3])^2 >= 0.1)
    @NLconstraint(m, (x[1]-x[3])^2 >= 0.1)
    @NLconstraint(m, (x[1]-x[4])^2 >= 0.1)
    @NLconstraint(m, (x[2]-x[4])^2 >= 0.1)
    @NLconstraint(m, (x[3]-x[4])^2 >= 0.1)

    mktemp() do path,io
        out = stdout
        err = stderr
        redirect_stdout(io)
        redirect_stderr(io)

        status = try
            optimize!(m)
        catch e
            e
        end

        flush(io)
        redirect_stdout(out)
        redirect_stderr(err)
    end
    @test length(sort(unique(convert.(Int, JuMP.value.(x))))) == 4
    @test all(1 .<= JuMP.value.(x) .<= 5)
    @test JuMP.value(x[1]) ≈ 1.0
end

@testset "bruteforce full strong w/o restart" begin
    println("==================================")
    println("Bruteforce full strong w/o restart")
    println("==================================")
    juniper_all_solutions = DefaultTestSolver(
        branch_strategy=:StrongPseudoCost,
        all_solutions = true,
        list_of_solutions = true,
        strong_branching_perc = 100,
        strong_branching_nsteps = 100,
        strong_restart = false
    )

    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        juniper_all_solutions...)
    )

    @variable(m, 1 <= x[1:4] <= 5, Int)

    @objective(m, Min, x[1])

    @constraint(m, x[1] >= 0.9)
    @constraint(m, x[1] <= 1.1)
    @NLconstraint(m, (x[1]-x[2])^2 >= 0.1)
    @NLconstraint(m, (x[2]-x[3])^2 >= 0.1)
    @NLconstraint(m, (x[1]-x[3])^2 >= 0.1)
    @NLconstraint(m, (x[1]-x[4])^2 >= 0.1)
    @NLconstraint(m, (x[2]-x[4])^2 >= 0.1)
    @NLconstraint(m, (x[3]-x[4])^2 >= 0.1)

    status = solve(m)

    list_of_solutions = Juniper.getsolutions(internalmodel(m))
    @test length(unique(list_of_solutions)) == Juniper.getnsolutions(internalmodel(m))

    @test status == MOI.LOCALLY_SOLVED
    @test Juniper.getnsolutions(internalmodel(m)) == 24
end

@testset "bruteforce approx time limit" begin
    println("==================================")
    println("Bruteforce  approx time limit")
    println("==================================")
    juniper_all_solutions = DefaultTestSolver(
        branch_strategy=:StrongPseudoCost,
        strong_branching_approx_time_limit=2,
        all_solutions = true,
        list_of_solutions = true,
        strong_restart = true
    )

    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        juniper_all_solutions...)
    )

    @variable(m, 1 <= x[1:4] <= 5, Int)

    @objective(m, Min, x[1])

    @constraint(m, x[1] >= 0.9)
    @constraint(m, x[1] <= 1.1)
    @NLconstraint(m, (x[1]-x[2])^2 >= 0.1)
    @NLconstraint(m, (x[2]-x[3])^2 >= 0.1)
    @NLconstraint(m, (x[1]-x[3])^2 >= 0.1)
    @NLconstraint(m, (x[1]-x[4])^2 >= 0.1)
    @NLconstraint(m, (x[2]-x[4])^2 >= 0.1)
    @NLconstraint(m, (x[3]-x[4])^2 >= 0.1)

    status = solve(m)

    list_of_solutions = Juniper.getsolutions(internalmodel(m))
    @test length(unique(list_of_solutions)) == Juniper.getnsolutions(internalmodel(m))

    @test status == MOI.LOCALLY_SOLVED
    @test Juniper.getnsolutions(internalmodel(m)) == 24
end

@testset "bruteforce time limit reliable" begin
    println("==================================")
    println("Bruteforce  time reliable")
    println("==================================")
    juniper_all_solutions = DefaultTestSolver(
        branch_strategy=:Reliability,
        strong_branching_time_limit=0.02,
        reliability_branching_perc=100,
        all_solutions = true,
        list_of_solutions = true,
        strong_restart = true
    )

    optimizer = optimizer_with_attributes(Juniper.Optimizer, juniper_all_solutions...)

    m = Model(optimizer)
    moi_optimizer = JuMP.backend(m).optimizer.model

    # default => 1
    @test MOI.supports(moi_optimizer, MOI.NumberOfThreads()) == true
    MOI.set(moi_optimizer, MOI.NumberOfThreads(), nothing)

    @variable(m, 1 <= x[1:4] <= 5, Int)

    @objective(m, Min, x[1])

    @constraint(m, x[1] >= 0.9)
    @constraint(m, x[1] <= 1.1)
    @NLconstraint(m, (x[1]-x[2])^2 >= 0.1)
    @NLconstraint(m, (x[2]-x[3])^2 >= 0.1)
    @NLconstraint(m, (x[1]-x[3])^2 >= 0.1)
    @NLconstraint(m, (x[1]-x[4])^2 >= 0.1)
    @NLconstraint(m, (x[2]-x[4])^2 >= 0.1)
    @NLconstraint(m, (x[3]-x[4])^2 >= 0.1)

    status = solve(m)

    list_of_solutions = Juniper.getsolutions(internalmodel(m))
    @test length(unique(list_of_solutions)) == Juniper.getnsolutions(internalmodel(m))
    @test MOI.get(moi_optimizer, MOI.NumberOfThreads()) == 1
    # all solutions are saved => nsolutions should equal length(solutions)
    @test MOI.get(moi_optimizer, MOI.ResultCount()) == Juniper.getnsolutions(internalmodel(m))

    @test status == MOI.LOCALLY_SOLVED
    @test Juniper.getnsolutions(internalmodel(m)) == 24
end

@testset "bruteforce PseudoCost" begin
    println("==================================")
    println("Bruteforce PseudoCost")
    println("==================================")
    juniper_all_solutions = DefaultTestSolver(
        branch_strategy=:PseudoCost,
        all_solutions = true,
        list_of_solutions = true,
        strong_restart = true
    )

    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        juniper_all_solutions...)
    )

    @variable(m, 1 <= x[1:4] <= 5, Int)

    @objective(m, Min, x[1])

    @constraint(m, x[1] >= 0.9)
    @constraint(m, x[1] <= 1.1)
    @NLconstraint(m, (x[1]-x[2])^2 >= 0.1)
    @NLconstraint(m, (x[2]-x[3])^2 >= 0.1)
    @NLconstraint(m, (x[1]-x[3])^2 >= 0.1)
    @NLconstraint(m, (x[1]-x[4])^2 >= 0.1)
    @NLconstraint(m, (x[2]-x[4])^2 >= 0.1)
    @NLconstraint(m, (x[3]-x[4])^2 >= 0.1)

    status = solve(m)
    println("Status: ", status)
    list_of_solutions = Juniper.getsolutions(internalmodel(m))
    @test length(unique(list_of_solutions)) == Juniper.getnsolutions(internalmodel(m))

    @test status == MOI.LOCALLY_SOLVED
    @test Juniper.getnsolutions(internalmodel(m)) == 24
end

@testset "bruteforce Reliability" begin
    println("==================================")
    println("Bruteforce Reliability")
    println("==================================")
    juniper_all_solutions = DefaultTestSolver(
        branch_strategy=:Reliability,
        all_solutions = true,
        list_of_solutions = true,
    )

    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        juniper_all_solutions...)
    )

    @variable(m, 1 <= x[1:4] <= 5, Int)

    @objective(m, Min, x[1])

    @constraint(m, x[1] >= 0.9)
    @constraint(m, x[1] <= 1.1)
    @NLconstraint(m, (x[1]-x[2])^2 >= 0.1)
    @NLconstraint(m, (x[2]-x[3])^2 >= 0.1)
    @NLconstraint(m, (x[1]-x[3])^2 >= 0.1)
    @NLconstraint(m, (x[1]-x[4])^2 >= 0.1)
    @NLconstraint(m, (x[2]-x[4])^2 >= 0.1)
    @NLconstraint(m, (x[3]-x[4])^2 >= 0.1)

    status = solve(m)
    println("Status: ", status)
    list_of_solutions = Juniper.getsolutions(internalmodel(m))
    @test length(unique(list_of_solutions)) == Juniper.getnsolutions(internalmodel(m))

    @test status == MOI.LOCALLY_SOLVED
    @test Juniper.getnsolutions(internalmodel(m)) == 24
end

@testset "no integer" begin
    println("==================================")
    println("no integer")
    println("==================================")

    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        juniper_strong_restart...)
    )

    println("Create variables/constr/obj")
    @variable(m, 1 <= x <= 5, start = 2.7)
    @variable(m, -2 <= y <= 2)

    @objective(m, Min, -x-y)

    @NLconstraint(m, y==2*cos(2*x))
    println("before solve")
    status = solve(m)
    println("Status: ", status)

    @test status == MOI.LOCALLY_SOLVED
end

@testset "infeasible cos" begin
    println("==================================")
    println("Infeasible cos")
    println("==================================")
    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        juniper_strong_restart...)
    )

    @variable(m, 1 <= x <= 5, Int)
    @variable(m, -2 <= y <= 2, Int)

    @objective(m, Min, -x-y)

    @NLconstraint(m, y==2*cos(2*x))

    status = solve(m)

    println("Status: ", status)

    @test status == MOI.LOCALLY_INFEASIBLE
    @test JuMP.termination_status(m) == MOI.LOCALLY_INFEASIBLE
    @test JuMP.primal_status(m) == MOI.INFEASIBLE_POINT
    @test JuMP.dual_status(m) == MOI.INFEASIBLE_POINT
    @test isnan(getobjgap(m))
end

@testset "infeasible int reliable" begin
    println("==================================")
    println("Infeasible int reliable")
    println("==================================")
    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        juniper_reliable_restart...)
    )

    @variable(m, 1 <= x <= 5, Int)
    @variable(m, -2 <= y <= 2, Int)

    @objective(m, Min, -x-y)

    @NLconstraint(m, y >= sqrt(2))
    @NLconstraint(m, y <= sqrt(3))

    status = solve(m)
    println("Status: ", status)

    @test status == MOI.LOCALLY_INFEASIBLE
    @test isnan(getobjgap(m))
end

@testset "infeasible sin with different bounds" begin
    println("==================================")
    println("Infeasible  sin with different bounds")
    println("==================================")
    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver(
            branch_strategy=:MostInfeasible,
            feasibility_pump = true,
            time_limit = 1,
            mip_solver=optimizer_with_attributes(Cbc.Optimizer)
        )...)
    )

    @variable(m, x, Int, start=3)
    @variable(m, y >= 2, Int)

    @constraint(m, 0 <= x <= 5)

    @objective(m, Min, -x-y)

    @NLconstraint(m, y==sin(x))

    status = solve(m)

    @test status == MOI.LOCALLY_INFEASIBLE
end

@testset "infeasible relaxation" begin
    println("==================================")
    println("Infeasible relaxation")
    println("==================================")

    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver(;debug=true)...)
    )

    @variable(m, 0 <= x[1:10] <= 2, Int)

    @objective(m, Min, sum(x))

    @constraint(m, sum(x[1:5]) <= 20)
    @NLconstraint(m, x[1]*x[2]*x[3] >= 10)

    status = solve(m)

    debug1 = internalmodel(m).debugDict

    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver(;debug=true)...)
    )

    @variable(m, 0 <= x[1:10] <= 2, Int)

    @objective(m, Min, sum(x))

    @constraint(m, sum(x[1:5]) <= 20)
    @NLconstraint(m, x[1]*x[2]*x[3] >= 10)

    status = solve(m)

    debug2 = internalmodel(m).debugDict
    opts = internalmodel(m).options

    # should be deterministic
    @test debug1[:relaxation][:nrestarts] == debug2[:relaxation][:nrestarts] == opts.num_resolve_root_relaxation
    for i=1:debug1[:relaxation][:nrestarts]
        @test debug1[:relaxation][:restarts][i] == debug2[:relaxation][:restarts][i]
    end

    println("Status: ", status)

    @test status == MOI.LOCALLY_INFEASIBLE
end

@testset "infeasible relaxation 2" begin
    println("==================================")
    println("Infeasible relaxation 2")
    println("==================================")

    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        juniper_strong_no_restart...)
    )

    @variable(m, x[1:3], Int)
    @variable(m, y)

    @objective(m, Max, sum(x))

    @NLconstraint(m, x[1]^2+x[2]^2+x[3]^2+y^2 <= 3)
    @NLconstraint(m, x[1]^2*x[2]^2*x[3]^2*y^2 >= 10)

    status = solve(m)
    println("Status: ", status)

    @test status == MOI.LOCALLY_INFEASIBLE
end

@testset "infeasible integer" begin
    println("==================================")
    println("Infeasible integer")
    println("==================================")
    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        juniper_strong_no_restart...)
    )

    @variable(m, 0 <= x[1:10] <= 2, Int)

    @objective(m, Min, sum(x))

    @constraint(m, sum(x[1:5]) <= 20)
    @NLconstraint(m, x[1]*x[2]*x[3] >= 7)
    @NLconstraint(m, x[1]*x[2]*x[3] <= 7.5)

    status = solve(m)
    println("Status: ", status)

    @test status == MOI.LOCALLY_INFEASIBLE
end

@testset "infeasible in strong" begin
    println("==================================")
    println("Infeasible in strong")
    println("==================================")
    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        juniper_strong_no_restart...)
    )

    @variable(m, 0 <= x[1:5] <= 2, Int)

    @objective(m, Min, sum(x))

    @NLconstraint(m, x[3]^2 <= 2)
    @NLconstraint(m, x[3]^2 >= 1.2)

    status = solve(m)
    println("Status: ", status)

    @test status == MOI.LOCALLY_INFEASIBLE
end

@testset "One Integer small Reliable" begin
    println("==================================")
    println("One Integer small Reliable")
    println("==================================")
    m = Model()

    @variable(m, x >= 0, Int, start=2)
    @variable(m, y >= 0)
    @variable(m, 0 <= u <= 10, Int)
    @variable(m, w == 1)

    myf(x,y) = -3x-y
    function ∇f(g,x,y)
        g[1] = -3
        g[2] = -1
    end

    register_args = [:myf, 2, myf, ∇f]
    JuMP.register(m, register_args...)

    @NLobjective(m,Min,myf(x,y))

    JuMP.set_optimizer(m, optimizer_with_attributes(
        Juniper.Optimizer,
        juniper_reliable_restart...
        )
    )

    @constraint(m, 3x + 10 <= 20)
    @NLconstraint(m, y^2 <= u*w)

    status = solve(m)
    println("Obj: ", JuMP.objective_value(m))
    println("x: ", JuMP.value(x))
    println("y: ", JuMP.value(y))

    @test status == MOI.LOCALLY_SOLVED
    @test isapprox(JuMP.objective_value(m), -12.162277, atol=opt_atol)
    @test isapprox(JuMP.value(x), 3, atol=sol_atol)
    @test isapprox(JuMP.value(y), 3.162277, atol=sol_atol)
end

@testset "One Integer small Strong" begin
    println("==================================")
    println("One Integer small Strong")
    println("==================================")

    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        juniper_strong_no_restart...)
    )

    @variable(m, x >= 0, Int)
    @variable(m, y >= 0)
    @variable(m, 0 <= u <= 10, Int)
    @variable(m, w == 1)

    @objective(m, Min, -3x - y)

    @constraint(m, 3x + 10 <= 20)
    @NLconstraint(m, y^2 <= u*w)

    status = solve(m)
    println("Obj: ", JuMP.objective_value(m))
    println("x: ", JuMP.value(x))
    println("y: ", JuMP.value(y))

    @test status == MOI.LOCALLY_SOLVED
    @test isapprox(JuMP.objective_value(m), -12.162277, atol=opt_atol)
    @test isapprox(JuMP.value(x), 3, atol=sol_atol)
    @test isapprox(JuMP.value(y), 3.162277, atol=sol_atol)
end

@testset "One Integer small MostInfeasible" begin
    println("==================================")
    println("One Integer small MostInfeasible")
    println("==================================")

    m = Model()

    @variable(m, x >= 0, Int)
    @variable(m, y >= 0)
    @variable(m, 0 <= u <= 10, Int)
    @variable(m, w == 1)

    special_minimizer_fct(x,y) = -3x-y
    register_args = [:special_minimizer_fct, 2, special_minimizer_fct]
    JuMP.register(m, register_args...; autodiff=true)

    @NLobjective(m,Min,special_minimizer_fct(x, y))

    @constraint(m, 3x + 10 <= 20)
    @NLconstraint(m, y^2 <= u*w)

    JuMP.set_optimizer(m, optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver(
            branch_strategy=:MostInfeasible,
        )...)
    )

    status = solve(m)
    println("Obj: ", JuMP.objective_value(m))
    println("x: ", JuMP.value(x))
    println("y: ", JuMP.value(y))

    @test status == MOI.LOCALLY_SOLVED
    @test isapprox(JuMP.objective_value(m), -12.162277, atol=opt_atol)
    @test isapprox(JuMP.value(x), 3, atol=sol_atol)
    @test isapprox(JuMP.value(y), 3.162277, atol=sol_atol)
end

@testset "One Integer small PseudoCost" begin
    println("==================================")
    println("One Integer small PseudoCost")
    println("==================================")

    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        juniper_pseudo...)
    )

    @variable(m, x >= 0, Int)
    @variable(m, y >= 0)
    @variable(m, 0 <= u <= 10, Int)
    @variable(m, w == 1)

    @objective(m, Min, -3x - y)

    @constraint(m, 3x + 10 <= 20)
    @NLconstraint(m, y^2 <= u*w)

    status = solve(m)
    println("Obj: ", JuMP.objective_value(m))
    println("x: ", JuMP.value(x))
    println("y: ", JuMP.value(y))

    @test status == MOI.LOCALLY_SOLVED
    @test isapprox(JuMP.objective_value(m), -12.162277, atol=opt_atol)
    @test isapprox(JuMP.value(x), 3, atol=sol_atol)
    @test isapprox(JuMP.value(y), 3.162277, atol=sol_atol)
end

@testset "Three Integers Small Strong" begin
    println("==================================")
    println("Three Integers Small Strong")
    println("==================================")

    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        juniper_strong_no_restart...)
    )

    @variable(m, x >= 0, Int)
    @variable(m, y >= 0, Int)
    @variable(m, 0 <= u <= 10, Int)
    @variable(m, w == 1)

    @objective(m, Min, -3x - y)

    @constraint(m, 3x + 10 <= 20)
    @NLconstraint(m, y^2 <= u*w)

    status = solve(m)
    println("Obj: ", JuMP.objective_value(m))

    @test status == MOI.LOCALLY_SOLVED
    @test isapprox(JuMP.objective_value(m), -12, atol=opt_atol)
    @test isapprox(JuMP.value(x), 3, atol=sol_atol)
    @test isapprox(JuMP.value(y), 3, atol=sol_atol)
end

@testset "Three Integers Small MostInfeasible" begin
    println("==================================")
    println("Three Integers Small MostInfeasible")
    println("==================================")

    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        juniper_mosti...)
    )

    @variable(m, x >= 0, Int)
    @variable(m, y >= 0, Int)
    @variable(m, 0 <= u <= 10, Int)
    @variable(m, w == 1)

    @objective(m, Min, -3x - y)

    @constraint(m, 3x + 10 <= 20)
    @NLconstraint(m, y^2 <= u*w)

    status = solve(m)
    println("Obj: ", JuMP.objective_value(m))

    @test status == MOI.LOCALLY_SOLVED
    @test isapprox(JuMP.objective_value(m), -12, atol=opt_atol)
    @test isapprox(JuMP.value(x), 3, atol=sol_atol)
    @test isapprox(JuMP.value(y), 3, atol=sol_atol)
end

@testset "Three Integers Small PseudoCost" begin
    println("==================================")
    println("Three Integers Small PseudoCost")
    println("==================================")

    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        juniper_pseudo...)
    )

    @variable(m, x >= 0, Int)
    @variable(m, y >= 0, Int)
    @variable(m, 0 <= u <= 10, Int)
    @variable(m, w == 1)

    @objective(m, Min, -3x - y)

    @constraint(m, 3x + 10 <= 20)
    @NLconstraint(m, y^2 <= u*w)

    status = solve(m)
    println("Obj: ", JuMP.objective_value(m))

    @test status == MOI.LOCALLY_SOLVED
    @test isapprox(JuMP.objective_value(m), -12, atol=opt_atol)
    @test isapprox(JuMP.value(x), 3, atol=sol_atol)
    @test isapprox(JuMP.value(y), 3, atol=sol_atol)
end

@testset "Knapsack Max" begin
    println("==================================")
    println("KNAPSACK")
    println("==================================")

    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver(;traverse_strategy=:DBFS,
            incumbent_constr=true,mip_solver=optimizer_with_attributes(Cbc.Optimizer),
            strong_branching_time_limit=1)...)
    )

    v = [10,20,12,23,42]
    w = [12,45,12,22,21]
    @variable(m, x[1:5], Bin)

    @objective(m, Max, dot(v,x))

    @NLconstraint(m, sum(w[i]*x[i]^2 for i=1:5) <= 45)

    status = solve(m)
    println("Obj: ", JuMP.objective_value(m))

    @test status == MOI.LOCALLY_SOLVED
    @test isapprox(JuMP.objective_value(m), 65, atol=opt_atol)
    @test isapprox(JuMP.objective_bound(m), 65, atol=opt_atol)
    @test isapprox(JuMP.value.(x), [0,0,0,1,1], atol=sol_atol)
end

@testset "Knapsack Max Reliable" begin
    println("==================================")
    println("KNAPSACK Reliable no restart")
    println("==================================")

    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver(;branch_strategy=:Reliability,
              strong_restart=false,strong_branching_time_limit=1,gain_mu=0.5)...)
    )

    v = [10,20,12,23,42]
    w = [12,45,12,22,21]
    @variable(m, x[1:5], Bin)

    @objective(m, Max, dot(v,x))

    @NLconstraint(m, sum(w[i]*x[i]^2 for i=1:5) <= 45)

    status = solve(m)
    println("Obj: ", JuMP.objective_value(m))

    @test status == MOI.LOCALLY_SOLVED
    @test isapprox(JuMP.objective_value(m), 65, atol=opt_atol)
    @test isapprox(JuMP.objective_bound(m), 65, atol=opt_atol)
    @test isapprox(JuMP.value.(x), [0,0,0,1,1], atol=sol_atol)
end

@testset "Integer at root" begin
    println("==================================")
    println("INTEGER AT ROOT")
    println("==================================")
    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver()...)
    )

    @variable(m, x[1:6] <= 1, Int)
    @constraint(m, x[1:6] .== 1)
    @objective(m, Max, sum(x))
    @NLconstraint(m, x[1]*x[2]*x[3]+x[4]*x[5]*x[6] <= 100)
    status = solve(m)
    @test status == MOI.LOCALLY_SOLVED
    @test isapprox(JuMP.objective_value(m), 6, atol=opt_atol)
    # objective bound doesn't exists anymore in Ipopt
    #@test isapprox(JuMP.objective_bound(m), 0, atol=opt_atol) # Ipopt return 0
end

@testset "Knapsack Max with epsilon" begin
    println("==================================")
    println("KNAPSACK with epsilon")
    println("==================================")

    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver(;traverse_strategy=:DBFS,obj_epsilon=0.5)...)
    )

    v = [10,20,12,23,42]
    w = [12,45,12,22,21]
    @variable(m, x[1:5], Bin)

    @objective(m, Max, dot(v,x))

    @NLconstraint(m, sum(w[i]*x[i]^2 for i=1:5) <= 45)

    status = solve(m)
    println("Obj: ", JuMP.objective_value(m))

    @test status == MOI.LOCALLY_SOLVED
    @test isapprox(JuMP.objective_value(m), 65, atol=opt_atol)
    @test isapprox(JuMP.value.(x), [0,0,0,1,1], atol=sol_atol)
end

@testset "Knapsack Max with epsilon too strong" begin
    println("==================================")
    println("KNAPSACK with epsilon too strong")
    println("==================================")

    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver(;traverse_strategy=:DBFS,obj_epsilon=0.1)...)
    )

    v = [10,20,12,23,42]
    w = [12,45,12,22,21]
    # set all to 1 which is infeasible otherwise incumbent solution would be found
    @variable(m, x[1:5], Bin, start=1)

    @objective(m, Max, dot(v,x))

    @NLconstraint(m, sum(w[i]*x[i]^2 for i=1:5) <= 45)

    status = solve(m)

    @test status == MOI.LOCALLY_INFEASIBLE
end

@testset "Batch.mod Restart" begin
    println("==================================")
    println("BATCH.MOD RESTART")
    println("==================================")

    m = batch_problem()

    JuMP.set_optimizer(m, optimizer_with_attributes(
        Juniper.Optimizer,
        juniper_strong_restart...)
    )

    status = solve(m)
    @test status == MOI.LOCALLY_SOLVED

    juniper_val = JuMP.objective_value(m)

    println("Solution by Juniper")
    println("obj: ", juniper_val)

    @test isapprox(juniper_val, 285506.5082, atol=opt_atol, rtol=opt_rtol)
end

@testset "Batch.mod Restart 2 Levels" begin
    println("==================================")
    println("BATCH.MOD RESTART 2 LEVELS")
    println("==================================")

    m = batch_problem()

    JuMP.set_optimizer(m, optimizer_with_attributes(
        Juniper.Optimizer,
        juniper_strong_restart_2...)
    )

    status = solve(m)
    @test status == MOI.LOCALLY_SOLVED

    juniper_val = JuMP.objective_value(m)

    println("Solution by Juniper")
    println("obj: ", juniper_val)

    @test isapprox(juniper_val, 285506.5082, atol=opt_atol, rtol=opt_rtol)
end

@testset "cvxnonsep_nsig20r.mod restart" begin
    println("==================================")
    println("cvxnonsep_nsig20r.MOD RESTART")
    println("==================================")

    m = cvxnonsep_nsig20r_problem()

    JuMP.set_optimizer(m, optimizer_with_attributes(
        Juniper.Optimizer,
        juniper_strong_restart...)
    )

    status = solve(m)
    @test status == MOI.LOCALLY_SOLVED

    juniper_val = JuMP.objective_value(m)

    println("Solution by Juniper")
    println("obj: ", juniper_val)

    @test isapprox(juniper_val, 80.9493, atol=opt_atol, rtol=opt_rtol)
end

@testset "cvxnonsep_nsig20r.mod no restart" begin
    println("==================================")
    println("cvxnonsep_nsig20r.MOD NO RESTART")
    println("==================================")

    m = cvxnonsep_nsig20r_problem()

    JuMP.set_optimizer(m, optimizer_with_attributes(
        Juniper.Optimizer,
        juniper_strong_no_restart...)
    )

    status = solve(m)
    @test status == MOI.LOCALLY_SOLVED

    juniper_val = JuMP.objective_value(m)

    println("Solution by Juniper")
    println("obj: ", juniper_val)

    @test isapprox(juniper_val, 80.9493, atol=opt_atol, rtol=opt_rtol)
end

@testset "Knapsack solution limit and table print test" begin
    println("==================================")
    println("Knapsack solution limit and table print test")
    println("==================================")

    juniper_one_solution = DefaultTestSolver(
        log_levels=[:Table],
        branch_strategy=:MostInfeasible,
        solution_limit=1
    )

    m = Model()
    v = [10,20,12,23,42]
    w = [12,45,12,22,21]
    @variable(m, x[1:5], Bin)

    @objective(m, Max, dot(v,x))

    @NLconstraint(m, sum(w[i]*x[i]^2 for i=1:5) <= 45)

    JuMP.set_optimizer(m, optimizer_with_attributes(
        Juniper.Optimizer,
        juniper_one_solution...)
    )

    status = solve(m)
    @test status == MOI.SOLUTION_LIMIT

    # maybe 2 found at the same time
    @test Juniper.getnsolutions(internalmodel(m)) <= 2
    @test Juniper.getnsolutions(internalmodel(m)) >= 1
end

@testset "bruteforce obj_epsilon" begin
    println("==================================")
    println("Bruteforce obj_epsilon")
    println("==================================")
    juniper_obj_eps = DefaultTestSolver(
        branch_strategy=:StrongPseudoCost,
        obj_epsilon=0.4
    )

    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        juniper_obj_eps...)
    )

    @variable(m, 1 <= x[1:4] <= 5, Int)

    @objective(m, Min, x[1])

    @constraint(m, x[1] >= 0.9)
    @constraint(m, x[1] <= 1.1)
    @NLconstraint(m, (x[1]-x[2])^2 >= 0.1)
    @NLconstraint(m, (x[2]-x[3])^2 >= 0.1)
    @NLconstraint(m, (x[1]-x[3])^2 >= 0.1)
    @NLconstraint(m, (x[1]-x[4])^2 >= 0.1)
    @NLconstraint(m, (x[2]-x[4])^2 >= 0.1)
    @NLconstraint(m, (x[3]-x[4])^2 >= 0.1)

    status = solve(m)

    # maybe 2 found at the same time
    @test status == MOI.LOCALLY_SOLVED
end

@testset "bruteforce best_obj_stop not reachable" begin
    println("==================================")
    println("Bruteforce best_obj_stop not reachable")
    println("==================================")
    juniper_best_obj_stop = DefaultTestSolver(
        branch_strategy=:StrongPseudoCost,
        best_obj_stop=0.8
    )
    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        juniper_best_obj_stop...)
    )

    @variable(m, 1 <= x[1:4] <= 5, Int)

    @objective(m, Min, x[1])

    @constraint(m, x[1] >= 0.9)
    @constraint(m, x[1] <= 1.1)
    @NLconstraint(m, (x[1]-x[2])^2 >= 0.1)
    @NLconstraint(m, (x[2]-x[3])^2 >= 0.1)
    @NLconstraint(m, (x[1]-x[3])^2 >= 0.1)
    @NLconstraint(m, (x[1]-x[4])^2 >= 0.1)
    @NLconstraint(m, (x[2]-x[4])^2 >= 0.1)
    @NLconstraint(m, (x[3]-x[4])^2 >= 0.1)

    status = solve(m)

    # not possible to reach but should be solved anyway
    @test status == MOI.LOCALLY_SOLVED
end

@testset "bruteforce best_obj_stop reachable" begin
    println("==================================")
    println("Bruteforce best_obj_stop reachable")
    println("==================================")
    juniper_one_solution = DefaultTestSolver(
        branch_strategy=:StrongPseudoCost,
        best_obj_stop=1
    )

    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        juniper_one_solution...)
    )

    @variable(m, 1 <= x[1:4] <= 5, Int)

    @objective(m, Min, x[1])

    @constraint(m, x[1] >= 0.9)
    @constraint(m, x[1] <= 1.1)
    @NLconstraint(m, (x[1]-x[2])^2 >= 0.1)
    @NLconstraint(m, (x[2]-x[3])^2 >= 0.1)
    @NLconstraint(m, (x[1]-x[3])^2 >= 0.1)
    @NLconstraint(m, (x[1]-x[4])^2 >= 0.1)
    @NLconstraint(m, (x[2]-x[4])^2 >= 0.1)
    @NLconstraint(m, (x[3]-x[4])^2 >= 0.1)

    status = solve(m)

    # reachable and should break
    @test status == MOI.OBJECTIVE_LIMIT
    # maybe 2 found at the same time
    @test Juniper.getnsolutions(internalmodel(m)) <= 2
    @test Juniper.getnsolutions(internalmodel(m)) >= 1
end

# this test has a lot "Only almost locally solved" warnings (in mumps at least)
@testset "Sum 1/x = 2" begin
    println("==================================")
    println("Sum 1/x = 2")
    println("==================================")
    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        juniper_pseudo...)
    )

    @variable(m, 1 <= x[1:11], Int)
    @constraint(m, [i=2:11], x[i-1] <= x[i] - 1)
    @NLconstraint(m, sum(1 / x[i] for i in 1:11) == 2)

    status = solve(m)

    # in mumps/on travis this returns ALMOST_LOCALLY_SOLVED but with ma27/locally it is LOCALLY_SOLVED
    @test Juniper.state_is_optimal(status; allow_almost=true)
    @test isapprox(2, sum(1/v for v in JuMP.value.(x)), atol=opt_atol)
end

# this test has a lot "Only almost locally solved" warnings
@testset "Sum 1/x = 2 don't allow almost" begin
    println("==================================")
    println("Sum 1/x = 2 don't allow almost")
    println("==================================")
    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver(
                branch_strategy=:PseudoCost,
                allow_almost_solved = false
            )...)
    )

    @variable(m, 1 <= x[1:11], Int)
    @constraint(m, [i=2:11], x[i-1] <= x[i] - 1)
    @NLconstraint(m, sum(1 / x[i] for i in 1:11) == 2)

    status = solve(m)

    # even without almost this should still be solveable
    @test status == MOI.LOCALLY_SOLVED
    @test isapprox(2, sum(1/v for v in JuMP.value.(x)), atol=opt_atol)
end

#this test has an expression where a variable will be dereferenced twice
@testset "Nested variable reference" begin
    println("==================================")
    println("Nested variable reference")
    println("==================================")

    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver()...)
    )

    x = @variable(m, x[i=1:5,j=1:5], Bin)
    xrowsum = @NLexpression(m, xrowsum[i=1:5], sum(x[i,j] for j in 1:5))
    @NLobjective(m, Max, sum(sum(x[i,j] + xrowsum[i] for i=1:5) for j = 1:5))

    optimize!(m)

    @test JuMP.termination_status(m) == MOI.LOCALLY_SOLVED
end

@testset "Expr dereferencing for @NLexpression" begin
    # See issue 184
    m = Model()
    set_optimizer(m, optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver(
            mip_solver=optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0))...
    ))

    @variable(m, 0 <= a_var <= 1)
    @variable(m, bin_var, Int)
    b_expr = @NLexpression(m, bin_var/1)
    @NLconstraint(m, 1.1 >= b_expr)
    an_expr = @NLexpression(m, a_var / 1)

    @NLobjective(m, Max, bin_var + an_expr)

    optimize!(m)
    @test JuMP.objective_value(m) ≈ 2.0
    @test JuMP.value(bin_var) ≈ 1
    @test JuMP.value(a_var) ≈ 1
end

end


text:
using AntennaExpts, YagiAnt_GP

yagiant_gp()

ext: jl


text:
module BifurcationKit
	using Printf, Dates, LinearMaps, BlockArrays, RecipesBase, StructArrays, Requires
	using Setfield: setproperties, @set, @set!, Lens, get, set, @lens
	using Parameters: @with_kw, @unpack, @with_kw_noshow
	using RecursiveArrayTools: VectorOfArray
	using DocStringExtensions
	using DataStructures: CircularBuffer
	using ForwardDiff

	# we put this here to be used in LinearBorderSolver
	abstract type AbstractContinuationIterable end
	abstract type AbstractContinuationState end

	include("BorderedArrays.jl")
	include("LinearSolver.jl")
	include("EigSolver.jl")
	include("LinearBorderSolver.jl")
	include("Preconditioner.jl")
	include("Newton.jl")
	include("ContParameters.jl")
	include("Results.jl")

	include("events/Event.jl")

	include("Continuation.jl")
	include("events/EventDetection.jl")
	include("events/BifurcationDetection.jl")

	include("Bifurcations.jl")
	include("Predictor.jl")

	include("DeflationOperator.jl")
	include("BorderedProblem.jl")

	include("Utils.jl")

	include("codim2/codim2.jl")
	include("codim2/MinAugFold.jl")
	include("codim2/MinAugHopf.jl")

	include("BifurcationPoints.jl")

	include("bifdiagram/BranchSwitching.jl")
	include("NormalForms.jl")
	include("codim2/BifurcationPoints.jl")
	include("codim2/NormalForms.jl")
	include("bifdiagram/BifurcationDiagram.jl")

	include("DeflatedContinuation.jl")

	include("periodicorbit/Sections.jl")
	include("periodicorbit/PeriodicOrbits.jl")
	include("periodicorbit/PeriodicOrbitTrapeze.jl")
	include("periodicorbit/PeriodicOrbitCollocation.jl")
	# include("periodicorbit/PeriodicOrbitMIRK.jl")
	# include("periodicorbit/PeriodicOrbitFDAdapt.jl")
	include("periodicorbit/PeriodicOrbitUtils.jl")
	include("periodicorbit/Flow.jl")
	include("periodicorbit/FlowDE.jl")
	include("periodicorbit/StandardShooting.jl")
	include("periodicorbit/PoincareShooting.jl")
	include("periodicorbit/ShootingDE.jl")
	include("periodicorbit/Floquet.jl")

	include("wave/WaveProblem.jl")
	include("wave/EigSolver.jl")

	include("plotting/Recipes.jl")
	include("Diffeqwrap.jl")

	using Requires

	function __init__()
		# if Plots.jl is available, then we allow plotting of solutions
		@require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
			using .Plots
			include("plotting/PlotCont.jl")
		end
		@require AbstractPlotting="537997a7-5e4e-5d89-9595-2241ea00577e" begin
			using .AbstractPlotting: @recipe, inline!, layoutscene, Figure, Axis, lines!
			include("plotting/RecipesMakie.jl")
		end

		@require GLMakie="e9467ef8-e4e7-5192-8a1a-b1aee30e663a" begin
			@info "Loading GLMakie code"
			using .GLMakie: @recipe, inline!, layoutscene, Figure, Axis, lines!, PointBased, Point2f0, scatter!
			include("plotting/RecipesMakie.jl")
		end

		@require JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" begin
			using .JLD2
			"""
			Save solution / data in JLD2 file
			- `filename` is for example "example.jld2"
			- `sol` is the solution
			- `p` is the parameter
			- `i` is the index of the solution to be saved
			"""
			function saveToFile(iter::AbstractContinuationIterable, sol, p, i::Int64, br::ContResult)
				if iter.contParams.saveToFile == false; return nothing; end
				filename = iter.filename
				# this allows to save two branches forward/backward in case
				# bothside = true is passed to continuation
				fd = iter.contParams.ds >=0 ? "fw" : "bw"

					# create a group in the JLD format
					jldopen(filename*".jld2", "a+") do file
					if haskey(file, "sol-$fd-$i")
						delete!(file, "sol-$fd-$i")
					end
					mygroup = JLD2.Group(file, "sol-$fd-$i")
						mygroup["sol"] = sol
						mygroup["param"] = p
					end

				jldopen(filename*"-branch.jld2", "a+") do file
					if haskey(file, "branch"*fd)
						delete!(file, "branch"*fd)
					end
					file["branch"*fd] = br
					end
			end

			# final save of branch, in case bothsided = true is used
			function saveToFile(iter::AbstractContinuationIterable, br::ContResult)
				if iter.contParams.saveToFile == false; return nothing; end
				filename = iter.filename

				jldopen(filename*"-branch.jld2", "a+") do file
					if haskey(file, "branchfw")
						delete!(file, "branchfw")
					end
					if haskey(file, "branchbw")
						delete!(file, "branchbw")
					end
					if haskey(file, "branch")
						delete!(file, "branch")
					end
					file["branch"] = br
				end
			end
		end

	end

	# linear solvers
	export DefaultLS, GMRESIterativeSolvers, GMRESKrylovKit,
			DefaultEig, EigArpack, EigIterativeSolvers, EigKrylovKit, EigArnoldiMethod, geteigenvector, AbstractEigenSolver

	# bordered nonlinear problems
	export BorderedProblem, JacobianBorderedProblem, LinearSolverBorderedProblem

	# preconditioner based on deflation
	export PrecPartialSchurKrylovKit, PrecPartialSchurArnoldiMethod

	# bordered linear problems
	export MatrixBLS, BorderingBLS, MatrixFreeBLS, LSFromBLS, BorderedArray

	# nonlinear deflation
	export DeflationOperator, DeflatedProblem, DeflatedLinearSolver, scalardM

	# predictors for continuation
	export SecantPred, BorderedPred, NaturalPred, MultiplePred, PolynomialPred

	# newton methods
	export NewtonPar, newton, newtonDeflated, newtonPALC, newtonFold, newtonHopf, newtonBordered

	# continuation methods
	export ContinuationPar, ContResult, GenericBifPoint, continuation, continuation!, continuationFold, continuationHopf, continuationPOTrap, continuationBordered, eigenvec, eigenvals, getSolx, getSolp

	# events
	export ContinuousEvent, DiscreteEvent, PairOfEvents, SetOfEvents, SaveAtEvent, FoldDetectEvent, BifDetectEvent

	# iterators for continuation
	export ContIterable, iterate, ContState, solution, getx, getp

	# codim2 Fold continuation
	export FoldPoint, FoldProblemMinimallyAugmented, FoldLinearSolveMinAug, foldPoint

	# codim2 Hopf continuation
	export HopfPoint, HopfProblemMinimallyAugmented, HopfLinearSolveMinAug

	# normal form
	export computeNormalForm, predictor

	# automatic bifurcation diagram
	export bifurcationdiagram, bifurcationdiagram!, Branch, BifDiagNode, getBranch, getBranchesFromBP

	# Periodic orbit computation
	export generateSolution, getPeriod, getAmplitude, getMaximum, getPeriodicOrbit, sectionSS, sectionPS, guessFromHopf

	# Periodic orbit computation based on Trapeze method
	export PeriodicOrbitTrapProblem, continuationPOTrap, continuationPOTrapBPFromPO

	# Periodic orbit computation based on Shooting
	export Flow, ShootingProblem, PoincareShootingProblem, continuationPOShooting, AbstractShootingProblem, SectionPS, SectionSS

	# Periodic orbit computation based on Collocation
	export PeriodicOrbitOCollProblem

	# Floquet multipliers computation
	export FloquetQaD

	# waves
	export TWProblem
end


text:
#!/usr/bin/julia

function find_sum(num_list::Array{T}, try_sum::T) where {T<:Number}
    sorted_list = sort(num_list)
    unique!(sorted_list)
    skip = false
    found_nums = zeros(T, 2)
    found_prod = zero(T)
    counter = 0
    for i in sorted_list
        for j in Iterators.reverse(sorted_list)
            counter += 1
            if i >= j
                skip = true
                break
            elseif i + j < try_sum
                break
            elseif i + j == try_sum
                found_nums[:] = [i, j]
                found_prod = i * j
                skip = true
                break
            end
        end
        skip && break
    end
    return found_nums, found_prod, counter
end

function find_triplet_sum(num_list::Array{T}, try_sum::T) where {T<:Number}
    sorted_list = sort(num_list)
    unique!(sorted_list)
    skip = false
    found_nums = zeros(T, 3)
    found_prod = zero(T)
    counter = 0
    for i in sorted_list
        skip_inner = false
        for j in sorted_list[2:end]
            outer_sum = i + j
            for k in Iterators.reverse(sorted_list)
                counter += 1
                inner_sum = outer_sum + k
                if inner_sum < try_sum
                    break
                elseif j >= k
                    skip_inner = true
                    break
                elseif inner_sum == try_sum
                    found_nums[:] = [i, j, k]
                    skip = true
                    break
                end
            end
            (skip_inner || skip) && break
        end
        skip && break
    end
    found_prod = prod(found_nums)
    return found_nums, found_prod, counter
end

function main()
    num_list = open("input/01.txt", "r") do io
        map(s->tryparse(Int, s), readlines(io))
    end

    @time find_sum(num_list, 2020)
    @time find_sum(num_list, 2020)
    @info "Found answer" find_sum(num_list, 2020)

    @time find_triplet_sum(num_list, 2020)
    @time find_triplet_sum(num_list, 2020)
    @info "Found answer" find_triplet_sum(num_list, 2020)
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end


text:
"""
    FeedForward(NeuralNetwork)

Feed forward process for the neural network.

    Usage:
```julia-repl
FeedForward(NeuralNetwork)
```
Arguments:
- `NeuralNetwork` : The NeuralNetwork structure.
"""
function FeedForward(NeuralNetwork)
    NeuralNetwork.layer1 = sigmoid(NeuralNetwork.input * NeuralNetwork.w₁)
    NeuralNetwork.output = sigmoid(NeuralNetwork.layer1 * NeuralNetwork.w₂)
end


text:
function check_domterm()
    global isDomterm

    isDomterm = true
    try
        if Sys.isunix()
            run(`domterm is-domterm`)
        else
            println("Domterm can only run in Unix or WSL")
        end
    catch err
        println("Not running in domterm\n")
        isDomterm = false
    end
end

"""add a CSS style rule to domterm

    Example ::
    ```Julia

       dt_addstyle("table {border: 1px solid}")
       dt_addstyle("table td{border-color: red}")

    ```
    The above functions set all tables in the current DomTerm session
    to have 1px-width border and red for the inner-border color.

"""
function dt_addst(rule1::String, rules...)
    args = rule1*" "*join(rules, ' ')
    r = String(read(`domterm add-style $args`))
    return nothing
end

function dt_settings(args...)
    arg = args.join(' ')
    r = run(`domterm $arg`)
    return nothing
end

function dt_reversevideo(turn_on::Bool=true)
    arg = turn_on ? "on" : "off"
    run(`domterm reverse-video $arg`)
    return nothing
end

function dt_status()
    s = read(`domterm status`, String)
    return s
end

function dt_list()
    s = read(`domterm list`, String)
    return s
end

function dt_listst()
    s = read(`domterm list-stylesheets`, String)
    return s
end

function dt_printst(i)
    s = read(`domterm print-stylesheet $i`, String)
    return s
end

using Base64: Base64EncodePipe

const img_style="display:block; margin:auto;"

function imagecat(mime::AbstractString, fn::AbstractString; attr="")
    at = img_style
    if !isempty(attr); at *= attr; end
    x = open(fn) do io
        read(io)
    end # x is Vector{UInt8}
    buf = IOBuffer()
    b64 = Base64EncodePipe(buf)
    write(b64, x)
    close(b64)
    uri = String(take!(buf))
    print("\e]72;<img style=\" $(at) \" src='data:$(mime);base64,", uri, "'/>\a")
end

svgcat(fn::AbstractString; attr="") = imagecat("image/svg+xml", fn; attr)
pngcat(fn::AbstractString; attr="") = imagecat("image/png", fn; attr)
jpegcat(fn::AbstractString; attr="") = imagecat("image/jpeg", fn; attr)
gifcat(fn::AbstractString; attr="") = imagecat("image/gif", fn; attr)

raw_htmlcat(html_content::String) = print("\e]72;"*html_content*"\a")
hput(html_content::String) = raw_htmlcat(html_content)

function htmlcat(fn::AbstractString)
    x = String(open(io -> read(io), fn))

    # remove elements not recognized by domterm
    x = replace(x, r"<\?xml.*?>"s=>"")
    x = replace(x, r"<!DOCTYPE.*?>"s=>"")
    x = replace(x, r"<!doctype.*?>"s=>"")
    x = replace(x, r"<mfenced.*?>"s=>"")
    x = replace(x, r"</mfenced>"s=>"")
    #x = replace(x, r"\n"=>"") # strip extra spaces off the output

    raw_htmlcat(x)
end

dt_clear() = println("\e[7J")

# Copied from DomTerm documentation
# "\e]721;" key ";" html-text "\a"
#
# Replace previously-inserted HTML. Looks for the latest element (in document
# order) whose class includes can-replace-children and that has a
# replace-key="key" attribute. Essentially sets the innerHTML to html-text
# (after safety-scrubbing).

const dt_cmd = Dict(
    :reversevid => dt_reversevideo,
    :clear => dt_clear,
    :status => dt_status,
    :list => dt_list,
    :listst => dt_listst,
    :printst => dt_printst,
    :addst => dt_addst,
    :set => dt_settings,
    :jpegcat => jpegcat,
    :pngcat => pngcat,
    :svgcat => svgcat,
    :gifcat => gifcat,
    :htmlcat => htmlcat
)
macro dt(cmd, args...)
    if isempty(args)
        expr = Expr(:call, dt_cmd[cmd])
    else
        expr = Expr(:call, dt_cmd[cmd], args...)
    end

    r = eval(expr)
    r !== nothing && println(r)
    return nothing
end

# get the size in pixel of the current terminal
iolock_begin() = ccall(:jl_iolock_begin, Cvoid, ())
iolock_end() = ccall(:jl_iolock_end, Cvoid, ())

function tty_set_raw!(io::IO, isRaw::Bool)
    iolock_begin()
    ccall(:jl_tty_set_mode, Int32, (Ptr{Cvoid},Int32), io.handle::Ptr{Cvoid}, isRaw)
    iolock_end()
end

function dt_winsizex()

    h = w = -1

    tty_set_raw!(stdin, true)
    write(stdout, "\e[14t")
    out = Vector{UInt8}()
    while true
        ch = read(stdin, Char)
        push!(out, ch)
        ch == 't' && break
    end
    tty_set_raw!(stdin, false)

    vh = Vector{UInt8}()
    vw = Vector{UInt8}()

    i = findfirst([UInt8(';')], out).start
    i == nothing && error("malformed xterm report")

    i += 1
    while out[i] != UInt8(';')
        push!(vw, out[i])
        i += 1
    end
    i += 1
    while out[i] != UInt8('t')
        push!(vh, out[i])
        i += 1
    end

    w = parse(Int, String(vw), base=10)
    h = parse(Int, String(vh), base=10)
    return h, w
end


text:
using LinearAlgebra
using LinearMaps
using MAT
using SparseArrays
using Arpack

using PyCall, SparseArrays

function scipyCSC_to_julia(A)
    m, n = A.shape
    colPtr = Int[i+1 for i in PyArray(A."indptr")]
    rowVal = Int[i+1 for i in PyArray(A."indices")]
    nzVal = Vector{Float64}(PyArray(A."data"))
    B = SparseMatrixCSC{Float64,Int}(m, n, colPtr, rowVal, nzVal)
    return PyCall.pyjlwrap_new(B)
end

function read_arxiv(file::String)
    I = Int64[]
    J = Int64[]
    open(file) do f
        for line in eachline(f)
            if line[1] == '#'; continue; end
            data = split(line, ",")
            push!(I, parse(Int64, data[1]))
            push!(J, parse(Int64, data[2]))
        end
    end
    I .+= 1
    J .+= 1
    n = max(maximum(I), maximum(J))
    A = sparse(I, J, 1, n, n)
    A = max.(A, A')
    A = min.(A, 1)
    return A
end

function main(PyA, k::Int64)
    m, n = PyA.shape
    colPtr = Int[i+1 for i in PyArray(PyA."indptr")]
    rowVal = Int[i+1 for i in PyArray(PyA."indices")]
    nzVal = Vector{Float64}(PyArray(PyA."data"))
    A = SparseMatrixCSC{Float64,Int}(m, n, colPtr, rowVal, nzVal)
    d = vec(sum(A, dims=2))
    τ = sum(d) / length(d)
    N = size(A)[1]

    # normalized regularized laplacian
    D = Diagonal(1.0 ./ sqrt.(d .+ τ))
    Aop = LinearMap{Float64}(X -> A * X .+ (τ / N) * sum(X), N, N, isposdef=true, issymmetric=true)
    NRL = I + D * Aop * D

    (Λ, V) = eigs(NRL, nev=k, tol=1e-6, ncv=2*k+1, which=:LM)

    # axis rotation (not necessary, but could be helpful)
    piv = qr(V', Val(true)).jpvt[1:k]
    piv_svd = svd(V[piv,:]', full=false)
    SCDM_V = V * (piv_svd.U * piv_svd.Vt)

    # save

    return SCDM_V
end

#A = read_arxiv(ARGS[1])
#embed = main(A, 128)
#matwrite("$(ARGS[2])_spectral_embedding.mat", Dict("V" => embed), compress=true)

max_stars_repo_name: amazon-research/gnn-tail-generalization


text:
abstract type AbstractDiffusion end
abstract type AbstractStaticDiffusion <: AbstractDiffusion end
abstract type AbstractDynamicDiffusion <: AbstractDiffusion end
isstatic(diffusion::AbstractStaticDiffusion) = true
isdynamic(diffusion::AbstractStaticDiffusion) = false
isstatic(diffusion::AbstractDynamicDiffusion) = false
isdynamic(diffusion::AbstractDynamicDiffusion) = true

estimate_global_diffusion(diffusion::AbstractDynamicDiffusion, d, q, Eltype) = NaN

struct DynamicDiffusion <: AbstractDynamicDiffusion end
initial_diffusion(diffusion::DynamicDiffusion, d, q, Eltype) = one(Eltype)
estimate_local_diffusion(kind::DynamicDiffusion, integ) = local_scalar_diffusion(integ)

struct DynamicMVDiffusion <: AbstractDynamicDiffusion end
initial_diffusion(diffusionmodel::DynamicMVDiffusion, d, q, Eltype) =
    kron(Diagonal(ones(Eltype, d)), Diagonal(ones(Eltype, q + 1)))
estimate_local_diffusion(kind::DynamicMVDiffusion, integ) = local_diagonal_diffusion(integ)

Base.@kwdef struct FixedDiffusion{T<:Number} <: AbstractStaticDiffusion
    initial_diffusion::T = 1.0
    calibrate::Bool = true
end
initial_diffusion(diffusionmodel::FixedDiffusion, d, q, Eltype) =
    diffusionmodel.initial_diffusion * one(Eltype)
estimate_local_diffusion(kind::FixedDiffusion, integ) = local_scalar_diffusion(integ)
function estimate_global_diffusion(rule::FixedDiffusion, integ)
    @unpack d, measurement, m_tmp = integ.cache
    # sol_diffusions = integ.sol.diffusions

    v, S = measurement.μ, measurement.Σ
    e, _ = m_tmp.μ, m_tmp.Σ.mat
    _S = copy!(m_tmp.Σ.mat, S.mat)
    if _S isa Diagonal
        e .= v ./ _S.diag
    else
        S_chol = cholesky!(_S)
        ldiv!(e, S_chol, v)
    end
    diffusion_t = dot(v, e) / d

    if integ.success_iter == 0
        # @assert length(sol_diffusions) == 0
        global_diffusion = diffusion_t
        return global_diffusion
    else
        # @assert length(sol_diffusions) == integ.success_iter
        diffusion_prev = integ.cache.global_diffusion
        global_diffusion =
            diffusion_prev + (diffusion_t - diffusion_prev) / integ.success_iter
        # @info "compute diffusion" diffusion_prev global_diffusion
        return global_diffusion
    end
end

Base.@kwdef struct FixedMVDiffusion{T} <: AbstractStaticDiffusion
    initial_diffusion::T = 1.0
    calibrate::Bool = true
end
function initial_diffusion(diffusionmodel::FixedMVDiffusion, d, q, Eltype)
    initdiff = diffusionmodel.initial_diffusion
    @assert initdiff isa Number || length(initdiff) == d
    return kron(Diagonal(initdiff .* ones(Eltype, d)), Diagonal(ones(Eltype, q + 1)))
end
estimate_local_diffusion(kind::FixedMVDiffusion, integ) = local_diagonal_diffusion(integ)
function estimate_global_diffusion(kind::FixedMVDiffusion, integ)
    @unpack q, measurement = integ.cache

    v, S = measurement.μ, measurement.Σ
    S_11 = diag(S)[1]

    Σ_ii = v .^ 2 ./ S_11
    Σ = Diagonal(Σ_ii)
    Σ_out = kron(Σ, I(q + 1))

    if integ.success_iter == 0
        # @assert length(diffusions) == 0
        return Σ_out
    else
        # @assert length(diffusions) == integ.success_iter
        diffusion_prev = integ.cache.global_diffusion
        diffusion = diffusion_prev + (Σ_out - diffusion_prev) / integ.success_iter
        return diffusion
    end
end

"""
    local_scalar_diffusion(integ)

Computes a scalar local diffusion estimate:
```math
σ² = zᵀ ⋅ (H*Q*H')⁻¹ ⋅ z
```
"""
function local_scalar_diffusion(integ)
    @unpack d, R, H, Qh, measurement, m_tmp = integ.cache
    z = measurement.μ
    e, HQH = m_tmp.μ, m_tmp.Σ
    X_A_Xt!(HQH, Qh, H)
    if HQH.mat isa Diagonal
        e .= z ./ HQH.mat.diag
        σ² = dot(e, z) / d
        return σ²
    else
        C = cholesky!(HQH.mat)
        ldiv!(e, C, z)
        σ² = dot(z, e) / d
        return σ²
    end
end

"""
    local_diagonal_diffusion(integ)

Computes a diagonal local diffusion estimate:
```math
Σᵢᵢ = zᵢ² / (H*Q*H')ᵢᵢ
```
"""
function local_diagonal_diffusion(integ)
    @unpack q, H, Qh, measurement, m_tmp = integ.cache
    z = measurement.μ
    HQH = X_A_Xt!(m_tmp.Σ, Qh, H)
    Q0_11 = diag(HQH)[1]

    Σ_ii = z .^ 2 ./ Q0_11
    # Σ_ii .= max.(Σ_ii, eps(eltype(Σ_ii)))
    Σ = Diagonal(Σ_ii)

    Σ_out = kron(Σ, I(q + 1))
    return Σ_out
end


text:
const unit_excludes = String[
    "number_threads",                 # EAGO won't support number of threads in near future
    "raw_status_string",              # TODO: ADD internal status states to EAGO
    "variablenames",                  # SEEMS LIKE A FALSE PASS

    "solve_qcp_edge_cases",           # Need to update nonconvex...
    "solve_qp_zero_offdiag",
     "solve_qp_edge_cases",

    "solve_zero_one_with_bounds_3",   # GLPK has a non-standard return code
    "solve_result_index",             # TODO: Should throw error when querying for multiple results... (expected behavior?)
]

const contlinear_excludes = String[
    "partial_start",                 # EAGO doesn't support VariablePrimalStart
]

const intlinear_excludes = String[
    "indicator1",    # doesn't currently support indicator sets
    "indicator2",    # can't check using Cbc until https://github.com/jump-dev/Cbc.jl/issues/151 is resolved
    "indicator3",
    "indicator4",
    "int2",
]

const contconic_excludes = String[
    "dualexp",  # Not directly bridged to common cones
    "dualpow",
    "logdet",
    "rootdet",
]

const contquadratic_excludes = String[
    "ncqcp",
    #"qp",
    "socp",
    "qcp",
]

const nlp_excludes = String[
    "feasibility_sense_with_objective_and_hessian",
    "feasibility_sense_with_no_objective_and_no_hessian",
    "feasibility_sense_with_objective_and_no_hessian",
    "hs071_no_hessian",
    "feasibility_sense_with_no_objective_and_with_hessian",
    "nlp_objective_and_moi_objective",
    "hs071",
    "hs071_hessian_vector_product_test"
]

function test_moi(T::Type{<:Real}; solver_options...)

    optimizer = MOIU.CachingOptimizer(MOIU.UniversalFallback(MOIU.Model{T}()), EAGO.Optimizer(; solver_options...))

    tol = 2sqrt(sqrt(eps(T)))
    config = MOIT.TestConfig{T}(
    atol = tol,
    rtol = tol,
    solve = true,
    query = true,
    modify_lhs = false,
    duals = false,
    infeas_certificates = false,
    )

    @testset "unit tests" begin
        MOIT.unittest(MOIB.full_bridge_optimizer(optimizer, T), config, unit_excludes)
    end

    @testset "continuous linear tests" begin
        MOIT.contlineartest(MOIB.full_bridge_optimizer(optimizer, T), config, contlinear_excludes)
    end

    @testset "mixed-integer linear tests" begin
        MOIT.intlineartest(MOIB.full_bridge_optimizer(optimizer, T), config, intlinear_excludes)
    end

    @testset "continuous conic tests" begin
        MOIT.contconictest(MOIB.full_bridge_optimizer(optimizer, T), config, contconic_excludes)
    end

    @testset "continuous quadratic tests" begin
        MOIT.contquadratictest(MOIB.full_bridge_optimizer(optimizer, T), config, contquadratic_excludes)
    end

    @testset "nlp" begin
        MOIT.nlptest(MOIB.full_bridge_optimizer(optimizer, T), config, nlp_excludes)
    end
end

# Test with mip_solver = Cbc as it supports SOS1 & SOS2 constraints
# TODO: Use bridges for SOS1 & SOS2 constraint if unsupported
# Need to test with GLPK as well to ensure subsolver supports constraint
# coefficient modification.
default_options = (:verbosity => 0,)
test_moi(Float64; default_options...)


text:
aswconsts = Constants{Float64}
aswrefs = ReferenceValues{Float64}
aswweight = Weight{Float64}
aswsensor = Sensor{Float64}
aswengine = Engine{Float64}
aswstrut = Strut{Float64}
aswjoint = Joint{Float64}
aswjangle = Jangle{Float64}
aswground = Ground{Float64}
aswbeamspanwise = SpanwiseVariables{Float64}
aswbeam = Beam{Float64}
aswgeominput = Configuration{Float64}

aswconstraints = Constraints
aswparameters = Parameters
operatingpoint = OperatingPoint

stabderivs_s = StabilityDerivatives
solutionparameters = SolutionParameters
solutiongeometry = SolutionGeometry
solutionspanwise = SolutionSpanwise
solution = Solution


text:
############################################################
# judiPDEfull ################################################
############################################################

# Authors: Philipp Witte (pwitte@eos.ubc.ca), Henryk Modzelewski (hmodzelewski@eos.ubc.ca)
# Date: January 2017

export judiPDEfull, judiPDEfullException, subsample

############################################################

# Type for linear operator representing  Pr*A(m)^-1*Ps,
# i.e. it includes source and receiver projections

struct judiPDEfull{DDT<:Number,RDT<:Number} <: joAbstractLinearOperator{DDT,RDT}
    name::String
    m::Integer
    n::Integer
    info::Info
    model::Model
    srcGeometry::Geometry
    recGeometry::Geometry
    options::Options
    fop::Function  # forward
    fop_T::Union{Function, Nothing}  # transpose
end

mutable struct judiPDEfullException <: Exception
    msg :: String
end

############################################################
## Constructor

# Set up info structure for linear operators
function judiModeling(model::Model, srcGeometry::Geometry, recGeometry::Geometry; options=Options(), DDT::DataType=Float32, RDT::DataType=DDT)
    ntComp = get_computational_nt(srcGeometry, recGeometry, model)
    info = Info(prod(model.n), get_nsrc(srcGeometry), ntComp)
    return judiModeling(info, model, srcGeometry, recGeometry; options=options, DDT=DDT, RDT=RDT)
end

function judiModeling(info::Info,model::Model, srcGeometry::Geometry, recGeometry::Geometry; options=Options(), DDT::DataType=Float32, RDT::DataType=DDT)
# JOLI wrapper for nonlinear forward modeling
    (DDT == Float32 && RDT == Float32) || throw(judiPDEfullException("Domain and range types not supported"))

    # Determine dimensions
    m, n = n_samples(recGeometry, info), n_samples(srcGeometry, info)

    srcnum = 1:info.nsrc

    return F = judiPDEfull{Float32,Float32}("Proj*F*Proj'", m, n, info, model, srcGeometry, recGeometry, options,
                              src -> time_modeling(model, srcGeometry, process_input_data(src, srcGeometry, info),
                                                   recGeometry, nothing, nothing, srcnum, 'F', 1, options),
                              rec -> time_modeling(model, srcGeometry, nothing, recGeometry,
                                                   process_input_data(rec, recGeometry, info), nothing, srcnum, 'F', -1, options),
                              )
end

############################################################
## overloaded Base functions

# conj(judiPDEfull)
conj(A::judiPDEfull{DDT,RDT}) where {DDT,RDT} =
    judiPDEfull{DDT,RDT}("conj("*A.name*")",A.m,A.n,A.info,A.model,A.srcGeometry,A.recGeometry,A.options,
        A.fop,
        A.fop_T
        )

# transpose(judiPDEfull)
transpose(A::judiPDEfull{DDT,RDT}) where {DDT,RDT} =
    judiPDEfull{DDT,RDT}("Proj*F'*Proj'",A.n,A.m,A.info,A.model,A.srcGeometry,A.recGeometry,A.options,
        A.fop_T,
        A.fop
        )

# adjoint(judiPDEfull)
adjoint(A::judiPDEfull{DDT,RDT}) where {DDT,RDT} =
    judiPDEfull{DDT,RDT}("Proj*F'*Proj'",A.n,A.m,A.info,A.model,A.srcGeometry,A.recGeometry,A.options,
        A.fop_T,
        A.fop
        )

############################################################
## overloaded Base *(...judiPDEfull...)

# *(judiPDEfull,vec)
function *(A::judiPDEfull{ADDT,ARDT}, v::AbstractVector{vDT}) where {ADDT,ARDT,vDT}
    A.n == size(v,1) || throw(judiPDEfullException("Shape mismatch: A:$(size(A)), v: $(size(v))"))
    jo_check_type_match(ADDT,vDT,join(["DDT for *(judiPDEfull,judiVector):",A.name,typeof(A),vDT]," / "))
    V = A.fop(v)
    jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiPDEfull,judiVector):",A.name,typeof(A),eltype(V)]," / "))
    return V
end

# *(judiPDEfull,judiVector)
function *(A::judiPDEfull{ADDT,ARDT}, v::judiVector{vDT, AT}) where {ADDT, ARDT,vDT, AT}
    A.n == size(v,1) || throw(judiPDEfullException("Shape mismatch: A:$(size(A)), v: $(size(v))"))
    if compareGeometry(A.srcGeometry,v.geometry) == false && compareGeometry(A.recGeometry,v.geometry) == false
        throw(judiPDEfullException("geometry mismatch"))
    end
    jo_check_type_match(ADDT,vDT,join(["DDT for *(judiPDEfull,judiVector):",A.name,typeof(A),vDT]," / "))
    V = A.fop(v)
    jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiPDEfull,judiVector):",A.name,typeof(A),eltype(V)]," / "))
    return V
end

# *(num,judiPDEfull)
function *(a::Number,A::judiPDEfull{ADDT,ARDT}) where {ADDT,ARDT}
    return judiPDEfull{ADDT,ARDT}("(N*"*A.name*")",A.m,A.n,A.info,A.model,A.srcGeometry,A.recGeometry,A.options,
        v1 -> jo_convert(ARDT,a, false)*A.fop(v1),
        v2 -> jo_convert(ADDT,a, false)*A.fop_T(v2)
        )
end

############################################################
## overloaded Bases +(...judiPDEfull...), -(...judiPDEfull...)
# -(judiPDEfull)
-(A::judiPDEfull{DDT,RDT}) where {DDT,RDT} =
    judiPDEfull{DDT,RDT}("(-"*A.name*")",A.m,A.n,A.info,A.model,A.srcGeometry,A.recGeometry,A.options,
        v1 -> -A.fop(v1),
        v2 -> -A.fop_T(v2)
        )

############################################################
## Additional overloaded functions

# Subsample Modeling operator
function subsample(F::judiPDEfull{ADDT,ARDT}, srcnum) where {ADDT,ARDT}

    srcGeometry = subsample(F.srcGeometry,srcnum)       # Geometry of subsampled data container
    recGeometry = subsample(F.recGeometry,srcnum)

    info = Info(F.info.n, length(srcnum), F.info.nt[srcnum])
    return judiModeling(info, F.model, srcGeometry, recGeometry; options=F.options, DDT=ADDT, RDT=ARDT)
end

getindex(F::judiPDEfull,a) = subsample(F,a)

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
using Test
@testset "Generic functions" begin
    import Random: randn, randexp, bitrand, randstring, randsubseq, shuffle, randperm, randcycle

    @info "Testing Generic functions"
    r = PCG.PCGStateOneseq(123)
    @test randn(r) == -0.4957408739873887
    @test randn(r, Float64, 8) == [
        -0.38659391873108484, -0.12391575881839223, 0.305709506865529, 1.0234147128314572,
        1.1925044460767074, 0.8448640519165571, -0.2982275988025108, 0.9615482011555472
    ]
    @test randexp(r) == 1.0279939998988223
    @test bitrand(r, 3, 3) == [false false false; false true false; false true true]
    # randstring method changed after a julia release
    randstringres = ""
    if VERSION >= v"1.5"
        randstringres = "6j22eD3r"
    else
        randstringres = "KuSFMMEc"
    end
    @test randstring(r) == randstringres
    a = 1:100
    @test randsubseq(r, a, 0.1) == [1, 11, 12, 17, 19, 21, 27, 38, 44, 54, 58, 78, 80]
    a = 1:10
    @test shuffle(r, a) == [3, 7, 5, 10, 2, 6, 1, 4, 9, 8]
    @test randperm(r, 10) == [2, 10, 1, 6, 4, 3, 7, 8, 9, 5]
    @test randcycle(r, 10) == [8, 4, 5, 1, 10, 2, 3, 9, 7, 6]
end


text:
@testset "Test fast_mul function" begin
    n_off = 200
    n_diag = 50
    Ne = 500

    ψ_off = randn(Ne, n_off)
    ψ_diag = randn(Ne, n_off)
    λ = 0.1
    ψ_aug = vcat(ψ_off, Matrix(√λ*I, n_off, n_off))
    F = qr(ψ_aug)
    Q1 = Matrix(F.Q)[1:Ne,:]

    out = fast_mul(deepcopy(ψ_diag), F.Q, Ne, n_off)

    @test norm(out - Q1*Q1'*ψ_diag)<1e-10
end

@testset "Test fast_mul2 function" begin
    n_off = 200
    n_diag = 50
    Ne = 500

    ψ_off = randn(Ne, n_off)
    ψ_diag = randn(Ne, n_off)
    F = qr(ψ_off)
    Q1 = Matrix(F.Q)[1:Ne,:]

    out = fast_mul2(deepcopy(ψ_diag), F.Q, Ne, n_off)

    @test norm(out - Q1*Q1'*ψ_diag)<1e-10
end

@testset "Test optimization for Nx=1 and p=0" begin
   Nx = 1
   Ne = 10
   p = 0
   γ = 2.0

   X = reshape([-1.5438
      -1.5518
       0.8671
      -0.1454
      -0.3862
       1.3162
      -0.7965
       0.1354
       0.4178
       0.8199],1,Ne)

   S = RadialMap(Nx, p; γ = γ);
   W = create_weights(S, X)
   center_std!(S, X);
   compute_weights(S, X, W)

   # With λ = δ = 0.0

   λ = 0.0
   δ = 0.0
   S.C[1] = optimize(S.C[1], W, λ, δ);

   @test norm(extractcoeff(S.C[1]) - [0.092087627168250; 1.061773632748180])<1e-10

   # With λ = 0.1 and δ = 1e-8

   λ = 0.1
   δ = 1e-8
   S.C[1] = optimize(S.C[1], W, λ, δ);

   @test norm(extractcoeff(S.C[1]) - [ 0.092087627168250; 1.061773632748180])<1e-10

end

@testset "Test optimization for Nx=2 and p=0" begin

    Nx = 2
    Ne = 10
    p = 0
    γ = 2.0
    λ = 0.0

    X = Matrix([-0.8544   -0.7938;0.3573    0.5410; 2.7485   -0.5591; -1.5130    1.9766; 0.4340    0.5447; -0.2298   -0.1379;-0.8271   0.6199;-0.8320   -0.0056;0.4979    1.1072;2.3156   -0.1856]')

     S = RadialMap(Nx, p; γ = γ)

    center_std!(S, X)
    W = create_weights(S, X)
   compute_weights(S, X, W)

    # With λ = δ = 0.0

    λ = 0.0
    δ = 0.0

    S.C[2] = optimize(S.C[2], W, λ, δ)

    @test norm(extractcoeff(S.C[2]) -[0.357284296467942; -0.513627363149479; 1.411806803694895 ])<1e-10

    # With λ = 0.1  & δ = 1e-8
    λ = 0.1
    δ = 1e-8

    S.C[2] = optimize(S.C[2], W, λ, δ)

    @test norm(extractcoeff(S.C[2]) -[0.353742957402480; -0.512879943932609; 1.411791355362390 ])<1e-10
end

@testset "Test optimization for Nx=3 and p=0" begin

   Nx = 3
   Ne = 10
   p = 0
   γ = 2.0

   X = Matrix([0.1733    0.2760    1.0093; 0.6163    0.4598  0.0510;0.8639   -0.2175   -0.4404;-1.4169    0.7961   -0.8485;0.2000   -1.5180   -0.2404;0.1958   -1.0749    0.6029;-0.1511   -3.0722   -1.5163;0.8928    0.5214   -0.0683;-0.3642   -0.9910    0.7824;-0.8262   -0.2531   -1.4207]')
   S = RadialMap(Nx, p; γ = γ)
   center_std!(S, X)
   W = create_weights(S, X)
  compute_weights(S, X, W)

   # With λ = δ = 0.0
   λ = 0.0
   δ = 0.0

   S.C[3] = optimize(S.C[3], W, λ, δ)

   @test norm(extractcoeff(S.C[3]) -[-0.609710825728394; -0.275729138896713; 0.156805424543943; 1.366651306860556])<1e-10

   # With λ = 0.1  & δ = 1e-8
   λ = 0.1
   δ = 1e-8

   S.C[3] = optimize(S.C[3], W, λ, δ)

   @test norm(extractcoeff(S.C[3]) -[-0.603717368357810; -0.273039744501052; 0.158055928017947; 1.366632953261741])<1e-10

end

@testset "Test optimization for Nx=1 and p=3" begin
   Nx = 1
   Ne = 10
   p = 3
   γ = 2.0

   X = reshape([-1.5438
      -1.5518
       0.8671
      -0.1454
      -0.3862
       1.3162
      -0.7965
       0.1354
       0.4178
       0.8199],1,Ne)

   S = RadialMap(Nx, p; γ = γ)
   W = create_weights(S, X)
   center_std!(S, X)
  compute_weights(S, X, W)

   # With λ = δ = 0.0
   λ = 0.0
   δ = 0.0
   S.C[1] = optimize(S.C[1], W, λ, δ);

   @test norm(extractcoeff(S.C[1]) - [0.453363732161572; 2.304194295945015; 0.0; 0.0; 0.0;2.552281803302206])<1e-10

   # With λ = 0.1 and δ = 1e-8
   λ = 0.1
   δ = 1e-8
   S.C[1] = optimize(S.C[1], W, λ, δ);

   @test norm(extractcoeff(S.C[1]) - [0.450882345511096; 2.295315495868833; 0.000000043115815; 0.000000032750391; 0.000000035316480; 2.547661799360081])<1e-10

end

@testset "Test optimization for Nx=2 and p=2" begin

    Nx = 2
    Ne = 10
    p = 2
    γ = 2.0
    λ = 0.0

    X = Matrix([-0.8544   -0.7938;0.3573    0.5410;
                2.7485   -0.5591; -1.5130    1.9766;
                0.4340    0.5447; -0.2298   -0.1379;
                 -0.8271   0.6199;-0.8320   -0.0056;
                  0.4979    1.1072;2.3156   -0.1856]')

    S = RadialMap(Nx, p; γ = γ)

    center_std!(S, X);

    W = create_weights(S, X)
   compute_weights(S, X, W)
    λ = 0.0
    δ = 0.0

    ψ_off, ψ_mono, dψ_mono = rearrange(W,2)

    @test norm(ψ_off -[-0.854400000000000   0.156738618846197   0.137603564500746;
       0.357300000000000   0.140640215047531   0.156654816189471;
       2.748500000000000   0.058401570295014   0.104057315307723;
      -1.513000000000000   0.151170210915668   0.116608188942393;
       0.434000000000000   0.138617339050529   0.156745732145583;
      -0.229800000000000   0.152477989231576   0.151338989275322;
      -0.827100000000000   0.156747096365829   0.138351021862248;
      -0.832000000000000   0.156746902701754   0.138217735507039;
       0.497900000000000   0.136859341544330   0.156712800799488;
       2.315600000000000   0.073104270826721   0.119634493493623])<1e-10

    @test norm(ψ_mono -[-0.507201580119279   0.215069670075122   0.076485784936572   0.000988033380379;
      -0.016861853328784   0.748419931509968   0.499210902664930   0.126296856325961;
      -0.344698250357007   0.297187601194573   0.119362800788456   0.002977317564677;
      -0.000006458746850   0.987394188705160   0.937353523618463   1.123806421603995;
      -0.016625169238807   0.749707195139661   0.500789097335070   0.127438446122333;
      -0.140628787073782   0.471189825582940   0.233355359282062   0.016320116571698;
      -0.012393614690139   0.775096387582616   0.532827724212160   0.152299949965802;
      -0.099762356857753   0.528810174417060   0.278806935049185   0.025944471088712;
      -0.001358786238840   0.901159801018090   0.726875126285091   0.395949949965802;
      -0.157948250357007   0.450499667316035   0.218036687265869   0.013696473624777])<1e-10

      @test norm(dψ_mono -[0.758353059754083   0.319323853074808   0.153626107072438   0.004883698987942
       0.064354624108671   0.348381106647456   0.426538543776643   0.307522804512740
       0.621958562984279   0.378284965307517   0.213076838459543   0.013276098094082
       0.000046704774842   0.035615595830791   0.131731605768779   0.958660103131430
       0.063583925230259   0.347436551236228   0.426538543776643   0.309553737400797
       0.348279343219786   0.434770043962370   0.327283964601450   0.059220575922185
       0.049425557686157   0.327622035258933   0.425094484357949   0.352020159021600
       0.270920987266371   0.434770043962370   0.359163756509705   0.087695615816903
       0.006927621078729   0.190131276118064   0.355549818399670   0.647979840978400
       0.378041437015721   0.432547368859386   0.314948152204838   0.050949232751226])<1e-10

    # Rescale variable

    ψ_off, ψ_mono, dψ_mono = rearrange(W,2)

    no = size(ψ_off,2)
    nd = size(ψ_mono,2)
    nx = size(ψ_off,2)+size(ψ_mono,2)+1
    nlog = size(dψ_mono,2)

    x = zeros(nx)

    μψ = deepcopy(mean(ψ_mono, dims=1)[1,:])
    σψ = deepcopy(std(ψ_mono, dims=1, corrected=false)[1,:])

    ψ_mono .-= μψ'
    ψ_mono ./= σψ'
    dψ_mono ./= σψ'

    A = ψ_mono'*ψ_mono
    rmul!(A,1/Ne)

    μψ_off = deepcopy(mean(ψ_off, dims = 1)[1,:])
    σψ_off = deepcopy(std(ψ_off, dims = 1, corrected = false)[1,:])
    ψ_off .-= μψ_off'
    ψ_off ./= σψ_off'

    # Verify that ψ_off, ψ_mono and dψ_mono are correctly rescaled
    @test norm(ψ_off -[ -0.804643054198426   0.720356396590189   0.000623052929297;
       0.111611046705843   0.248725359467419   1.070112263485765;
       1.919770497132754  -2.160600297987687  -1.882578752010958;
      -1.302658198917046   0.557220098094831  -1.178004367363986;
       0.169609470027917   0.189461651315464   1.075216055914641;
      -0.332337771187114   0.595533762677487   0.771695150386939;
      -0.783999547592264   0.720604760429600   0.042583423646065;
      -0.787704792367729   0.720599086700126   0.035101063896698;
       0.217928886589593   0.137958024156060   1.073367373218454;
       1.592423463806471  -1.729858841443491  -1.008115264102931])<1e-10

    @test norm(ψ_mono -[-2.326960601058417  -1.624755258466316  -1.277802709623370  -0.600690875072745;
       0.695935005892621   0.555916671737632   0.330653865772685  -0.219729086471929;
      -1.325143748753542  -1.289005413932899  -1.114656931934708  -0.594643086552582;
       0.799846835694096   1.532994011692314   1.997773585133276   2.812882817478138;
       0.697394139747648   0.561179816554494   0.336658849834394  -0.216258444258148;
      -0.067075820810003  -0.577574682660793  -0.680918616973539  -0.554078532076930;
       0.723481252328697   0.664986832870473   0.458564877903041  -0.140674915908932;
       0.184861653602686  -0.341986390319673  -0.507976714922384  -0.524818730374682;
       0.791509871317472   1.180413498950799   1.196909505920374   0.600065740739791;
      -0.173848587961260  -0.662169086426033  -0.739205711109769  -0.562054887501982])<1e-10

    @test norm(dψ_mono -[4.675171126462628   1.305597115808049   0.584542795566616   0.014847340017250;
       0.396739851731833   1.424401477250050   1.622966548768615   0.934925689100038;
       3.834312630667317   1.546667293731231   0.810751070876959   0.040361771475244;
       0.000287930287880   0.145618996973497   0.501234865401183   2.914502418618957;
       0.391988569862249   1.420539539569087   1.622966548768615   0.941100097313083;
       2.147107483013941   1.777613886251173   1.245305809396379   0.180041404866467;
       0.304703643285317   1.339525313221805   1.617471945373052   1.070205802369169;
       1.670200918858584   1.777613886251173   1.366607475104451   0.266610745107435;
       0.042708094370782   0.777376457581746   1.352856547439218   1.969977479587423;
       2.330587828729104   1.768526189933186   1.198368407926584   0.154895005638792])<1e-10

    # Verify that we get the same QR decomposition
    #Assemble reduced QR to solve least square problem
    F = qr(vcat(ψ_off, Matrix(√λ*I, no, no)))
    Q1 = F.Q[1:Ne,1:no]
    # Asqrt = ψ_mono - Q1*(Q1'*ψ_mono)
    Asqrt = ψ_mono - fast_mul(ψ_mono, F.Q, Ne, no)

    # A = (Asqrt*Asqrt')/Ne
    A = BLAS.gemm('T', 'N', 1/Ne, Asqrt, Asqrt)

    @test norm(Q1 - [ -0.254450475470134   0.038542438904502  -0.287807026572883;
       0.035294511962590  -0.342003510157944   0.197040760785064;
       0.607084735573325   0.334774229448390  -0.420694384780400;
      -0.411936692127055   0.650442320408613   0.454822398694172;
       0.053635223802228  -0.337693395873161   0.262404719972692;
      -0.105094430945516  -0.271842813033772  -0.232724768964479;
      -0.247922425493313   0.019468149870942  -0.298845707008466;
      -0.249094126771204   0.022854173498846  -0.296968816757815;
       0.068915164956764  -0.332045729496864   0.316658455434137;
       0.503568514512315   0.217504136430448   0.306114369197955])<1e-10

    @test norm(Asqrt - [ -1.950781401302997  -1.300838981851897  -0.991890513984422  -0.458060391024002;
       0.099892510115221   0.085352068075599   0.033990870033547  -0.016912053174583;
       0.149656589279390   0.450550985993256   0.548653802955680   0.551333139636105;
       0.144073370024957   0.204556328914631   0.251633283508030   0.327703732540707;
      -0.003122789715729  -0.027881050368865  -0.081782795321125  -0.121731102269892;
       0.083258235300422  -0.268551915835203  -0.230658522708182   0.178298456149084;
       1.110570901959247   1.016081154773970   0.783074186747397   0.062561590138456;
       0.570158295426148   0.004447219686410  -0.190152589533267  -0.332197943172162;
       0.005871968596686   0.493585572789772   0.676896851023931   0.601539422000769;
      -0.209577679683306  -0.657301382177627  -0.799764572721539  -0.792534850824440])<1e-10

    @test norm(A - [0.546802799943131   0.389242527702498   0.297060133111796   0.108653878210649;
       0.389242527702498   0.372531651490579   0.331069886098672   0.174530422370543;
       0.297060133111796   0.331069886098672   0.315641016351602   0.196073264491252;
       0.108653878210649   0.174530422370543   0.196073264491252   0.177230269927555])<1e-10

    # Solve nonlinear optimization problem for g_mono
    g_mono = ones(nd)
    # Construct loss function, the convention is flipped inside the solver
    lhd = LHD(A, dψ_mono, λ, δ)
    g_mono,_,_,_ = projected_newton(g_mono, lhd, "TrueHessian")

    g_mono .+= δ

    @test norm(g_mono - [0.586712518793935; 0.421099709077361; 0; 1.154877772114433] )<1e-10

    # Approximate off-diagonal components
    @test norm(F.R -[3.162277660168379  -2.987894843243091  -1.096427127338353;
                      0  -1.035608229844348  -2.909339089462357;
                      0                   0   0.577575637438866])<1e-10

    g_off = -F.R\(Q1'*(ψ_mono*g_mono))

    @test norm(g_off -   [18.950658982438039; 21.898788492005636;-7.359525277156965])<1e-10

    g_mono ./= σψ
    g_off ./= σψ_off

    @test norm(g_mono - [3.617024276646359; 1.721720943628389; 0; 3.511040095485572] )<1e-10

    @test norm(g_off - 100*[ 0.143299653441423; 6.415635026543544; -4.131451837044035])<1e-10

    x[no+1] = -dot(μψ, g_mono)-dot(μψ_off, g_off)

    @test norm(x[no+1] -(-32.224538086041441))<1e-10

    x[no+2:nx] .= deepcopy(g_mono)

    if Nx>1
    x[1:no] .= deepcopy(g_off)
    end

    @test norm(x - [  14.32996534414215; 641.5635026543499; -413.14518370440067; -32.22453808604119; 3.61702427664633; 1.721720943628398; 0.0; 3.511040095485544])<1e-10

     S.C[2] = optimize(S.C[2], W, λ, δ);
     @test norm(extractcoeff(S.C[2]) - [  14.32996534414215; 641.5635026543499; -413.14518370440067; -32.22453808604119; 3.61702427664633; 1.721720943628398; 0.0; 3.511040095485544])<1e-10

       # With λ=0.1 and δ =1e-8
       λ = 0.1
       δ = 1e-8
       S.C[2] = optimize(S.C[2], W, λ, δ);
       @test norm(extractcoeff(S.C[2]) - [  1.765903393601119; 67.455814562648044; -48.756672231027366;  -3.584068900212547; 2.685855024209524; 1.468924592466743; 0.477896050052401; 1.308323191869906])<1e-10

 end

 @testset "Test optimization for Nx=2 and p=3" begin
     Nx = 2
     Ne = 10
     p = 3
     γ = 2.0

     X = Matrix([ -0.8544   -0.7938;0.3573    0.5410;
       2.7485   -0.5591; -1.5130    1.9766;
       0.4340    0.5447; -0.2298   -0.1379;
        -0.8271   0.6199;-0.8320   -0.0056;
         0.4979    1.1072;2.3156   -0.1856]')

     S = RadialMap(Nx, p; γ = γ)

     center_std!(S, X);

     W = create_weights(S, X)
    compute_weights(S, X, W)

     # With λ = δ = 0.0
     λ = 0.0
     δ = 0.0

     S.C[2] = optimize(S.C[2], W, λ, δ);
     @test norm(extractcoeff(S.C[2]) - [100*0.090840649611966; 100*1.373434239444840; 100*0.414826633064425; 100*(-0.287690260374867); -29.015187078639887; 4.209897756012423; 0.108369156771108; 0; 0; 8.850495668947518])<1e-10

     # With λ = 0.1 and δ = 1e-8
     λ = 0.1
     δ = 1e-8

     xopt = optimize(S.C[2], W, λ, δ);
     @test  norm(extractcoeff(S.C[2]) - [3.855369324206674; 28.717863158966477; 49.828564728403187; -23.158217471101509; -12.460695228234130; 3.669743272444025; 0.000000038099920; 0.794735311810544; 0.486353838959244; 5.376840166095511])<1e-10
 end

 @testset "Test optimization for Nx=3 and p=3" begin

     Nx = 3
     Ne = 10
     p = 3
     γ = 2.0

     X = Matrix([0.1733    0.2760    1.0093;
                   0.6163    0.4598  0.0510;
                 0.8639   -0.2175   -0.4404;
                -1.4169    0.7961   -0.8485;
                 0.2000   -1.5180   -0.2404;
                 0.1958   -1.0749    0.6029;
                -0.1511   -3.0722   -1.5163;
                 0.8928    0.5214   -0.0683;
                -0.3642   -0.9910    0.7824;
                -0.8262   -0.2531   -1.4207]')

      S = RadialMap(Nx, p; γ = γ)

     center_std!(S, X);
     W = create_weights(S, X)
    compute_weights(S, X, W)

     # With λ = δ = 0.0
     λ = 0.0
     δ = 0.0

     S.C[3] = optimize(S.C[3], W, λ, δ);
     x_matlab  = [0.002236347046586e5;
     0.825719944026571e5    ;
     -1.071956974125768e5    ;
     0.593296179017564e5    ;
     -0.011049031761152e5    ;
     0.833628110516284e5      ;
     -1.085535253851930e5       ;
     0.696061056463158e5      ;
     -1.562689086584463e+04 ;
           0.0             ;
           0.0             ;
           0.0             ;
           0.0             ;
           1000*2.263679644184734  ]
     @test norm(extractcoeff(S.C[3]) - x_matlab)/norm(x_matlab)<1e-10

     # With λ = 0.1 and δ = 1e-8
     λ = 0.1
     δ = 1e-8
     #
     S.C[3] = optimize(S.C[3], W, λ, δ);

     @test norm(extractcoeff(S.C[3]) - [   1.051136807671425;
     0.638509720542414;
     -23.285653280072371;
     -0.955334980052370;
     -2.821327423813651;
     -15.674943048773969;
     5.749217387570198;
     13.386252213401695;
     6.134578583234069;
     5.778639058642066;
     0.000000041429348;
     0.000000029234255;
     0.000000036290793;
     3.914178426750072])<1e-10
 end

## Test optimization with SparseRadialMap

@testset "Test optimization with SparseRadialMap Nx=1 and p=0" begin

    Nx = 1
    Ne = 10
    γ = 2.0

    X  =  reshape([0.017374310426544
       0.735136768751132
      -0.984272698022620
       0.134246840483926
      -0.720332223266792
       0.570079159353393
      -0.178637174884214
      -0.567524930231864
       0.028803259857213
      -0.614918061442867],1,Ne)

      S = SparseRadialMap(Nx, [0]; γ = γ)
      center_std!(S, X)

      S.C[1], error = optimize(S.C[1], X, nothing, S.λ, S.δ)

    @test norm(extractcoeff(S.C[1])- [0.295725384806749; 1.871626642209824])<1e-14
end

@testset "Test optimization with SparseRadialMap Nx=1 and p=2" begin

    Nx = 1
    Ne = 10
    γ = 2.0

    X = reshape([0.017374310426544
       0.735136768751132
      -0.984272698022620
       0.134246840483926
      -0.720332223266792
       0.570079159353393
      -0.178637174884214
      -0.567524930231864
       0.028803259857213
      -0.614918061442867],1,Ne)

      S = SparseRadialMap(Nx, [2]; γ = γ)
      center_std!(S, X)

      S.C[1], _ = optimize(S.C[1], X, nothing, S.λ, S.δ)

    @test norm(extractcoeff(S.C[1]) - [0.033932167701879; 4.133761421140692; 0.000000039870398; 0.000000040266888; 3.328673620931276])<1e-14
end

@testset "Test optimization with SparseRadialMap Nx=2 and p=[0 0]" begin

    Nx = 2
    Ne = 10
    γ = 2.0

    X = Matrix([  -0.737609918489704   0.413557038884469
       0.602049917255313   2.288861842327617
       1.099780136322616   0.087381084349074
       0.200615477691937  -1.421682449747359
      -0.427674259885219  -2.022946349735899
      -0.087812697126117  -0.198683683159131
       0.135174851422381   0.384267545498855
       0.753036725634817   0.470902988284163
      -0.228309526653063   1.272982451997852
      -0.684323467801375   0.567990940582801]')

      S = SparseRadialMap(Nx, [[-1], [0; 0]])
      center_std!(S, X)

      S.C[2], error = optimize(S.C[2], X, nothing, S.λ, S.δ)

    @test norm(extractcoeff(S.C[2]) - [-0.350857378599181; -0.139765206111085; 0.877501808357206])<1e-14
end

@testset "Test optimization with SparseRadialMap Nx=2 and p=[0 1]" begin
    Nx = 2
    Ne = 10
    γ = 2.0

    X = Matrix([  -0.737609918489704   0.413557038884469
           0.602049917255313   2.288861842327617
           1.099780136322616   0.087381084349074
           0.200615477691937  -1.421682449747359
          -0.427674259885219  -2.022946349735899
          -0.087812697126117  -0.198683683159131
           0.135174851422381   0.384267545498855
           0.753036725634817   0.470902988284163
          -0.228309526653063   1.272982451997852
          -0.684323467801375   0.567990940582801]')

    S = SparseRadialMap(Nx, [[-1], [0; 1]]; γ = γ)
    center_std!(S, X)

    S.C[2], _ = optimize(S.C[2], X, nothing, S.λ, S.δ)

    @test norm(extractcoeff(S.C[2]) - [-0.296920845324240; -0.450160011091217; 0.685091255311441;
    1.423017573574061;
    0.566245516556744])<1e-14
end

@testset "Test optimization with SparseRadialMap Nx=2 and p=[0 2]" begin

Nx = 2
Ne = 10
γ = 2.0

X = Matrix([  -0.737609918489704   0.413557038884469
   0.602049917255313   2.288861842327617
   1.099780136322616   0.087381084349074
   0.200615477691937  -1.421682449747359
  -0.427674259885219  -2.022946349735899
  -0.087812697126117  -0.198683683159131
   0.135174851422381   0.384267545498855
   0.753036725634817   0.470902988284163
  -0.228309526653063   1.272982451997852
  -0.684323467801375   0.567990940582801]')

  S = SparseRadialMap(Nx, [[-1], [0; 2]]; γ = γ)
  center_std!(S, X)

  S.C[2], _ = optimize(S.C[2], X, nothing, S.λ, S.δ)

@test norm(extractcoeff(S.C[2]) - [-0.290599917218335; -0.244106257901270; 0.928468014306679;
   0.000000037435467;
   1.275553022078607;
   0.695142104297181])<1e-14
end

@testset "Test optimization with SparseRadialMap Nx=2 and p=[2 2]" begin

Nx = 2
Ne = 10
γ = 2.0

X = Matrix([  -0.737609918489704   0.413557038884469
   0.602049917255313   2.288861842327617
   1.099780136322616   0.087381084349074
   0.200615477691937  -1.421682449747359
  -0.427674259885219  -2.022946349735899
  -0.087812697126117  -0.198683683159131
   0.135174851422381   0.384267545498855
   0.753036725634817   0.470902988284163
  -0.228309526653063   1.272982451997852
  -0.684323467801375   0.567990940582801]')

  S = SparseRadialMap(Nx, [[-1], [2; 2]]; γ = γ)
  center_std!(S, X)

  S.C[2], _ = optimize(S.C[2], X, nothing, S.λ, S.δ)

@test norm(extractcoeff(S.C[2]) - [   2.205922038593491;
  21.094127932392510;
 -15.821105702007772;
  -2.010486812456307;
   0.999304051913047;
   0.000000037435467;
   1.303627375210608;
   0.764098525416240])<1e-12
end

@testset "Test optimization with SparseRadialMap Nx=2 and p=[2 0]" begin

Nx = 2
Ne = 10
γ = 2.0

X = Matrix([  -0.737609918489704   0.413557038884469
   0.602049917255313   2.288861842327617
   1.099780136322616   0.087381084349074
   0.200615477691937  -1.421682449747359
  -0.427674259885219  -2.022946349735899
  -0.087812697126117  -0.198683683159131
   0.135174851422381   0.384267545498855
   0.753036725634817   0.470902988284163
  -0.228309526653063   1.272982451997852
  -0.684323467801375   0.567990940582801]')

  S = SparseRadialMap(Nx, [[-1], [2; 0]]; γ = γ)
  center_std!(S, X)

  S.C[2], _ = optimize(S.C[2], X, nothing, S.λ, S.δ)

  @test norm(extractcoeff(S.C[2]) - [ 2.175874123786137; 21.548828045380638; -15.447825995278832;
                     -2.188592622357425;
                      0.933194195948753])<1e-12
end

@testset "Test optimization with SparseRadialMap Nx=3 and p=[0 0 0]" begin

Nx = 3
Ne = 10
γ = 2.0

X = Matrix([  0.375413787969851  -0.991365932946799   2.278410854329425
   1.121729788801830   0.216910057645809   1.225483370520654
  -0.015678780143914  -0.433156720743754   0.831025086529337
  -0.674285107908432   0.572290138313797  -1.493962274420239
   2.973707136592813   0.264962909286999   0.967684185813051
  -0.369888286517882   0.515287247888629   0.671466599596567
   0.435525083368351  -1.123591346893271  -1.469294876880461
  -1.130590197154928  -0.559193704135775  -0.001325508177795
  -2.003193706192794  -0.846920007816157   0.559462142237872
   0.277688627801369   1.179218444917551   0.729724225819224]')

  S = SparseRadialMap(Nx, [[-1], [-1; -1], [0; 0; 0]]; γ = γ)
  center_std!(S, X)

  S.C[3], _ = optimize(S.C[3], X, nothing, S.λ, S.δ)

  @test norm(extractcoeff(S.C[3]) - [-0.238818177855404;  0.112131135966394; -0.369644592955786;  0.946375346756640])<1e-12

end

@testset "Test optimization with SparseRadialMap Nx=3 and p=[0 -1 0]" begin

Nx = 3
Ne = 10
γ = 2.0

X = Matrix([  0.375413787969851  -0.991365932946799   2.278410854329425
   1.121729788801830   0.216910057645809   1.225483370520654
  -0.015678780143914  -0.433156720743754   0.831025086529337
  -0.674285107908432   0.572290138313797  -1.493962274420239
   2.973707136592813   0.264962909286999   0.967684185813051
  -0.369888286517882   0.515287247888629   0.671466599596567
   0.435525083368351  -1.123591346893271  -1.469294876880461
  -1.130590197154928  -0.559193704135775  -0.001325508177795
  -2.003193706192794  -0.846920007816157   0.559462142237872
   0.277688627801369   1.179218444917551   0.729724225819224]')

  S = SparseRadialMap(Nx, [[-1], [-1; -1], [0; -1; 0]]; γ = γ)
  center_std!(S, X)

  S.C[3], _ = optimize(S.C[3], X, nothing, S.λ, S.δ)

  @test norm(extractcoeff(S.C[3]) - [-0.221895676794685; -0.383494876901790; 0.943249178908141])<1e-12

end

@testset "Test optimization with SparseRadialMap Nx=3 and p=[2 -1 0]" begin

Nx = 3
Ne = 10
γ = 2.0

X = Matrix([  0.375413787969851  -0.991365932946799   2.278410854329425
   1.121729788801830   0.216910057645809   1.225483370520654
  -0.015678780143914  -0.433156720743754   0.831025086529337
  -0.674285107908432   0.572290138313797  -1.493962274420239
   2.973707136592813   0.264962909286999   0.967684185813051
  -0.369888286517882   0.515287247888629   0.671466599596567
   0.435525083368351  -1.123591346893271  -1.469294876880461
  -1.130590197154928  -0.559193704135775  -0.001325508177795
  -2.003193706192794  -0.846920007816157   0.559462142237872
   0.277688627801369   1.179218444917551   0.729724225819224]')

  S = SparseRadialMap(Nx, [[-1], [-1; -1], [2; -1; 0]]; γ = γ)
  center_std!(S, X)

  S.C[3], _ = optimize(S.C[3], X, nothing, S.λ, S.δ)

  @test norm(extractcoeff(S.C[3]) - [   0.262639648863806;
                       16.549041782004032;
                      -11.747435767576015;
                       -1.410877218960467;
                        0.996852513304659])<1e-12

end

@testset "Test optimization with SparseRadialMap Nx=3 and p=[2 -1 2]" begin

Nx = 3
Ne = 10
γ = 2.0

X = Matrix([  0.375413787969851  -0.991365932946799   2.278410854329425
   1.121729788801830   0.216910057645809   1.225483370520654
  -0.015678780143914  -0.433156720743754   0.831025086529337
  -0.674285107908432   0.572290138313797  -1.493962274420239
   2.973707136592813   0.264962909286999   0.967684185813051
  -0.369888286517882   0.515287247888629   0.671466599596567
   0.435525083368351  -1.123591346893271  -1.469294876880461
  -1.130590197154928  -0.559193704135775  -0.001325508177795
  -2.003193706192794  -0.846920007816157   0.559462142237872
   0.277688627801369   1.179218444917551   0.729724225819224]')

  S = SparseRadialMap(Nx, [[-1], [-1; -1], [2; -1; 2]]; γ = γ)
  center_std!(S, X)

  S.C[3], _  = optimize(S.C[3], X, nothing, S.λ, S.δ)

  @test norm(extractcoeff(S.C[3]) - [     0.246247362219524;
                         20.601841748108630;
                        -16.322738975189164;
                         -1.294380861832026;
                          0.898035319369626;
                          0.000000043296784;
                          2.067302695773216;
                          0.631236879169851])<1e-12

end

@testset "Test optimization with SparseRadialMap Nx=3 and p=[-1 -1 2]" begin

Nx = 3
Ne = 10
γ = 2.0

X = Matrix([  0.375413787969851  -0.991365932946799   2.278410854329425
   1.121729788801830   0.216910057645809   1.225483370520654
  -0.015678780143914  -0.433156720743754   0.831025086529337
  -0.674285107908432   0.572290138313797  -1.493962274420239
   2.973707136592813   0.264962909286999   0.967684185813051
  -0.369888286517882   0.515287247888629   0.671466599596567
   0.435525083368351  -1.123591346893271  -1.469294876880461
  -1.130590197154928  -0.559193704135775  -0.001325508177795
  -2.003193706192794  -0.846920007816157   0.559462142237872
   0.277688627801369   1.179218444917551   0.729724225819224]')

  S = SparseRadialMap(Nx, [[-1], [-1; -1], [-1; -1; 2]]; γ = γ)
  center_std!(S, X)

  S.C[3], _  = optimize(S.C[3], X, nothing, S.λ, S.δ)

  @test norm(extractcoeff(S.C[3]) - [       -0.1178422702775955;
                             0.9302021941755624;
                          4.3296784157577845e-8;
                             1.4379052718065835;
                             0.6286862437766934])<1e-12

end

@testset "SparseRadialMap Test with optimization and evaluation of the resulting map" begin

      X = Matrix([6.831108232125667   6.831108465893829  10.748176564682273  17.220877261555913;
      7.341399771164814   7.341399033871628  11.440770534993137  18.015116463525178;
      6.877822405336048   6.877822863906546  10.793381797612751  17.379545045967195;
      6.971116690210811   6.971117274254750  10.941318835214862  17.571340756941474;
      7.800189373932513   7.800189247976281  11.929084958885008  18.852191934527564;
      6.706165535739543   6.706164938380805  10.596868885136731  17.100665224055952;
      7.162689725595848   7.162690676866085  11.236535720904731  17.691824497167673;
      7.700262639909202   7.700262776963255  11.870709137884843  18.589199307190668;
      6.941836510989174   6.941836045832242  10.822143995520504  17.705463568541791;
      7.284400905037024   7.284402833890172  11.438244192350945  17.713180862215779;
      6.114110264862537   6.114110399157427   9.744642327011279  16.427244631267527;
      6.884587321648968   6.884587349033709  10.880879551471059  17.360443337222691;
      6.964593674795581   6.964592742775720  11.015485574206627  17.379079978320874;
      7.174466645033212   7.174468203234895  11.165089908005559  17.802087293640490;
      7.625148752606352   7.625148374481078  11.784628799983244  18.470370177818385;
      7.878948627898106   7.878948112320454  12.118366730171301  18.770876038639738;
      7.406860218451342   7.406861450644357  11.470902695133232  18.293040279245396;
      7.308212626757207   7.308212927902813  11.381148951581917  17.962547003679237;
      7.519983057000856   7.519981781231936  11.597728250168640  18.407564589808729;
      7.813106349722561   7.813107415959238  11.973288464744240  18.802480702520111]')

      order = [[-1], [1; 1], [-1; 1; 0], [-1; 1; 1; 0]]
      T = SparseRadialMap(4, order; λ = 0.0)

      optimize(T, X, nothing; start = 2)

      @test norm(T(X[:,1]) -   [6.831108232125667;
                                         0.07966241537360474;
                                        -0.18141931240603526;
                                         -2.0970508045942893])<1e-8
end


text:
# Defaults used throughout the package. Their role is to provide default values
# for various arguments and keyword arguments and a single point of control
# for the default behaviour of package functions.
#
# Note: Constants that have local scope (such as preprocessing.jl constants)
#       should not be added here.

const DEFAULT_LANGUAGE = Languages.English()
const DEFAULT_HASH_FUNCTION = hash
const DEFAULT_CARDINALITY = 100
const DEFAULT_FLOAT_TYPE = eltype(1.0)
const DEFAULT_DTM_TYPE = Int  # can be anything <:Real
const DEFAULT_NGRAM_COMPLEXITY = 1  # default ngram complexity
const DEFAULT_CORPUS_SPARSITY = 0.05  # if a term is present in less docs of a corpus than this percent, it is sparse
const DEFAULT_DOC_SPARSITY = 0.05  # if a term is present less than this percent in a document, it is sparse
const BM25_KAPPA = 2  # default value for BM25 κ parameter
const BM25_BETA = 0.75  # default value for BM25 β parameter

# Tokenization
# Regular expression on which to split text into tokens; It used by the `tokenize_stringanalysis` function
const DEFAULT_TOKENIZATION_REGEX = r"(,|\n|\r|\:|\\|\/|;|\.|\[|\]|\{|\}|\'|\`|\"|\"|\?|\!|\=|\~|\&|\s+)"
const DEFAULT_WT_TOKENIZER = WordTokenizers.nltk_word_tokenize
const DEFAULT_TOKENIZER = :default  # can be :default of :stringanalysis