text:
function setup_update_matrices()
    options = Dict(:verbose => false,
                   :eps_abs => 1e-08,
                   :eps_rel => 1e-08,
                   :polish => false,
                   :check_termination => 1)

    seed!(1)

    n = 5
    m = 8
    p = 0.7
    Pt = sprandn(n, n, p)
    Pt_new = copy(Pt)
    P = Pt * Pt' + sparse(I, n, n)
    #  PtI = findall(!iszero, P)
    #  (Pti, Ptj) = (getindex.(PtI, 1), getindex.(PtI, 2))
    Ptx = copy(Pt.nzval)

    Pt_newx = Ptx + 0.1 * randn(length(Ptx))
    # Pt_new = sparse(Pi, Pj, Pt_newx)
    P_new = Pt_new * Pt_new' + sparse(I, n, n)
    q = randn(n)
    A = sprandn(m, n, p)

    (Ai, Aj, _) = findnz(A)
    Ax = copy(A.nzval)
    A_newx = Ax + randn(length(Ax))
    A_new = sparse(Ai, Aj, A_newx)

    # A_new = copy(A)
    # A_new.nzval += randn(length(A_new.nzval))
    l = zeros(m)
    u = 30 .+ randn(m)

    problem = Dict()
    problem[:P] = P
    problem[:P_new] = P_new
    problem[:q] = q
    problem[:A] = A
    problem[:A_new] = A_new
    problem[:l] = l
    problem[:u] = u
    problem[:m] = m
    problem[:n] = n
    return problem, options
end

tol = 1e-5

# This unit test relies on a precomputed reference solution for a randomly generated problem (done in Julia v1.0).
# However, in newer versions of Julia the random number stream changed leading to a different problem being generated.
if VERSION < v"1.1"

        @testset "update_matrices" begin

            @testset "solve" begin
                problem, options = setup_update_matrices()

                (n, m, P, q, A, l, u) = (problem[:n], problem[:m], problem[:P], problem[:q], problem[:A], problem[:l], problem[:u])

                model = OSQP.Model()
                OSQP.setup!(model; P=P, q=q, A=A, l=l, u=u, options...)
                results = OSQP.solve!(model)

                # # Solve with Gurobi
                # using Gurobi
                # env = Gurobi.Env()
                # setparam!(env, "OutputFlag", 1)
                # model = Gurobi.Model(env, "model")
                # add_cvars!(model, q); update_model!(model)
                # add_constrs!(model, A, repmat(['<'], m), u); update_model!(model)
                # add_constrs!(model, A, repmat(['>'], m), l); update_model!(model)
                # add_qpterms!(model, P); update_model!(model)
                # optimize(model)
                # x_test = get_solution(model)
                # obj_test= get_objval(model)
                # y_test = -Gurobi.get_dblattrarray(model, "Pi", 1, 2 * m)
                # y_test = y_test[m+1:end] + y_test[1:m]
                        #
                # println("x_test = $(x_test)")
                # println("y_test = $(y_test)")
                # println("obj_test = $(obj_test)")

                x_test = [-0.324865, 0.598681, -0.066646, -0.00653471, -0.0736556]
                y_test = [-2.72542e-10, -1.927e-12, -0.547374, -2.25907e-10, -0.645086, -4.09874e-9, -1.4083e-11, -1.26234e-11]
                obj_test = -0.6736977821770016

                @test isapprox(results.x, x_test, atol=tol)
                @test isapprox(results.y, y_test, atol=tol)
                @test isapprox(results.info.obj_val, obj_test, atol=tol)
                @test results.info.status == :Solved

            end

                @testset "update_P" begin

                problem, options = setup_update_matrices()

                (n, m, P, q, A, l, u) = (problem[:n], problem[:m],
                             problem[:P], problem[:q], problem[:A], problem[:l], problem[:u])
                P_new = problem[:P_new]

                model = OSQP.Model()
                OSQP.setup!(model; P=P, q=q, A=A, l=l, u=u, options...)

                # Update matrix
                Pnew_triu = triu(P_new)
                Pnew_triu_idx = collect(1 : length(Pnew_triu.nzval))
                OSQP.update!(model, Px=Pnew_triu.nzval, Px_idx=Pnew_triu_idx)
                results = OSQP.solve!(model)

                # # Solve with Gurobi
                # using Gurobi
                # env = Gurobi.Env()
                # setparam!(env, "OutputFlag", 1)
                # model = Gurobi.Model(env, "model")
                # add_cvars!(model, q); update_model!(model)
                # add_constrs!(model, A, repmat(['<'], m), u); update_model!(model)
                # add_constrs!(model, A, repmat(['>'], m), l); update_model!(model)
                # add_qpterms!(model, P_new); update_model!(model)
                # optimize(model)
                # x_test = get_solution(model)
                # obj_test= get_objval(model)
                # y_test = -Gurobi.get_dblattrarray(model, "Pi", 1, 2 * m)
                # y_test = y_test[m+1:end] + y_test[1:m]
                        #
                # println("x_test = $(x_test)")
                # println("y_test = $(y_test)")
                # println("obj_test = $(obj_test)")

                x_test = [-0.324865, 0.598681, -0.066646, -0.00653471, -0.0736556]
                y_test = [-2.72542e-10, -1.927e-12, -0.547374, -2.25907e-10, -0.645086, -4.09874e-9, -1.4083e-11, -1.26234e-11]
                obj_test = -0.6736977821770016

                @test isapprox(results.x, x_test, atol=tol)
                @test isapprox(results.y, y_test, atol=tol)
                @test isapprox(results.info.obj_val, obj_test, atol=tol)
                @test results.info.status == :Solved
            end

                @testset "update_P_allind" begin

                problem, options = setup_update_matrices()

                (n, m, P, q, A, l, u) = (problem[:n], problem[:m],
                             problem[:P], problem[:q], problem[:A], problem[:l], problem[:u])
                P_new = problem[:P_new]

                model = OSQP.Model()
                OSQP.setup!(model; P=P, q=q, A=A, l=l, u=u, options...)

                # Update matrix
                Pnew_triu = triu(P_new)
                OSQP.update!(model, Px=Pnew_triu.nzval)
                results = OSQP.solve!(model)

                # # Solve with Gurobi
                # using Gurobi
                # env = Gurobi.Env()
                # setparam!(env, "OutputFlag", 1)
                # model = Gurobi.Model(env, "model")
                # add_cvars!(model, q); update_model!(model)
                # add_constrs!(model, A, repmat(['<'], m), u); update_model!(model)
                # add_constrs!(model, A, repmat(['>'], m), l); update_model!(model)
                # add_qpterms!(model, P_new); update_model!(model)
                # optimize(model)
                # x_test = get_solution(model)
                # obj_test= get_objval(model)
                # y_test = -Gurobi.get_dblattrarray(model, "Pi", 1, 2 * m)
                # y_test = y_test[m+1:end] + y_test[1:m]
                        #
                # println("x_test = $(x_test)")
                # println("y_test = $(y_test)")
                # println("obj_test = $(obj_test)")

                x_test = [-0.324865, 0.598681, -0.066646, -0.00653471, -0.0736556]
                y_test = [-2.72542e-10, -1.927e-12, -0.547374, -2.25907e-10, -0.645086, -4.09874e-9, -1.4083e-11, -1.26234e-11]
                obj_test = -0.6736977821770016

                @test isapprox(results.x, x_test, atol=tol)
                @test isapprox(results.y, y_test, atol=tol)
                @test isapprox(results.info.obj_val, obj_test, atol=tol)
                @test results.info.status == :Solved
            end

                @testset "update_A" begin

                problem, options = setup_update_matrices()

                (n, m, P, q, A, l, u) = (problem[:n], problem[:m],
                             problem[:P], problem[:q], problem[:A], problem[:l], problem[:u])
                A_new = problem[:A_new]

                model = OSQP.Model()
                OSQP.setup!(model; P=P, q=q, A=A, l=l, u=u, options...)

                # Update matrix
                A_new_idx = collect(1 : length(A_new.nzval))
                OSQP.update!(model, Ax=A_new.nzval, Ax_idx=A_new_idx)
                results = OSQP.solve!(model)

                # # Solve with Gurobi
                # using Gurobi
                # env = Gurobi.Env()
                # setparam!(env, "OutputFlag", 1)
                # model = Gurobi.Model(env, "model")
                # add_cvars!(model, q); update_model!(model)
                # add_constrs!(model, A_new, repmat(['<'], m), u); update_model!(model)
                # add_constrs!(model, A_new, repmat(['>'], m), l); update_model!(model)
                # add_qpterms!(model, P); update_model!(model)
                # optimize(model)
                # x_test = get_solution(model)
                # obj_test= get_objval(model)
                # y_test = -Gurobi.get_dblattrarray(model, "Pi", 1, 2 * m)
                # y_test = y_test[m+1:end] + y_test[1:m]
                        #
                # println("x_test = $(x_test)")
                # println("y_test = $(y_test)")
                # println("obj_test = $(obj_test)")
                        #
                x_test = [-0.0318085, 0.067069, -0.0242966, -0.0593736, -0.0321274]
                y_test = [-2.16989e-10, -1.01848, -0.516013, -0.0227263, -3.19721e-10, -1.04482, -3.4596e-10, -4.51608e-10]
                obj_test = -0.02187014865840703

                @test isapprox(results.x, x_test, atol=tol)
                @test isapprox(results.y, y_test, atol=tol)
                @test isapprox(results.info.obj_val, obj_test, atol=tol)
                @test results.info.status == :Solved
            end

            @testset "update_A_allind" begin

                problem, options = setup_update_matrices()

                (n, m, P, q, A, l, u) = (problem[:n], problem[:m],
                             problem[:P], problem[:q], problem[:A], problem[:l], problem[:u])
                A_new = problem[:A_new]

                model = OSQP.Model()
                OSQP.setup!(model; P=P, q=q, A=A, l=l, u=u, options...)

                # Update matrix
                OSQP.update!(model, Ax=A_new.nzval)
                results = OSQP.solve!(model)

                # # Solve with Gurobi
                # using Gurobi
                # env = Gurobi.Env()
                # setparam!(env, "OutputFlag", 1)
                # model = Gurobi.Model(env, "model")
                # add_cvars!(model, q); update_model!(model)
                # add_constrs!(model, A_new, repmat(['<'], m), u); update_model!(model)
                # add_constrs!(model, A_new, repmat(['>'], m), l); update_model!(model)
                # add_qpterms!(model, P); update_model!(model)
                # optimize(model)
                # x_test = get_solution(model)
                # obj_test= get_objval(model)
                # y_test = -Gurobi.get_dblattrarray(model, "Pi", 1, 2 * m)
                # y_test = y_test[m+1:end] + y_test[1:m]
                        #
                # println("x_test = $(x_test)")
                # println("y_test = $(y_test)")
                # println("obj_test = $(obj_test)")
                        #
                x_test = [-0.0318085, 0.067069, -0.0242966, -0.0593736, -0.0321274]
                y_test = [-2.16989e-10, -1.01848, -0.516013, -0.0227263, -3.19721e-10, -1.04482, -3.4596e-10, -4.51608e-10]
                obj_test = -0.02187014865840703

                @test isapprox(results.x, x_test, atol=tol)
                @test isapprox(results.y, y_test, atol=tol)
                @test isapprox(results.info.obj_val, obj_test, atol=tol)
                @test results.info.status == :Solved
            end

                @testset "update_P_A_indP_indA" begin

                problem, options = setup_update_matrices()

                (n, m, P, q, A, l, u) = (problem[:n], problem[:m],
                             problem[:P], problem[:q], problem[:A], problem[:l], problem[:u])
                P_new = problem[:P_new]
                A_new = problem[:A_new]

                model = OSQP.Model()
                OSQP.setup!(model; P=P, q=q, A=A, l=l, u=u, options...)

                # Update matrices P and A
                P_new_triu = triu(P_new)
                P_new_triu_idx = collect(1 : length(P_new_triu.nzval))
                A_new_idx = collect(1 : length(A_new.nzval))

                OSQP.update!(model, Px=P_new_triu.nzval, Px_idx=P_new_triu_idx,
                         Ax=A_new.nzval, Ax_idx=A_new_idx)
                results = OSQP.solve!(model)

                # # Solve with Gurobi
                # using Gurobi
                # env = Gurobi.Env()
                # setparam!(env, "OutputFlag", 1)
                # model = Gurobi.Model(env, "model")
                # add_cvars!(model, q); update_model!(model)
                # add_constrs!(model, A_new, repmat(['<'], m), u); update_model!(model)
                # add_constrs!(model, A_new, repmat(['>'], m), l); update_model!(model)
                # add_qpterms!(model, P_new); update_model!(model)
                # optimize(model)
                # x_test = get_solution(model)
                # obj_test= get_objval(model)
                # y_test = -Gurobi.get_dblattrarray(model, "Pi", 1, 2 * m)
                # y_test = y_test[m+1:end] + y_test[1:m]
                        #
                # println("x_test = $(x_test)")
                # println("y_test = $(y_test)")
                # println("obj_test = $(obj_test)")
                        #
                x_test = [-0.0318085, 0.067069, -0.0242966, -0.0593736, -0.0321274]
                y_test = [-2.16989e-10, -1.01848, -0.516013, -0.0227263, -3.19721e-10, -1.04482, -3.4596e-10, -4.51608e-10]
                obj_test = -0.02187014865840703

                @test isapprox(results.x, x_test, atol=tol)
                @test isapprox(results.y, y_test, atol=tol)
                @test isapprox(results.info.obj_val, obj_test, atol=tol)
                @test results.info.status == :Solved
            end

                @testset "update_P_A_indP" begin

                problem, options = setup_update_matrices()

                (n, m, P, q, A, l, u) = (problem[:n], problem[:m],
                             problem[:P], problem[:q], problem[:A], problem[:l], problem[:u])
                P_new = problem[:P_new]
                A_new = problem[:A_new]

                model = OSQP.Model()
                OSQP.setup!(model; P=P, q=q, A=A, l=l, u=u, options...)

                # Update matrices P and A
                P_new_triu = triu(P_new)
                P_new_triu_idx = collect(1 : length(P_new_triu.nzval))

                OSQP.update!(model, Px=P_new_triu.nzval, Px_idx=P_new_triu_idx,
                         Ax=A_new.nzval)
                results = OSQP.solve!(model)

                # # Solve with Gurobi
                # using Gurobi
                # env = Gurobi.Env()
                # setparam!(env, "OutputFlag", 1)
                # model = Gurobi.Model(env, "model")
                # add_cvars!(model, q); update_model!(model)
                # add_constrs!(model, A_new, repmat(['<'], m), u); update_model!(model)
                # add_constrs!(model, A_new, repmat(['>'], m), l); update_model!(model)
                # add_qpterms!(model, P_new); update_model!(model)
                # optimize(model)
                # x_test = get_solution(model)
                # obj_test= get_objval(model)
                # y_test = -Gurobi.get_dblattrarray(model, "Pi", 1, 2 * m)
                # y_test = y_test[m+1:end] + y_test[1:m]
                        #
                # println("x_test = $(x_test)")
                # println("y_test = $(y_test)")
                # println("obj_test = $(obj_test)")
                        #
                x_test = [-0.0318085, 0.067069, -0.0242966, -0.0593736, -0.0321274]
                y_test = [-2.16989e-10, -1.01848, -0.516013, -0.0227263, -3.19721e-10, -1.04482, -3.4596e-10, -4.51608e-10]
                obj_test = -0.02187014865840703

                @test isapprox(results.x, x_test, atol=tol)
                @test isapprox(results.y, y_test, atol=tol)
                @test isapprox(results.info.obj_val, obj_test, atol=tol)
                @test results.info.status == :Solved
            end

                @testset "update_P_A_indA" begin

                problem, options = setup_update_matrices()

                (n, m, P, q, A, l, u) = (problem[:n], problem[:m],
                             problem[:P], problem[:q], problem[:A], problem[:l], problem[:u])
                P_new = problem[:P_new]
                A_new = problem[:A_new]

                model = OSQP.Model()
                OSQP.setup!(model; P=P, q=q, A=A, l=l, u=u, options...)

                # Update matrices P and A
                P_new_triu = triu(P_new)
                A_new_idx = collect(1 : length(A_new.nzval))

                OSQP.update!(model, Px=P_new_triu.nzval,
                         Ax=A_new.nzval, Ax_idx=A_new_idx)
                results = OSQP.solve!(model)

                # # Solve with Gurobi
                # using Gurobi
                # env = Gurobi.Env()
                # setparam!(env, "OutputFlag", 1)
                # model = Gurobi.Model(env, "model")
                # add_cvars!(model, q); update_model!(model)
                # add_constrs!(model, A_new, repmat(['<'], m), u); update_model!(model)
                # add_constrs!(model, A_new, repmat(['>'], m), l); update_model!(model)
                # add_qpterms!(model, P_new); update_model!(model)
                # optimize(model)
                # x_test = get_solution(model)
                # obj_test= get_objval(model)
                # y_test = -Gurobi.get_dblattrarray(model, "Pi", 1, 2 * m)
                # y_test = y_test[m+1:end] + y_test[1:m]
                        #
                # println("x_test = $(x_test)")
                # println("y_test = $(y_test)")
                # println("obj_test = $(obj_test)")
                        #
                x_test = [-0.0318085, 0.067069, -0.0242966, -0.0593736, -0.0321274]
                y_test = [-2.16989e-10, -1.01848, -0.516013, -0.0227263, -3.19721e-10, -1.04482, -3.4596e-10, -4.51608e-10]
                obj_test = -0.02187014865840703

                @test isapprox(results.x, x_test, atol=tol)
                @test isapprox(results.y, y_test, atol=tol)
                @test isapprox(results.info.obj_val, obj_test, atol=tol)
                @test results.info.status == :Solved
            end

                @testset "update_P_A_allind" begin

                problem, options = setup_update_matrices()

                (n, m, P, q, A, l, u) = (problem[:n], problem[:m],
                             problem[:P], problem[:q], problem[:A], problem[:l], problem[:u])
                P_new = problem[:P_new]
                A_new = problem[:A_new]

                model = OSQP.Model()
                OSQP.setup!(model; P=P, q=q, A=A, l=l, u=u, options...)

                # Update matrices P and A
                P_new_triu = triu(P_new)

                OSQP.update!(model, Px=P_new_triu.nzval, Ax=A_new.nzval)
                results = OSQP.solve!(model)

                # # Solve with Gurobi
                # using Gurobi
                # env = Gurobi.Env()
                # setparam!(env, "OutputFlag", 1)
                # model = Gurobi.Model(env, "model")
                # add_cvars!(model, q); update_model!(model)
                # add_constrs!(model, A_new, repmat(['<'], m), u); update_model!(model)
                # add_constrs!(model, A_new, repmat(['>'], m), l); update_model!(model)
                # add_qpterms!(model, P_new); update_model!(model)
                # optimize(model)
                # x_test = get_solution(model)
                # obj_test= get_objval(model)
                # y_test = -Gurobi.get_dblattrarray(model, "Pi", 1, 2 * m)
                # y_test = y_test[m+1:end] + y_test[1:m]
                        #
                # println("x_test = $(x_test)")
                # println("y_test = $(y_test)")
                # println("obj_test = $(obj_test)")
                        #
                x_test = [-0.0318085, 0.067069, -0.0242966, -0.0593736, -0.0321274]
                y_test = [-2.16989e-10, -1.01848, -0.516013, -0.0227263, -3.19721e-10, -1.04482, -3.4596e-10, -4.51608e-10]
                obj_test = -0.02187014865840703

                @test isapprox(results.x, x_test, atol=tol)
                @test isapprox(results.y, y_test, atol=tol)
                @test isapprox(results.info.obj_val, obj_test, atol=tol)
                @test results.info.status == :Solved
            end

        end
end


text:
include("correct-guess.jl")
include("guess-vars.jl")
include("other.jl")
include("primal-project.jl")
include("gertz_init.jl")
include("lp.jl")

function linear_cons(nlp::Class_CUTEst, x0::Vector)
    a1 = eval_a(nlp, x0)
    d = randn(length(x0))
    a2 = eval_a(nlp, x0 + d)
    a3 = eval_a(nlp, x0 + 2 * d)

    return abs(a1 - 2 * a2 + a3) .< 1e-6
end

function mehrotra_init(nlp::Class_CUTEst, pars::Class_parameters, timer::class_advanced_timer)
    start_advanced_timer(timer, "INIT/x")
    x0 = suggested_starting_point(nlp)

    if pars.init.start_satisfying_bounds
      x = projection_onto_bounds_ipopt_style( nlp, pars, x0 )
    end
    pause_advanced_timer(timer, "INIT/x")

    a, J, g = eval_init(nlp, pars, timer, x)

    s, y = mehortra_guess( nlp, pars, timer, x, a, J, g )

    @assert(all(s .>= 0.0))

    li = linear_cons(nlp)
    nl = li .!= true
    ineq = ineq_cons(nlp)
    eq = ineq .!= true
    nl_eq = ineq .& nl
    nl_ineq = eq .& nl

    if pars.init.mehotra_scaling
      mu = mean(s .* y)
      conWeight = ((s - a) / mu)
      #v = conWeight[conWeight .> 0.0]
      #conWeight[conWeight .> 0.0] = min.(10.0,max.(v,1e-4))
    else
      mu = (1e-6 + norm(s,Inf) + norm(g,Inf))
      #mu = 1e-6 + norm(g,1) / length(s)
      conWeight = zeros(length(s))
      ais = cons_indicies(nlp)
      conWeight[ais] = 1.0

      #li = linear_cons(nlp, x)
    end

    #@show li
    #@show nl

    #@show li
    if true
      conWeight[nl_eq] *= pars.init.nl_eq_scale
      conWeight[nl_ineq] *= pars.init.nl_ineq_scale
      #@show li, nl_ineq
      conWeight[li] *= pars.init.linear_scale
      #mu *= pars.init.mu_scale
    end
    #*= pars.init.linear_scale

    @assert(all(conWeight .>= 0.0))

    iter_init = correct_guess3( nlp, pars, timer, x, a, J, g, y, mu, conWeight )
    #return correct_guess1( nlp, pars, timer, x0, a, J, g, s, y)
    #
    #@show li

    #@show li
    iter_init.frac_bd_predict[li] = pars.ls.fraction_to_boundary_linear
    iter_init.frac_bd[li] = pars.ls.fraction_to_boundary_linear

    @assert(all(iter_init.point.s .>= 0.0))
    @assert is_feasible(iter_init,pars.ls.comp_feas)

    return iter_init
end

function simple_init()

end

ext: jl


text:
println('x')
# x
println(Char(120))
# x
println(Char(0x110000))
println(isvalid(Char(0x110000)))
# false
println(Int('\0'))
# 0
println(Int('\t'))
# 9
println(Int('\n'))
# 10

v = [1, 2, 3]
println("v: $v")
# v: [1, 2, 3]

println(search("xylophone", 'o'))
# 4
println(search("xylophone", 'o', 5))
# 7                              ^ from
println(search("xylophone", 'o', 8))
# 0                              ^ from
println(repeat(".:Z:.", 10))
# .:Z:..:Z:..:Z:..:Z:..:Z:..:Z:..:Z:..:Z:..:Z:..:Z:.

max_stars_repo_name: dannypsnl/languages-learn

max_issues_repo_name: dannypsnl/languages-learn

max_forks_repo_name: dannypsnl/languages-learn


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.4'
#       jupytext_version: 1.2.4
#   kernelspec:
#     display_name: Julia 1.3.1
#     language: julia
#     name: julia-1.3
# ---

# # Example 3
#
# This example first reads in Global Ocean variables, generated by the [MITgcm](https://mitgcm.readthedocs.io/en/latest/index.html), which are split into a collection of subdomain arrays ( _tiles_ ) using `MeshArrays.jl`. It then writes them to a collection of `NetCDF` files ( _nctiles_ ) using `NCTiles.jl`.

# ### Packages & Helper Functions
#
# _These will be used throughout the notebook_

# +
if false
    using Pkg
    Pkg.add(PackageSpec(name="MeshArrays", rev="master"))
    Pkg.add(PackageSpec(name="MITgcmTools", rev="master"))
end

using NCTiles, MITgcmTools
include("helper_functions.jl");
# -

# ### File Paths & I/O Back-End
#
# _These will be used throughout the notebook_

# +
# File Paths
inputs = "input/"
get_testcases_if_needed(inputs)

outputs = "output/ex3/"
if ~ispath(outputs); mkpath(outputs); end

# I/O Back-End
nc=NCTiles.NCDatasets
# -

# ### Process Global Ocean Variables
#
# Here we process a two-dimensional field (`ETAN`), two three-dimensional tracer fields (`THETA`, `SALT`), and the three components of a vector field (`UVELMASS`, `VVELMASS` and `WVELMASS`). In each case, `flds` denotes the lazy representation of the processing chain (incl. all needed metadata) which the `write` function instantiates (i.e. outputs to files).
#
# _Note: on a `C-grid` these components are staggered in space._

# +
writedir=outputs
readme = readlines(inputs)

#2D example
(flds,savename,readme)=prep_nctiles_native(inputs,"state_2d_set1","ETAN",Float32)
write(flds,savename,README=readme);

# +
#3D examples
(flds,savename,readme)=prep_nctiles_native(inputs,"state_3d_set1","THETA",Float32);
write(flds,savename,README=readme);

(flds,savename,readme)=prep_nctiles_native(inputs,"state_3d_set1","SALT",Float32);
write(flds,savename,README=readme);

# +
#3D vector field example
(flds,savename,readme)=prep_nctiles_native(inputs,"trsp_3d_set1","UVELMASS",Float32);
write(flds,savename,README=readme);

(flds,savename,readme)=prep_nctiles_native(inputs,"trsp_3d_set1","VVELMASS",Float32);
write(flds,savename,README=readme);

(flds,savename,readme)=prep_nctiles_native(inputs,"trsp_3d_set1","WVELMASS",Float32);
write(flds,savename,README=readme);
# -


text:
"""
An [`Activity`](@ref)'s type. Available values are `GAME`, `STREAMING`,
`LISTENING`, `WATCHING`, and `COMPETING`.
More details [here](https://discordapp.com/developers/docs/topics/gateway#activity-object-activity-types).
"""
module ActivityType
    const GAME=0
    const STREAMING=1
    const LISTENING=2
    const WATCHING=3
    const CUSTOM=4
    const COMPETING=5
end
"""
Flags which indicate what an [`Activity`](@ref) payload contains.
More details [here](https://discordapp.com/developers/docs/topics/gateway#activity-object-activity-flags).
"""
module ActivityFlags
    const INSTANCE=1<<0
    const JOIN=1<<1
    const SPECTATE=1<<2
    const JOIN_REQUEST=1<<3
    const SYNC=1<<4
    const PLAY=1<<5
end
"""
[`AuditLog`](@ref) action types.
More details [here](https://discordapp.com/developers/docs/resources/audit-log#audit-log-entry-object-audit-log-events).
"""
module ActionType
    const GUILD_UPDATE=1
    const CHANNEL_CREATE=10
    const CHANNEL_UPDATE=11
    const CHANNEL_DELETE=12
    const CHANNEL_OVERWRITE_CREATE=13
    const CHANNEL_OVERWRITE_UPDATE=14
    const CHANNEL_OVERWRITE_DELETE=15
    const MEMBER_KICK=20
    const MEMBER_PRUNE=21
    const MEMBER_BAN_ADD=22
    const MEMBER_BAN_REMOVE=23
    const MEMBER_UPDATE=24
    const MEMBER_ROLE_UPDATE=25
    const ROLE_CREATE=30
    const ROLE_UPDATE=31
    const ROLE_DELETE=32
    const INVITE_CREATE=40
    const INVITE_UPDATE=41
    const INVITE_DELETE=42
    const WEBHOOK_CREATE=50
    const WEBHOOK_UPDATE=51
    const WEBHOOK_DELETE=52
    const EMOJI_CREATE=60
    const EMOJI_UPDATE=61
    const EMOJI_DELETE=62
    const MESSAGE_DELETE=72
end
"""
A [`DiscordChannel`](@ref)'s type. See full list
at https://discord.com/developers/docs/resources/channel#channel-object-channel-types
"""
module ChannelTypes
    const GUILD_TEXT = 0
    const DM = 1
    const GUILD_VOICE = 2
    const GROUP_DM = 3
    const GUILD_CATEGORY = 4
    const GUILD_NEWS = 5
    const GUILD_STORE = 6
    const GUILD_NEWS_THREAD = 10
    const GUILD_PUBLIC_THREAD = 11
    const GUILD_PRIVATE_THREAD = 12
    const GUILD_STAGE_VOICE = 13
end
"""
A [`Guild`](@ref)'s verification level.
More details [here](https://discordapp.com/developers/docs/resources/guild#guild-object-verification-level).
"""
module VerificationLevel
    const NONE = 0
    const LOW = 1
    const MEDIUM = 2
    const HIGH = 3
    const VERY_HIGH = 4
end

"""
A [`Guild`](@ref)'s default message notification level.
More details [here](https://discordapp.com/developers/docs/resources/guild#guild-object-default-message-notification-level).
"""
module MessageNotificationLevel
    const ALL_MESSAGES = 0
    const ONLY_MENTIONS = 1
end

"""
A [`Guild`](@ref)'s explicit content filter level.
More details [here](https://discordapp.com/developers/docs/resources/guild#guild-object-explicit-content-filter-level).
"""
module ExplicitContentFilterLevel
    const DISABLED = 0
    const MEMBERS_WITHOUT_ROLES = 1
    const ALL_MEMBERS = 2
end
"""
A [`Guild`](@ref)'s MFA level.
More details [here](https://discordapp.com/developers/docs/resources/guild#guild-object-mfa-level).
"""
module MFALevel
    const NONE = 0
    const ELEVATED = 1
end

module InteractionType
    const PING = 1
    const APPLICATIONCOMMAND = 2
    const MESSAGECOMPONENT = 3
end

module ApplicationCommandType
    const CHATINPUT = 1
    const UI = 2
    const MESSAGE = 3
end

module ComponentType
    const ACTIONROW = 1
    const BUTTON = 2
    const SELECTMENU = 3
end

module OptionType
    const SUB_COMMAND = 1
    const SUB_COMMAND_GROUP = 2
    const STRING = 3
    const INTEGER = 4
    const BOOLEAN = 5
    const USER = 6
    const CHANNEL = 7
    const ROLE = 8
    const MENTIONABLE = 9
    const NUMBER = 10
end

"""
A [`Message`](@ref)'s type.
More details [here](https://discordapp.com/developers/docs/resources/channel#message-object-message-types).
"""
module MessageType
    const DEFAULT = 0
    const RECIPIENT_ADD = 1
    const RECIPIENT_REMOVE = 2
    const CALL = 3
    const CHANNEL_NAME_CHANGE = 4
    const CHANNEL_ICON_CHANGE = 5
    const CHANNEL_PINNED_MESSAGE = 6
    const GUILD_MEMBER_JOIN = 7
    const USER_PREMIUM_GUILD_SUBSCRIPTION = 8
    const USER_PREMIUM_GUILD_SUBSCRIPTION_TIER_1 = 9
    const USER_PREMIUM_GUILD_SUBSCRIPTION_TIER_2 = 10
    const USER_PREMIUM_GUILD_SUBSCRIPTION_TIER_3 = 11
    const CHANNEL_FOLLOW_ADD = 12
    const GUILD_DISCOVERY_DISQUALIFIED = 14
    const GUILD_DISCOVERY_REQUALIFIED = 15
    const GUILD_DISCOVERY_GRACE_PERIOD_INITIAL_WARNING = 16
    const GUILD_DISCOVERY_GRACE_PERIOD_FINAL_WARNING = 17
    const THREAD_CREATED = 18
    const REPLY = 19
    const CHAT_INPUT_COMMAND = 20
    const THREAD_STARTER_MESSAGE = 21
    const GUILD_INVITE_REMINDER = 22
    const CONTEXT_MENU_COMMAND = 23
end

"""
A [`Message`](@ref)'s activity type.
More details [here](https://discordapp.com/developers/docs/resources/channel#message-object-message-activity-types).
"""
module MessageActivityType
    const JOIN = 1
    const SPECTATE = 2
    const LISTEN = 3
    const JOIN_REQUEST = 5
end


text:
using Replay

repl_script = """
2+2
print("")
display([1])
display([1 2; 3 4])
@time 1+1
using Pkg; Pkg.activate(".")
] add Example
rm Example
st
$CTRL_C
"""

replay(repl_script, stdout, color = :yes, julia_project=@__DIR__, use_ghostwriter=true)


text:
module admm_test

using Base.Test, ConstrainedLasso, GLMNet

info("Test lsq_constrsparsereg_admm: sum-to-zero constraint (single param value)")
# set up
srand(123)
n, p = 100, 20
# truth with sum constraint sum(β) = 0
β = zeros(p)
β[1:round(Int, p / 4)] = 0
β[(round(Int, p / 4) + 1):round(Int, p / 2)] = 1
β[(round(Int, p / 2) + 1):round(Int, 3p / 4)] = 0
β[(round(Int, 3p / 4) + 1):p] = -1
# generate data
X = randn(n, p)
y = X * β + randn(n)
# fit at a fixed parameter value
ρ = 2.0
β̂admm1 = lsq_constrsparsereg_admm(X, y, ρ; proj = x -> x - mean(x))
@test sum(β̂admm1)≈0.0 atol=1e-5

info("Test lsq_constrsparsereg_admm: sum-to-zero constraint (multiple param values)")
# fit at fixed parameter values
ρ = 1.0:2.0:20.0
β̂admm2 = lsq_constrsparsereg_admm(X, y, ρ; proj = x -> x - mean(x))
@testset "zero-sum for multiple param values" begin
for si in sum(β̂admm2, 1)
    @test si≈0.0 atol=1e-5
end
end

info("Test lsq_constrsparsereg_admm: non-negativity constraint")
# set up
n, p = 100, 20
β = zeros(p)
# truth
β[1:10] = 1:10
# generate data
srand(41)
X = randn(n, p)
y = X * β + randn(n)
# fit at the fixed parameter value; admmvaryscale=true
ρ = 3
β̂admm3 = lsq_constrsparsereg_admm(X, y, ρ; proj = x -> clamp.(x, 0, Inf),
        admmvaryscale = true)
@test all(β̂admm3 .>= 0)

end # end of module


text:
using NWNTools
using Test
using DelimitedFiles, StaticArrays

@testset "NWNTools.jl" begin
    @testset "Benchmark1" begin
        data = readdlm("benchmark1.txt",',', BigFloat)
        dims = SA[BigFloat("3.0"),BigFloat("5.0")]
        elecs = [Line(SVector{2}(data[1,:]),SVector{2}(data[2,:]))]
        grnds = [Line(SVector{2}(data[3,:]),SVector{2}(data[4,:]))]
        Rⱼ = BigFloat("20.0")
        V₊ = BigFloat("5.0")
        volts = [V₊]
        ρₑ = BigFloat("22.63676")/BigFloat("1000")
        D = BigFloat("0.06")

        lines = Line{BigFloat,2}[]
        for i ∈ 5:2:size(data)[1]
            push!(lines, Line(SVector{2}(data[i,:]),SVector{2}(data[i+1,:])))
        end
        props = repeat([WireProp(ρₑ,D)],length(lines))

        nwn_JDA = NWN_circuit{:JDA}(lines, props, elecs, volts, grnds, dims, Rⱼ,Rⱼ)
        R_JDA = sheet_resistance(nwn_JDA)
        nwn_MNR = NWN_circuit{:MNR}(lines, props, elecs, volts, grnds, dims, Rⱼ,Rⱼ)
        R_MNR = sheet_resistance(nwn_MNR)
        @test R_JDA≈(8*Rⱼ/3)
        numerator = 72*(3+5*sqrt(big(2)))*ρₑ^2+2*D^2*big(π)*Rⱼ*((81+69*sqrt(big(2)))*ρₑ+20*D^2*big(π)*Rⱼ)
        denominator = 15*D^2*big(π)*(2*(1+sqrt(big(2)))*ρₑ+D^2*big(π)*Rⱼ)
        R_MNR_analytic = numerator/denominator
        @test R_MNR≈R_MNR_analytic
end


text:
# using ErgodicControl
# using Test
#
# @testset "ErgodicControl.jl" begin
#     # Write your tests here.
# end

using ErgodicControl
using Base.Test

# write your own tests here
@test 1 == 1

em = ErgodicManagerR2("single gaussian", K=5, bins=100)
x0 = [0.4,0.1]
N = 40
h = 0.1

tm = TrajectoryManager(x0, h, N, ConstantInitializer([0.0,0.0]))
xd, ud = pto_trajectory(em, tm)

include("dubins.jl")

include("r2t.jl")

include("multi.jl")

include("r3.jl")

include("smc.jl")

include("se2.jl")

include("linear.jl")

include("max.jl")


text:
export SymmetricIntervalHull, ⊡,
       an_element

"""
    SymmetricIntervalHull{N, S<:LazySet{N}} <: AbstractHyperrectangle{N}

Type that represents the symmetric interval hull of a compact set.

### Fields

- `X`     -- compact set
- `cache` -- partial storage of already computed bounds, organized as mapping
    from dimension to tuples `(bound, valid)`, where `valid` is a flag
    indicating if the `bound` entry has been computed

### Notes

The symmetric interval hull can be computed with ``2n`` support vector queries
of unit vectors, where ``n`` is the dimension of the wrapped set (i.e., two
queries per dimension).
When asking for the support vector for a direction ``d``, one needs ``2k`` such
queries, where ``k`` is the number of non-zero entries in ``d``.

However, if one asks for many support vectors in a loop, the number of
computations may exceed ``2n``.
To be most efficient in such cases, this type stores the intermediately computed
bounds in the `cache` field.

The set `X` must be bounded. The flag `check_boundedness` (which defaults to `true`)
can be used to elide the boundedness check in the inner constructor. Misuse of this
flag can result in incorrect behavior.

The symmetric inteval hull of a set is a hyperrectangle, hence in particular convex.

An alias for this function is `⊡`.
"""
struct SymmetricIntervalHull{N, S<:LazySet{N}} <: AbstractHyperrectangle{N}
    X::S
    cache::Vector{N}

    # default constructor that initializes cache
    function SymmetricIntervalHull(X::S; check_boundedness::Bool=true) where {N, S<:LazySet{N}}
        @assert !check_boundedness || isbounded(X) "the symmetric interval hull is only defined for bounded sets"
        cache = fill(-one(N), dim(X))
        return new{N, S}(X, cache)
    end
end

"""
    ⊡

Alias for `SymmetricIntervalHull`.
"""
const ⊡ = SymmetricIntervalHull

isoperationtype(::Type{<:SymmetricIntervalHull}) = true
isconvextype(::Type{<:SymmetricIntervalHull}) = true

"""
    SymmetricIntervalHull(∅::EmptySet)

The symmetric interval hull of an empty set.

### Input

- `∅` -- an empty set

### Output

The empty set because it is absorbing for the symmetric interval hull.
"""
SymmetricIntervalHull(∅::EmptySet) = ∅

# --- AbstractHyperrectangle interface functions ---

"""
    radius_hyperrectangle(sih::SymmetricIntervalHull, i::Int)

Return the box radius of the symmetric interval hull of a set in a given
dimension.

### Input

- `sih` -- symmetric interval hull of a set
- `i`   -- dimension of interest

### Output

The radius in the given dimension.
If it was computed before, this is just a look-up, otherwise it requires two
support vector computations.
"""
function radius_hyperrectangle(sih::SymmetricIntervalHull, i::Int)
    return get_radius!(sih, i)
end

"""
    radius_hyperrectangle(sih::SymmetricIntervalHull)

Return the box radius of the symmetric interval hull of a set in every
dimension.

### Input

- `sih` -- symmetric interval hull of a set

### Output

The box radius of the symmetric interval hull of a set.

### Notes

This function computes the symmetric interval hull explicitly.
"""
function radius_hyperrectangle(sih::SymmetricIntervalHull)
    n = dim(sih)
    for i in 1:n
        get_radius!(sih, i, n)
    end
    return sih.cache
end

"""
    center(sih::SymmetricIntervalHull{N}, i::Int) where {N}

Return the center along a given dimension of the symmetric interval hull of a
set.

### Input

- `sih` -- symmetric interval hull of a set
- `i`   -- dimension of interest

### Output

The center along a given dimension of the symmetric interval hull of a set.
"""
@inline function center(sih::SymmetricIntervalHull{N}, i::Int) where {N}
    @boundscheck _check_bounds(sih, i)
    return zero(N)
end

# --- AbstractCentrallySymmetric interface functions ---

"""
    center(sih::SymmetricIntervalHull{N}) where {N}

Return the center of the symmetric interval hull of a set.

### Input

- `sih` -- symmetric interval hull of a set

### Output

The origin.
"""
function center(sih::SymmetricIntervalHull{N}) where {N}
    return zeros(N, dim(sih))
end

# --- LazySet interface functions ---

"""
    dim(sih::SymmetricIntervalHull)

Return the dimension of the symmetric interval hull of a set.

### Input

- `sih` -- symmetric interval hull of a set

### Output

The ambient dimension of the symmetric interval hull of a set.
"""
function dim(sih::SymmetricIntervalHull)
    return dim(sih.X)
end

"""
    σ(d::AbstractVector, sih::SymmetricIntervalHull)

Return the support vector of the symmetric interval hull of a set in a given
direction.

### Input

- `d`   -- direction
- `sih` -- symmetric interval hull of a set

### Output

The support vector of the symmetric interval hull of a set in the given
direction.
If the direction has norm zero, the origin is returned.

### Algorithm

For each non-zero entry in `d` we need to either look up the bound (if it has
been computed before) or compute it, in which case we store it for future
queries.
One such computation just asks for the support vector of the underlying set for
both the positive and negative unit vector in the respective dimension.
"""
function σ(d::AbstractVector, sih::SymmetricIntervalHull)
    n = length(d)
    @assert n == dim(sih) "cannot compute the support vector of a " *
        "$(dim(sih))-dimensional set along a vector of length $n"

    N = promote_type(eltype(d), eltype(sih))
    svec = similar(d)
    for i in eachindex(d)
        if d[i] == zero(N)
            svec[i] = zero(N)
        else
            svec[i] = sign(d[i]) * get_radius!(sih, i, n)
        end
    end
    return svec
end

# Faster support vector calculation for sev and SymmetricIntervalHull
function σ(d::SingleEntryVector, sih::SymmetricIntervalHull)
    N = promote_type(eltype(d), eltype(sih))
    @assert length(d) == dim(sih) "cannot compute the support vector of a " *
                                  "$(dim(sih))-dimensional set along a vector of length $(d.n)"

    entry = radius_hyperrectangle(sih, d.i)
    if d.v < zero(N)
        entry = -entry
    end
    return SingleEntryVector(d.i, d.n, entry)
end

# Faster support function calculation for sev and SymmetricIntervalHull
function ρ(d::SingleEntryVector, H::SymmetricIntervalHull)
    @assert length(d) == dim(H) "a $(d.n)-dimensional vector is " *
                                "incompatible with a $(dim(H))-dimensional set"
    return abs(d.v) * radius_hyperrectangle(H, d.i)
end

# --- SymmetricIntervalHull functions ---

"""
    get_radius!(sih::SymmetricIntervalHull{N},
                i::Int,
                n::Int=dim(sih)) where {N}

Compute the radius of the symmetric interval hull of a set in a given dimension.

### Input

- `sih` -- symmetric interval hull of a set
- `i`   -- dimension in which the radius should be computed
- `n`   -- (optional, default: `dim(sih)`) set dimension

### Output

The radius of the symmetric interval hull of a set in a given dimension.

### Algorithm

We ask for the support vector of the underlying set for both the positive and
negative unit vector in the dimension `i`.
"""
function get_radius!(sih::SymmetricIntervalHull{N},
                     i::Int,
                     n::Int=dim(sih)) where {N}
    if sih.cache[i] == -one(N)
        d = SingleEntryVector(i, n, one(N))
        right_bound = σ(d, sih.X)
        d = SingleEntryVector(i, n, -one(N))
        left_bound = σ(d, sih.X)
        sih.cache[i] = max(right_bound[i], abs(left_bound[i]))
    end
    return sih.cache[i]
end


text:
@testset "DataCollection" begin
    mutable struct Nested
        data::Vector{Float64}
    end

    function initialize()
        model = ABM(
            Agent3,
            GridSpace((10, 10));
            properties = Dict(
                :year => 0,
                :tick => 0,
                :flag => false,
                :container => Float64[],
                :deep => Nested([20.0, 52.1])
            ),
        )
        add_agent!((4, 3), model, 0.1)
        add_agent!((7, 5), model, 0.35)
        add_agent!((2, 9), model, 0.67)
        return model
    end

    function agent_step!(agent, model)
        agent.weight += 0.05
        if model.tick % 365 == 0
            agent.weight *= 2
        end
    end
    function model_step!(model)
        model.tick += 1
        model.flag = !model.flag
        if model.tick % 365 == 0
            model.year += 1
            model.deep.data[1] += 0.5
        end
    end

    x_position(agent) = first(agent.pos)
    model = initialize()

    @testset "DataFrame init" begin
        @test init_agent_dataframe(model, nothing) == DataFrame()
        @test collect_agent_data!(DataFrame(), model, nothing, 1) == DataFrame()

        @test init_model_dataframe(model, nothing) == DataFrame()
        @test collect_model_data!(DataFrame(), model, nothing, 1) == DataFrame()

        props = [:weight]
        @test sprint(
            show,
            "text/csv",
            describe(init_agent_dataframe(model, props), :eltype),
        ) ==
              "\"variable\",\"eltype\"\n\"step\",\"Int64\"\n\"id\",\"Int64\"\n\"weight\",\"Float64\"\n"
        props = [:year]
        @test sprint(
            show,
            "text/csv",
            describe(init_model_dataframe(model, props), :eltype),
        ) == "\"variable\",\"eltype\"\n\"step\",\"Int64\"\n\"year\",\"Int64\"\n"

        @test_throws ErrorException init_agent_dataframe(model, [:UNKNOWN])
    end

    @testset "aggname" begin
        @test aggname(:weight) == "weight"
        @test aggname(:weight, mean) == "mean_weight"
        @test aggname(x_position, length) == "length_x_position"
        @test aggname((x_position, length)) == "length_x_position"
        ypos(a) = a.pos[2] > 5
        @test aggname((x_position, length, ypos)) == "length_x_position_ypos"
    end

    @testset "Aggregate Collections" begin
        props = [:weight]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props, 1)
        # Expecting weight values of all three agents. ID and step included.
        @test size(df) == (3, 3)
        @test propertynames(df) == [:step, :id, :weight]
        @test mean(df[!, :weight]) ≈ 0.37333333333

        props = [(:weight, mean)]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props, 1)
        # Activate aggregation. Weight column is expected to be one value for this step,
        # renamed mean(weight). ID is meaningless and will therefore be dropped.
        @test size(df) == (1, 2)
        @test propertynames(df) == [:step, :mean_weight]
        @test df[1, aggname(:weight, mean)] ≈ 0.37333333333

        # Add a function as a property
        props = [:weight, x_position]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props, 1)
        @test size(df) == (3, 4)
        @test propertynames(df) == [:step, :id, :weight, :x_position]
        @test mean(df[!, :x_position]) ≈ 4.3333333

        props = [(:weight, mean), (x_position, mean)]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props, 1)
        @test size(df) == (1, 3)
        @test propertynames(df) == [:step, :mean_weight, :mean_x_position]
        @test df[1, aggname(x_position, mean)] ≈ 4.3333333

        xtest(agent) = agent.pos[1] > 5
        ytest(agent) = agent.pos[2] > 5

        props = [(:weight, mean, ytest)]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props, 1)
        @test size(df) == (1, 2)
        @test propertynames(df) == [:step, :mean_weight_ytest]
        @test df[1, aggname((:weight, mean, ytest))] ≈ 0.67

        props = [(:weight, mean), (:weight, mean, ytest)]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props, 1)
        @test size(df) == (1, 3)
        @test propertynames(df) == [:step, :mean_weight, :mean_weight_ytest]
        @test df[1, aggname(:weight, mean)] ≈ 0.37333333333
        @test df[1, aggname(:weight, mean, ytest)] ≈ 0.67

        props = [(:weight, mean, xtest), (:weight, mean, ytest)]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props, 1)
        @test size(df) == (1, 3)
        @test propertynames(df) == [:step, :mean_weight_xtest, :mean_weight_ytest]
        @test df[1, aggname(:weight, mean, xtest)] ≈ 0.35
        @test df[1, aggname(:weight, mean, ytest)] ≈ 0.67
    end

    @testset "High-level API for Collections" begin
        # Extract data from the model every year for five years,
        # with the average `weight` of all agents every six months.
        each_year(model, step) = step % 365 == 0
        six_months(model, step) = step % 182 == 0
        agent_data, model_data = run!(
            model,
            agent_step!,
            model_step!,
            365 * 5;
            when_model = each_year,
            when = six_months,
            mdata = [:flag, :year],
            adata = [(:weight, mean)],
        )

        @test size(agent_data) == (11, 2)
        @test propertynames(agent_data) == [:step, :mean_weight]
        @test maximum(agent_data[!, :step]) == 1820

        @test size(model_data) == (6, 3)
        @test propertynames(model_data) == [:step, :flag, :year]
        @test maximum(model_data[!, :step]) == 1825

        agent_data, model_data = run!(
            model,
            agent_step!,
            model_step!,
            365 * 5;
            when_model = [1, 365 * 5],
            when = false,
            mdata = [:flag, :year],
            adata = [(:weight, mean)],
        )
        @test size(agent_data) == (0, 2)
        @test size(model_data) == (2, 3)

        _, model_data = run!(
            model,
            agent_step!,
            model_step!,
            365 * 5;
            when_model = [1, 365 * 5],
            when = false,
            mdata = [:deep],
            obtainer = deepcopy,
        )
        @test model_data[1, :deep].data[1] < model_data[end, :deep].data[1]
    end

    @testset "Low-level API for Collections" begin
        # Generate three separate dataframes using the low level API.
        # User controls the evolution, as well as the identity of each `step`.
        # Note in this example daily data are labelled with a daily `step`,
        # and yearly data with a yearly `step`.
        model = initialize()
        model_props = [:flag, :year]
        agent_agg = [(:weight, mean)]
        agent_props = [:weight]
        daily_model_data = init_model_dataframe(model, model_props)
        daily_agent_aggregate = init_agent_dataframe(model, agent_agg)
        yearly_agent_data = init_agent_dataframe(model, agent_props)

        for year in 1:5
            for day in 1:365
                step!(model, agent_step!, model_step!, 1)
                collect_model_data!(daily_model_data, model, model_props, day * year)
                collect_agent_data!(daily_agent_aggregate, model, agent_agg, day * year)
            end
            collect_agent_data!(yearly_agent_data, model, agent_props, year)
        end

        @test size(daily_model_data) == (1825, 3)
        @test propertynames(daily_model_data) == [:step, :flag, :year]
        @test maximum(daily_model_data[!, :step]) == 1825

        @test size(daily_agent_aggregate) == (1825, 2)
        @test propertynames(daily_agent_aggregate) == [:step, :mean_weight]
        @test maximum(daily_agent_aggregate[!, :step]) == 1825

        @test size(yearly_agent_data) == (15, 3)
        @test propertynames(yearly_agent_data) == [:step, :id, :weight]
        @test maximum(yearly_agent_data[!, :step]) == 5

        @test dummystep(model) == nothing
        @test dummystep(model[1], model) == nothing
        tick = model.tick
        step!(model, agent_step!, 1)
        @test tick == model.tick
        stop(m, s) = m.year == 6
        step!(model, agent_step!, model_step!, stop)
        @test model.tick == 365 * 6
        step!(model, agent_step!, model_step!, 1, false)
    end

    @testset "Mixed-ABM collections" begin
        model = ABM(Union{Agent3, Agent4}, GridSpace((10,10)); warn = false)
        add_agent_pos!(Agent3(1, (6,8), 54.65),model)
        add_agent_pos!(Agent4(2, (10,7), 5),model)

        # Expect position type (both agents have it)
        props = [:pos]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props)
        @test size(df) == (2, 4)
        @test propertynames(df) == [:step, :id, :agent_type, :pos]
        @test df[1, :pos] == model[1].pos
        @test df[2, :pos] == model[2].pos

        # Expect weight for Agent3, but missing for Agent4
        props = [:weight]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props)
        @test size(df) == (2, 4)
        @test df[1, :weight] == model[1].weight
        @test ismissing(df[2, :weight])

        # Expect similar output, but using functions
        wpos(a) = a.pos[1] + a.weight
        props = [wpos]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props)
        @test size(df) == (2, 4)
        @test df[1, :wpos] == model[1].pos[1] + model[1].weight
        @test ismissing(df[2, :wpos])

        add_agent_pos!(Agent3(3, (2,4), 19.81),model)
        add_agent_pos!(Agent4(4, (4,1), 3),model)

        props = [:pos, :weight, :p, wpos]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props)
        @test size(df) == (4, 7)
        @test typeof(df.pos) <: Vector{Tuple{Int, Int}}
        @test typeof(df.weight) <: Vector{Union{Missing, Float64}}
        @test typeof(df.p) <: Vector{Union{Missing, Int}}
        @test typeof(df.wpos) <: Vector{Union{Missing, Float64}}

        # Expect something completely unknown to fail
        props = [:UNKNOWN]
        @test_throws ErrorException init_agent_dataframe(model, props)

        # Aggregates should behave in a similiar fashion
        pos1(a) = a.pos[1]
        props = [(:pos, length), (pos1, sum)]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props)
        @test size(df) == (1, 3)
        @test typeof(df.length_pos) <: Vector{Int}
        @test df[1, :length_pos] == 4
        @test typeof(df.sum_pos1) <: Vector{Int}
        @test df[1, :sum_pos1] == 22

        # Expect aggregate to filter out Agent4's
        a3(a) = a isa Agent3
        props = [(pos1, sum, a3)]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props)
        @test size(df) == (1, 2)
        @test typeof(df.sum_pos1_a3) <: Vector{Int}
        @test df[1, :sum_pos1_a3] == 8

        # When aggregating, missing data must be handled explicitly.
        props = [(:weight, sum)]
        df = init_agent_dataframe(model, props)
        @test_throws ErrorException collect_agent_data!(df, model, props)

        # Filtering out missings makes this work.
        props = [(:weight, sum, a3)]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props)
        @test size(df) == (1, 2)
        @test typeof(df.sum_weight_a3) <: Vector{Float64}
        @test df[1, :sum_weight_a3] ≈ 74.46

        # Handle missmatches
        # In this example, weight exists in both agents, but they have different types
        mutable struct Agent3Int <: AbstractAgent
            id::Int
            pos::Dims{2}
            weight::Int
        end
        model = ABM(Union{Agent3, Agent3Int}, GridSpace((10,10)); warn = false)
        add_agent_pos!(Agent3(1, (6,8), 54.65),model)
        add_agent_pos!(Agent3Int(2, (10,7), 5),model)
        add_agent_pos!(Agent3(3, (2,4), 19.81),model)
        add_agent_pos!(Agent3Int(4, (4,1), 3),model)

        props = [:weight]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props)
        @test size(df) == (4, 4)
        @test typeof(df.weight) <: Vector{Union{Float64, Int}}

        # Expect a1.weight <: Float64, a2.weight <: Int64 to fail in aggregate
        props = [(:weight, sum)]
        @test_throws ErrorException init_agent_dataframe(model, props)

        # Promotion is the fix in this case
        fweight(a) = Float64(a.weight)
        props = [(fweight, sum)]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props)
        @test size(df) == (1, 2)
        @test df[1, :sum_fweight] ≈ 82.46
    end

    @testset "Observers" begin
        model = initialize()
        model_props = [:container]
        model_data = init_model_dataframe(model, model_props)
        push!(model.container, 50.0)
        collect_model_data!(model_data, model, model_props, 0; obtainer = copy)
        push!(model.container, 37.2)
        collect_model_data!(model_data, model, model_props, 1; obtainer = copy)
        model.container[1] += 21.9
        collect_model_data!(model_data, model, model_props, 2; obtainer = copy)
        @test model_data.container[1][1] ≈ 50.0
        @test model_data.container[3][1] ≈ 71.9
        @test length.(model_data.container) == [1, 2, 2]

        model = initialize()
        model_props = [:deep]
        model_data = init_model_dataframe(model, model_props)
        push!(model.deep.data, 17.5)
        collect_model_data!(model_data, model, model_props, 0; obtainer = deepcopy)
        push!(model.deep.data, 1.2)
        collect_model_data!(model_data, model, model_props, 1; obtainer = deepcopy)
        model.deep.data[1] += 0.9
        collect_model_data!(model_data, model, model_props, 2; obtainer = deepcopy)
        @test model_data[1,:deep].data[1] ≈ 20.0
        @test model_data[3,:deep].data[1] ≈ 20.9
        @test [length(d.data) for d in model_data[!,:deep]] == [3, 4, 4]

        model = initialize()
        agent_data, model_data = run!(
            model,
            agent_step!,
            model_step!,
            365 * 5;
            when_model = [1, 365 * 5],
            when = false,
            mdata = [:flag, :year, :container, :deep],
            adata = [(:weight, mean)],
            obtainer = deepcopy,
        )
        @test size(agent_data) == (0, 2)
        @test size(model_data) == (2, 5)
    end
end

@testset "Parameter scan" begin
    n = 10
    parameters = Dict(
        :density => [0.6, 0.7, 0.8],
        :griddims => (20, 20),
    )

    forest, agent_step!, forest_step! = Models.forest_fire()
    forest_initiation(; kwargs...) = Models.forest_fire(; kwargs...)[1]

    burnt(a) = a.status == :burnt
    unburnt(a) = a.status == :green
    @testset "Standard Scan" begin
        adata = [(unburnt, count), (burnt, count)]
        data, _ = paramscan(
            parameters,
            forest_initiation;
            n = n,
            agent_step! = agent_step!,
            model_step! = forest_step!,
            adata = adata,
            progress = false,
        )
        # 3 is the number of combinations of changing params
        @test size(data) == ((n + 1) * 3, 4)
        data, _ = paramscan(
            parameters,
            forest_initiation;
            n = n,
            agent_step! = agent_step!,
            model_step! = forest_step!,
            include_constants = true,
            adata = adata,
            progress = false,
        )
        # 3 is the number of combinations of changing params,
        # 5 is 3+2, where 2 is the number of constant parameters
        @test size(data) == ((n + 1) * 3, 5)

        adata = [:status]
        data, _ = paramscan(
            parameters,
            forest_initiation;
            n = n,
            agent_step! = agent_step!,
            model_step! = forest_step!,
            adata = adata,
            progress = false,
        )
        @test unique(data.step) == 0:10
        @test unique(data.density) == [0.6, 0.7, 0.8]
    end

    @testset "Scan with replicates" begin
        adata = [(unburnt, count), (burnt, count)]
        data, _ = paramscan(
            parameters,
            forest_initiation;
            n = n,
            agent_step! = dummystep,
            model_step! = forest_step!,
            replicates = 3,
            adata = adata,
            progress = false,
        )
        # the first 3 is the number of combinations of changing params
        @test size(data) == (((n + 1) * 3) * 3, 5)
    end
end

@testset "Issue #179 fix" begin
    # only ids sorted, not properties
    model = ABM(Agent2)
    for i in 1:5
        add_agent!(model, i * 0.2)
    end
    data, _ = run!(model, dummystep, 2; adata = [:weight])
    @test data[1, :id] == 1 && data[1, :weight] ≈ 0.2
    @test data[3, :id] == 3 && data[3, :weight] ≈ 0.6
    @test data[6, :id] == 1 && data[6, :weight] ≈ 0.2
end


text:
import UnicodePlots, Plots
using  Printf, LoopVectorization, LinearAlgebra, SparseArrays

include("CreateMeshFCFV.jl")
include("VisuFCFV.jl")
include("DiscretisationFCFV_Stokes.jl")
include("SolversFCFV_Stokes.jl")
include("EvalAnalDani.jl")

#--------------------------------------------------------------------#

function SetUpProblem!(mesh, P, Vx, Vy, Sxx, Syy, Sxy, VxDir, VyDir, SxxNeu, SyyNeu, SxyNeu, SyxNeu, sx, sy, R, eta, gbar)
    # Evaluate T analytic on cell faces
    etam, etai = eta[1], eta[2]
    for in=1:mesh.nf
        x        = mesh.xf[in]
        y        = mesh.yf[in]
        vx, vy, pre, sxx, syy, sxy = EvalAnalDani( x, y, R, etam, etai )
        VxDir[in] = vx
        VyDir[in] = vy
        # Stress at faces
        p          = pre
        # Pseudo-tractions
        pre, dVxdx, dVxdy, dVydx, dVydy = Tractions( x, y, R, etam, etai, 1 )
        SxxNeu[in] = - pre + etam*dVxdx
        SyyNeu[in] = - pre + etam*dVydy
        SxyNeu[in] =         etam*dVxdy
        SyxNeu[in] =         etam*dVydx
    end
    # Evaluate T analytic on barycentres
    for iel=1:mesh.nel
        x        = mesh.xc[iel]
        y        = mesh.yc[iel]
        vx, vy, pre, sxx, syy, sxy = EvalAnalDani( x, y, R, etam, etai )
        P[iel]   = pre
        Vx[iel]  = vx
        Vy[iel]  = vy
        Sxx[iel] = sxx
        Syy[iel] = syy
        Sxy[iel] = sxy
        sx[iel]  = 0.0
        sy[iel]  = 0.0
        out          = mesh.phase[iel] == 1.0
        mesh.ke[iel] = (out==1) * 1.0*etam + (out!=1) * 1.0*etai
        # Compute jump condition
        for ifac=1:mesh.nf_el
            # Face
            nodei  = mesh.e2f[iel,ifac]
            xF     = mesh.xf[nodei]
            yF     = mesh.yf[nodei]
            nodei  = mesh.e2f[iel,ifac]
            bc     = mesh.bc[nodei]
            dAi    = mesh.dA[iel,ifac]
            ni_x   = mesh.n_x[iel,ifac]
            ni_y   = mesh.n_y[iel,ifac]
            phase1 = Int64(mesh.phase[iel])
            pre1, dVxdx1, dVxdy1, dVydx1, dVydy1 = Tractions( xF, yF, R, etam, etai, phase1 )
            eta_face1 = eta[phase1]
            tL_x  = ni_x*(-pre1 + eta_face1*dVxdx1) + ni_y*eta_face1*dVxdy1
            tL_y  = ni_y*(-pre1 + eta_face1*dVydy1) + ni_x*eta_face1*dVydx1
            # From element 2
            ineigh = (mesh.e2e[iel,ifac]>0) * mesh.e2e[iel,ifac] + (mesh.e2e[iel,ifac]<1) * iel
            phase2 = Int64(mesh.phase[ineigh])
            pre2, dVxdx2, dVxdy2, dVydx2, dVydy2 = Tractions( xF, yF, R, etam, etai, phase2)
            eta_face2 = eta[phase2]
            tR_x = ni_x*(-pre2 + eta_face2*dVxdx2) + ni_y*eta_face2*dVxdy2
            tR_y = ni_y*(-pre2 + eta_face2*dVydy2) + ni_x*eta_face2*dVydx2
            gbar[iel,ifac,1] = 0.5 * (tL_x - tR_x)
            gbar[iel,ifac,2] = 0.5 * (tL_y - tR_y)
        end
    end

    # # detect cells at interface and set to geometric mean
    # for iel=1:mesh.nel

    #     for ifac=1:mesh.nf_el
    #         nodei  = mesh.e2f[iel,ifac]
    #         if mesh.bc[nodei] == 3
    #             # mesh.ke[iel] = (etam*etai)^(1/2)
    #             if Int64(mesh.phase[iel]) == 1
    #                 mesh.ke[iel] = 2.0/(1.0/etam+1.0/etai)
    #             else
    #                 mesh.ke[iel] = 0.5*(etam+etai)
    #             end
    #         end
    #     end
    # end

    return
end

#--------------------------------------------------------------------#

function ComputeError( mesh, Vxe, Vye, Txxe, Tyye, Txye, Pe, R, eta )
    etam, etai = eta[1], eta[2]
    eVx  = zeros(mesh.nel)
    eVy  = zeros(mesh.nel)
    eTxx = zeros(mesh.nel)
    eTyy = zeros(mesh.nel)
    eTxy = zeros(mesh.nel)
    eP   = zeros(mesh.nel)
    Vxa  = zeros(mesh.nel)
    Vya  = zeros(mesh.nel)
    Txxa = zeros(mesh.nel)
    Tyya = zeros(mesh.nel)
    Txya = zeros(mesh.nel)
    Pa   = zeros(mesh.nel)
    for iel=1:mesh.nel
        x         = mesh.xc[iel]
        y         = mesh.yc[iel]
        vx, vy, pre, sxx, syy, sxy = EvalAnalDani( x, y, R, etam, etai )
        pre, dVxdx, dVxdy, dVydx, dVydy = Tractions( x, y, R, etam, etai, mesh.phase[iel] )
        Pa[iel]   = pre
        Vxa[iel]  = vx
        Vya[iel]  = vy
        Txxa[iel] = 2.0*mesh.ke[iel]*dVxdx
        Tyya[iel] = 2.0*mesh.ke[iel]*dVydy
        Txya[iel] = 1.0*mesh.ke[iel]*(dVxdy+dVydx)
        eVx[iel]  = Vxe[iel]  - Vxa[iel]
        eVy[iel]  = Vye[iel]  - Vya[iel]
        eTxx[iel] = Txxe[iel] - Txxa[iel]
        eTyy[iel] = Tyye[iel] - Tyya[iel]
        eTxy[iel] = Txye[iel] - Txya[iel]
        eP[iel]   = Pe[iel]   - Pa[iel]
    end
    errVx  = norm(eVx) /norm(Vxa)
    errVy  = norm(eVy) /norm(Vya)
    errTxx = norm(eTxx)/norm(Txxa)
    errTyy = norm(eTyy)/norm(Tyya)
    errTxy = norm(eTxy)/norm(Txya)
    errP   = norm(eP)  /norm(Pa)
    return errVx, errVy, errTxx, errTyy, errTxy, errP, Txxa, Tyya, Txya
end

#--------------------------------------------------------------------#

function StabParam(tau, dA, Vol, mesh_type, coeff)
    if mesh_type=="Quadrangles";        taui = coeff*tau  end
    if mesh_type=="UnstructTriangles";  taui = coeff*tau*dA  end
    return taui
end

#--------------------------------------------------------------------#

@views function main( n, solver)
    # ```This version include the jump condition derived from the analytical solution
    # The great thing is that pressure converges in L_infinity norm evem with quadrangles - this rocks
    # # Test #1: For a contrast of:  eta        = [10.0 1.0] (weak inclusion), tau = 20.0
    #     res = [1 2 4 8]*30 elements per dimension
    #     err = [4.839 3.944 2.811 1.798]
    # # Test #2: For a contrast of:  eta        = [1.0 100.0] (strong inclusion), tau = 1.0/5.0
    #     res = [1 2 4 8]*30 elements per dimension
    #     err = [4.169 2.658 1.628 0.83]
    # ```

    println("\n\n******** FCFV STOKES ********")

    # Create sides of mesh
    xmin, xmax = -3.0, 3.0
    ymin, ymax = -3.0, 3.0
    # n          = 16
    nx, ny     = 30*n, 30*n
    # solver     = 3
    R          = 1.0
    inclusion  = 1
    eta        = [1.0 100.0]
    # mesh_type  = "Quadrangles"
    mesh_type  = "UnstructTriangles"
    BC         = [1; 1; 1; 1] # S E N W --- 1: Dirichlet / 2: Neumann

    # Generate mesh
    if mesh_type=="Quadrangles"
        tau  = 4.0
        mesh = MakeQuadMesh( nx, ny, xmin, xmax, ymin, ymax, inclusion, R, BC )
    elseif mesh_type=="UnstructTriangles"
        tau  = 1.0/5.0
        # area = 1.0 # area factor: SETTING REPRODUCE THE RESULTS OF MATLAB CODE USING TRIANGLE
        # ninc = 29  # number of points that mesh the inclusion: SETTING REPRODUCE THE RESULTS OF MATLAB CODE USING TRIANGLE
        mesh = MakeTriangleMesh( nx, ny, xmin, xmax, ymin, ymax, inclusion, R, BC )
    end
    println("Number of elements: ", mesh.nel)

    # Source term and BCs etc...
    Pa     = zeros(mesh.nel)
    Vxa    = zeros(mesh.nel)
    Vya    = zeros(mesh.nel)
    Sxxa   = zeros(mesh.nel)
    Syya   = zeros(mesh.nel)
    Sxya   = zeros(mesh.nel)
    sex    = zeros(mesh.nel)
    sey    = zeros(mesh.nel)
    VxDir  = zeros(mesh.nf)
    VyDir  = zeros(mesh.nf)
    SxxNeu = zeros(mesh.nf)
    SyyNeu = zeros(mesh.nf)
    SxyNeu = zeros(mesh.nf)
    SyxNeu = zeros(mesh.nf)
    gbar   = zeros(mesh.nel,mesh.nf_el,2)
    println("Model configuration :")
    @time SetUpProblem!(mesh, Pa, Vxa, Vya, Sxxa, Syya, Sxya, VxDir, VyDir, SxxNeu, SyyNeu, SxyNeu, SyxNeu, sex, sey, R, eta, gbar)

    # Compute some mesh vectors
    println("Compute FCFV vectors:")
    @time ae, be, ze = ComputeFCFV(mesh, sex, sey, VxDir, VyDir, SxxNeu, SyyNeu, SxyNeu, SyxNeu, tau)

    # Assemble element matrices and RHS
    println("Compute element matrices:")
    @time Kuu_v, fu_v, Kup_v, fp = ElementAssemblyLoop(mesh, ae, be, ze, VxDir, VyDir, SxxNeu, SyyNeu, SxyNeu, SyxNeu, gbar, tau)

    # Assemble triplets and sparse
    println("Assemble triplets and sparse:")
    @time Kuu, fu, Kup = CreateTripletsSparse(mesh, Kuu_v, fu_v, Kup_v)

    # display(UnicodePlots.spy(Kuu))
    # display(UnicodePlots.spy(Kup))

    # Solve for hybrid variable
    println("Linear solve:")
    t_solve = @elapsed  Vxh, Vyh, Pe = StokesSolvers(mesh, Kuu, Kup, fu, fp, solver)
    println("Linear solver took: ", t_solve, " s")

    # # Reconstruct element values
    println("Compute element values:")
    @time Vxe, Vye, Txxe, Tyye, Txye = ComputeElementValues(mesh, Vxh, Vyh, Pe, ae, be, ze, VxDir, VyDir, tau)

    # # Compute discretisation errors
    err_Vx, err_Vy, err_Txx, err_Tyy, err_Txy, err_P, Txxa, Tyya, Txya = ComputeError( mesh, Vxe, Vye, Txxe, Tyye, Txye, Pe, R, eta )
    @printf("Error in Vx : %2.2e\n", err_Vx )
    @printf("Error in Vy : %2.2e\n", err_Vy )
    @printf("Error in Txx: %2.2e\n", err_Txx)
    @printf("Error in Tyy: %2.2e\n", err_Tyy)
    @printf("Error in Txy: %2.2e\n", err_Txy)
    @printf("Error in P  : %2.2e\n", err_P  )

    Perr = abs.(Pa.-Pe)
    Verr = sqrt.( (Vxe.-Vxa).^2 .+ (Vye.-Vya).^2 )
    println("L_inf P error: ", maximum(Perr), " --- L_inf V error: ", maximum(Verr))

    # # Visualise
    # println("Visualisation:")
    # # @time PlotMakie(mesh, Vxe, xmin, xmax, ymin, ymax, :jet1)
    # # @time PlotMakie( mesh, Verr, xmin, xmax, ymin, ymax, :jet1, minimum(Verr), maximum(Verr) )
    # # @time PlotMakie( mesh, Pe, xmin, xmax, ymin, ymax, :jet1, minimum(Pa), maximum(Pa) )
    # # @time PlotMakie( mesh, Perr, xmin, xmax, ymin, ymax, :jet1, minimum(Perr), maximum(Perr) )
    # # @time PlotMakie( mesh, Txxe, xmin, xmax, ymin, ymax, :jet1, -6.0, 2.0 )
    @time PlotMakie( mesh, Tyya, xmin, xmax, ymin, ymax, :jet1, -2.0, 6.0 )
    # # @time PlotMakie( mesh, log10.(mesh.ke), xmin, xmax, ymin, ymax, :jet1 )
    # # @time PlotMakie( mesh, mesh.phase, xmin, xmax, ymin, ymax, :jet1)

    ndof = 2*mesh.nf+mesh.nel
    print(t_solve)
    return t_solve, ndof

end

N = 1:1
ndofs = zeros(length(N) )
t0    = zeros(length(N) )
t1    = zeros(length(N) )
t3    = zeros(length(N) )

main( 1, 1 ) # warm up

for i=1:length(N)
    # t0[i], ndofs[i] = main( N[i], 0 )
    # t1[i], ndofs[i]          = main( N[i], 1 )
    t3[i] , ndofs[i]          = main( N[i], 3 )
end

# p = Plots.plot(  ndofs,  t0,  markershape=:dtriangle, label="Monolithic", legend=:bottomright, xlabel = "ndof", ylabel = "solve time" )
# p = Plots.plot( ndofs,  t1,  markershape=:dtriangle, label="Powell-Hestenes", legend=:bottomright, xlabel = "ndof", ylabel = "solve time" )
p = Plots.plot( ndofs,  t3,  markershape=:dtriangle, label="Modified PH", legend=:bottomright, xlabel = "ndof", ylabel = "solve time" )
display(p)


text:
""" Returns the baseline OptimalTreeLearners, with all valid training args. """
function base_classifier()
    IAI.OptimalTreeClassifier(
        hyperplane_config = (sparsity = :all,),
        random_seed = 1,
        max_depth = 5,
        cp = 1e-6,
        minbucket = 0.01,
        fast_num_support_restarts = 5,
        localsearch = true,
        ls_num_tree_restarts = 10,
        ls_num_hyper_restarts =  5)
end

function base_regressor()
    IAI.OptimalTreeRegressor(
        hyperplane_config = (sparsity = :all,),
        regression_sparsity = :all,
        regression_weighted_betas = true,
        random_seed = 1,
        max_depth = 5,
        cp = 1e-6,
        minbucket = 0.02,
        fast_num_support_restarts = 5,
        regression_lambda = 1e-5,
        ls_num_tree_restarts = 10,
        ls_num_hyper_restarts =  5)
end

function base_rf_regressor()
    IAI.Heuristics.RandomForestRegressor(
        num_trees = 20,
        hyperplane_config = (sparsity = :all,),
        regression_sparsity = 0,
        regression_weighted_betas = true,
        random_seed = 1,
        max_depth = 3,
        cp = 1e-6,
        minbucket = 0.02,
        fast_num_support_restarts = 5,
        regression_lambda = 1e-5,
        ls_num_tree_restarts = 10,
        ls_num_hyper_restarts =  5)
end

function base_rf_classifier()
    IAI.Heuristics.RandomForestClassifier(
        num_trees = 20,
        hyperplane_config = (sparsity = :all,),
        random_seed = 1,
        max_depth = 3,
        cp = 1e-6,
        minbucket = 0.02,
        fast_num_support_restarts = 5,
        ls_num_tree_restarts = 10,
        ls_num_hyper_restarts =  5)
end


text:
test_a1 = """
H ANDN920101
D alpha-CH chemical shifts (Andersen et al., 1992)
R LIT:1810048b PMID:1575719
A Andersen, N.H., Cao, B. and Chen, C.
T Peptide/protein structure analysis using the chemical shift index method:
upfield alpha-CH values reveal dynamic helices and aL sites
J Biochem. and Biophys. Res. Comm. 184, 1008-1014 (1992)
C BUNA790102    0.949
I    A/L     R/K     N/M     D/F     C/P     Q/S     E/T     G/W     H/Y     I/V
    4.35    4.38    4.75    4.76    4.65    4.37    4.29    3.97    4.63    3.95
    4.17    4.36    4.52    4.66    4.44    4.50    4.35    4.70    4.60    3.95
//
"""

test_a2 = """
H ALTS910101
D The PAM-120 matrix (Altschul, 1991)
R LIT:1713145 PMID:2051488
A Altschul, S.F.
T Amino acid substitution matrices from an information theoretic perspective
J J. Mol. Biol. 219, 555-565 (1991)
M rows = ARNDCQEGHILKMFPSTWYV, cols = ARNDCQEGHILKMFPSTWYV
    3.
    -3.      6.
    0.     -1.      4.
    0.     -3.      2.      5.
    -3.     -4.     -5.     -7.      9.
    -1.      1.      0.      1.     -7.      6.
    0.     -3.      1.      3.     -7.      2.      5.
    1.     -4.      0.      0.     -5.     -3.     -1.      5.
    -3.      1.      2.      0.     -4.      3.     -1.     -4.      7.
    -1.     -2.     -2.     -3.     -3.     -3.     -3.     -4.     -4.      6.
    -3.     -4.     -4.     -5.     -7.     -2.     -4.     -5.     -3.      1.      5.
    -2.      2.      1.     -1.     -7.      0.     -1.     -3.     -2.     -2.     -4.      5.
    -2.     -1.     -3.     -4.     -6.     -1.     -4.     -4.     -4.      1.      3.      0.      8.
    -4.     -4.     -4.     -7.     -6.     -6.     -6.     -5.     -2.      0.      0.     -6.     -1.      8.
    1.     -1.     -2.     -2.     -3.      0.     -1.     -2.     -1.     -3.     -3.     -2.     -3.     -5.      6.
    1.     -1.      1.      0.     -1.     -2.     -1.      1.     -2.     -2.     -4.     -1.     -2.     -3.      1.      3.
    1.     -2.      0.     -1.     -3.     -2.     -2.     -1.     -3.      0.     -3.     -1.     -1.     -4.     -1.      2.      4.
    -7.      1.     -5.     -8.     -8.     -6.     -8.     -8.     -5.     -7.     -5.     -5.     -7.     -1.     -7.     -2.     -6.     12.
    -4.     -6.     -2.     -5.     -1.     -5.     -4.     -6.     -1.     -2.     -3.     -6.     -4.      4.     -6.     -3.     -3.     -1.      8.
    0.     -3.     -3.     -3.     -2.     -3.     -3.     -2.     -3.      3.      1.     -4.      1.     -3.     -2.     -2.      0.     -8.     -3.      5.
//
"""

test_a3 = """
H TANS760101
D Statistical contact potential derived from 25 x-ray protein structures
R PMID:1004017
A Tanaka, S. and Scheraga, H.A.
T Medium- and long-range interaction parameters between amino acids
for predicting three-dimensional structures of proteins
J Macromolecules 9, 945-950 (1976)
M rows = ARNDCQEGHILKMFPSTWYV, cols = ARNDCQEGHILKMFPSTWYV
-2.6
-3.4 -4.3
-3.1 -4.1 -3.2
-2.8 -3.9 -3.1 -2.7
-4.2 -5.3 -4.9 -4.2 -7.1
-3.5 -4.5 -3.8 -3.2 -5.0 -3.4
-3.0 -4.2 -3.4 -3.3 -4.4 -3.6 -2.8
-3.8 -4.5 -4.0 -3.7 -5.4 -4.4 -3.8 -3.9
-4.0 -4.9 -4.4 -4.3 -5.6 -4.7 -4.5 -4.7 -4.9
-5.9 -6.2 -5.8 -5.4 -7.3 -5.9 -5.7 -6.3 -6.6 -8.2
-4.8 -5.1 -4.6 -4.3 -6.2 -5.0 -4.6 -5.2 -5.6 -7.5 -6.0
-3.1 -3.6 -3.3 -3.2 -4.4 -3.7 -3.8 -3.8 -4.1 -5.6 -4.6 -2.7
-4.6 -5.0 -4.2 -4.3 -6.2 -3.5 -4.6 -5.1 -5.4 -7.4 -6.3 -4.7 -5.8
-5.1 -5.8 -5.0 -4.9 -6.8 -5.3 -5.0 -5.6 -6.4 -8.0 -7.0 -4.9 -6.6 -7.1
-3.4 -4.2 -3.6 -3.3 -5.3 -4.0 -3.5 -4.2 -4.5 -6.0 -4.8 -3.6 -5.1 -5.2 -3.5
-2.9 -3.8 -3.1 -2.7 -4.6 -3.6 -3.2 -3.8 -4.3 -5.5 -4.4 -3.0 -4.1 -4.7 -3.4 -2.5
-3.3 -4.0 -3.5 -3.1 -4.8 -3.7 -3.3 -4.1 -4.5 -5.9 -4.8 -3.3 -4.6 -5.1 -3.6 -3.3 -3.1
-5.2 -5.8 -5.3 -5.1 -6.9 -5.8 -5.2 -5.8 -6.5 -7.8 -6.8 -5.0 -6.9 -7.4 -5.6 -5.0 -5.1 -6.8
-4.7 -5.6 -5.0 -4.7 -6.6 -5.2 -4.9 -5.4 -6.1 -7.4 -6.2 -4.9 -6.1 -6.6 -5.2 -4.7 -4.9 -6.8 -6.0
-4.3 -4.9 -4.3 -4.0 -6.0 -4.7 -4.2 -5.1 -5.3 -7.3 -6.2 -4.2 -6.0 -6.5 -4.7 -4.2 -4.4 -6.5 -5.9 -5.5
//
"""

test_index= """
I    A/L     R/K     N/M     D/F     C/P     Q/S     E/T     G/W     H/Y     I/V
    4.35    4.38    4.75    4.76    4.65    4.37    4.29    3.97    4.63    3.95
    4.17    4.36    4.52    4.66    4.44    4.50    4.35    4.70    4.60    3.95
"""

test_index_result = [
    4.35, 4.38, 4.75, 4.76, 4.65, 4.37, 4.29, 3.97, 4.63, 3.95, 4.17, 4.36,
    4.52, 4.66, 4.44, 4.50, 4.35, 4.70, 4.60, 3.95
]

test_lower_triangular_matrix = """
M rows = ARNDCQEGHILKMFPSTWYV, cols = ARNDCQEGHILKMFPSTWYV
-2.6
-3.4 -4.3
-3.1 -4.1 -3.2
-2.8 -3.9 -3.1 -2.7
-4.2 -5.3 -4.9 -4.2 -7.1
-3.5 -4.5 -3.8 -3.2 -5.0 -3.4
-3.0 -4.2 -3.4 -3.3 -4.4 -3.6 -2.8
-3.8 -4.5 -4.0 -3.7 -5.4 -4.4 -3.8 -3.9
-4.0 -4.9 -4.4 -4.3 -5.6 -4.7 -4.5 -4.7 -4.9
-5.9 -6.2 -5.8 -5.4 -7.3 -5.9 -5.7 -6.3 -6.6 -8.2
-4.8 -5.1 -4.6 -4.3 -6.2 -5.0 -4.6 -5.2 -5.6 -7.5 -6.0
-3.1 -3.6 -3.3 -3.2 -4.4 -3.7 -3.8 -3.8 -4.1 -5.6 -4.6 -2.7
-4.6 -5.0 -4.2 -4.3 -6.2 -3.5 -4.6 -5.1 -5.4 -7.4 -6.3 -4.7 -5.8
-5.1 -5.8 -5.0 -4.9 -6.8 -5.3 -5.0 -5.6 -6.4 -8.0 -7.0 -4.9 -6.6 -7.1
-3.4 -4.2 -3.6 -3.3 -5.3 -4.0 -3.5 -4.2 -4.5 -6.0 -4.8 -3.6 -5.1 -5.2 -3.5
-2.9 -3.8 -3.1 -2.7 -4.6 -3.6 -3.2 -3.8 -4.3 -5.5 -4.4 -3.0 -4.1 -4.7 -3.4 -2.5
-3.3 -4.0 -3.5 -3.1 -4.8 -3.7 -3.3 -4.1 -4.5 -5.9 -4.8 -3.3 -4.6 -5.1 -3.6 -3.3 -3.1
-5.2 -5.8 -5.3 -5.1 -6.9 -5.8 -5.2 -5.8 -6.5 -7.8 -6.8 -5.0 -6.9 -7.4 -5.6 -5.0 -5.1 -6.8
-4.7 -5.6 -5.0 -4.7 -6.6 -5.2 -4.9 -5.4 -6.1 -7.4 -6.2 -4.9 -6.1 -6.6 -5.2 -4.7 -4.9 -6.8 -6.0
-4.3 -4.9 -4.3 -4.0 -6.0 -4.7 -4.2 -5.1 -5.3 -7.3 -6.2 -4.2 -6.0 -6.5 -4.7 -4.2 -4.4 -6.5 -5.9 -5.5
"""

test_non_triangular_matrix = """
M rows = ARNDCQEGHILKMFPSTWYV, cols = ARNDCQEGHILKMFPSTWYV
  -0.94  1.26  0.55 0.76 -1.54  1.14 1.57 -0.78  0.44 -1.59 -1.64  1.91 -0.90 -1.49  0.28  0.20 -0.04  -0.92 -0.75 -1.45
   0.56  1.79  2.31 0.79 -0.67  2.54 0.72  1.09  0.94 -0.01  0.01  3.68  0.89 -0.05  1.37  0.83  1.35   0.00  0.33  0.44
   0.59  2.21  1.82 0.77 -0.90  0.46 3.06 -0.16  0.63 -0.33  0.20  2.43  0.99  0.63  0.54  0.24  0.63   0.11 -0.19  0.23
   0.66  0.76  0.76 1.19 -0.21  1.66 2.22  0.29  0.57  0.59  0.79  1.13  1.41  0.49  1.70  1.03  1.19   1.85  0.18  0.86
  -1.75  0.78 -1.00 0.32 -3.64  0.48 0.87 -1.67 -0.62 -2.77 -2.32  0.19 -1.22 -2.67 -1.62 -0.83 -1.14  -0.52 -1.94 -2.35
   0.33  2.15  1.22 1.26  1.37  1.17 2.56  0.92  1.02  0.11  0.00  2.58  0.79 -0.26  0.53  1.19  1.11   0.21  0.39  0.15
   0.82  1.05  2.18 2.11  0.01  2.42 2.58  1.15  0.97  0.20  0.31  1.31  1.25  0.12  2.00  1.09  1.13   0.58  0.31  0.39
  -0.40  0.95  0.03 0.14 -1.00  0.34 0.99 -1.32  0.13 -1.40 -1.36  1.58 -0.90 -1.41  0.82 -0.27  0.21  -0.59 -1.27 -1.09
  -0.75  2.19  0.13 0.68 -1.37  1.98 1.13  0.01  1.52 -0.83 -0.58  2.26 -0.82 -1.01  0.53 -0.17  0.02 -49.00 -0.61 -0.56
  -1.99  0.25 -0.20 1.00 -2.44 -0.12 0.88 -1.54 -0.05 -2.64 -2.33  0.75 -1.85 -2.46 -1.06 -0.59 -0.65  -1.82 -1.88 -2.45
  -2.02  0.34 -0.04 0.13 -2.29  0.24 0.73 -1.27 -0.46 -2.53 -2.44  0.67 -1.80 -2.28 -1.29 -0.40 -0.34  -1.76 -1.66 -2.26
   0.60  3.11  2.23 1.06  0.50  1.80 1.65  0.82  1.25  0.10  0.34  3.51  0.98 -0.21  1.15  2.09  1.30  -0.14  0.28  0.13
  -1.54 -0.06 -0.63 1.76 -2.51  0.14 0.72 -1.74  0.07 -2.27 -2.22  1.27 -1.77 -1.87  0.34 -0.02 -0.21  -0.93 -1.54 -1.81
  -2.12  0.33 -0.70 0.17 -2.30 -0.59 0.26 -1.60 -0.88 -2.53 -2.44 -0.42 -1.83 -2.68 -1.40 -0.82 -0.61  -1.63 -1.83 -2.25
   0.63  2.43 -0.19 1.31 -1.63  1.46 1.91  0.08  1.11 -0.20  0.47  1.94 -0.34  0.15  0.57  0.00  1.15   0.06  0.26 -0.06
  -0.41  0.88  1.02 1.04 -0.21  1.27 0.94  0.04  0.75 -0.48 -0.67  2.28  0.45 -0.92  0.75  0.50  0.96   0.22 -0.19 -0.54
  -0.32  1.48  0.35 0.43 -1.44  0.38 1.36 -0.38  0.20 -1.14 -1.00  1.38 -0.35 -0.97 -0.05 -0.16  0.29  -0.53 -0.76 -0.73
  -1.85  0.45 -0.03 0.80 -1.64 -0.23 0.11 -0.95  0.67 -1.58 -2.13  0.61 -1.75 -1.59 -1.07 -0.34 -0.40  -1.29 -1.27 -1.79
  -0.88 -0.20 -0.29 0.14 -1.31  0.09 0.71 -0.56 -0.57 -1.66 -1.38  1.40 -1.60 -1.97 -0.73 -0.32 -0.37  -1.40 -0.96 -1.38
  -1.74  0.85  0.24 0.72 -2.25  0.45 0.81 -1.29 -0.24 -2.46 -2.38  0.37 -1.21 -2.16 -1.00 -0.10 -0.57  -1.34 -1.52 -2.31
"""


text:
include("../src/quasirandom_module.jl")

module qtest
using ..quasirandom
using PyPlot
PyPlot.pygui(true)

# h = low.halton(1; base = 2)
# h = low.halton(10; base = 2)

n_points = 21
n_dim = 2
hseq = quasirandom.halton_filling_seq(n_points, n_dim)

mmc = quasirandom.morris_mitchell_criterion(hseq)

PyPlot.figure()
PyPlot.grid()
PyPlot.plot(hseq[:,1], hseq[:,2], linestyle = :none, marker = :., markersize = 10)

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
@testset "$(rpad("Symplecticity",80))" begin

    for s in 1:10
        g = TableauGauss(s)
        g̃ = get_symplectic_conjugate_coefficients(g)

        @test g.a ≈ g̃.a atol=1E-15
        @test g.b == g̃.b
        @test g.c == g̃.c

        @test compute_symplecticity_error(g) ≈ zeros(s,s)  atol=eps()
        @test check_symplecticity(g) == Array{Bool}(ones(s,s))
        @test issymplectic(g)
    end

    for s in 2:10
        A = TableauLobattoIIIA(s)
        B = TableauLobattoIIIB(s)
        E = TableauLobattoIIIE(s)
        Ã = get_symplectic_conjugate_coefficients(A)
        B̃ = get_symplectic_conjugate_coefficients(B)

        @test A.a ≈ B̃.a atol=1E-15
        @test A.b == Ã.b
        @test A.c == Ã.c

        @test B.a ≈ Ã.a atol=1E-15
        @test B.b == Ã.b
        @test B.c == Ã.c

        Â = symplecticize(A)
        B̂ = symplecticize(B)

        @test Â.a ≈ E.a atol=1E-15
        @test Â.b == E.b
        @test Â.c == E.c

        @test B̂.a ≈ E.a atol=1E-15
        @test B̂.b == E.b
        @test B̂.c == E.c

        @test  issymplectic(Â)
        @test  issymplectic(B̂)
    end

    for s in 2:10
        @test !issymplectic(TableauLobattoIIIA(s))
        @test !issymplectic(TableauLobattoIIIB(s))
        @test !issymplectic(TableauLobattoIIIC(s))
        @test !issymplectic(TableauLobattoIIIC̄(s))
        @test  issymplectic(TableauLobattoIIID(s))
        @test  issymplectic(TableauLobattoIIIE(s))
        @test !issymplectic(TableauRadauIA(s))
        @test !issymplectic(TableauRadauIIA(s))
    end

end


text:
using Revise
using Plots
using Pkg
# Pkg.activate(".")
using RadiativeTransfer
using RadiativeTransfer.Absorption
using RadiativeTransfer.Scattering
using RadiativeTransfer.vSmartMOM
using ForwardDiff

# Load parameters from file
parameters = vSmartMOM.parameters_from_yaml("src/vSmartMOM/ModelParameters/DefaultParameters.yaml")

# default_parameters

# Sets all the "specific" parameters
# parameters = vSmartMOM.default_parameters();
function runner(x, parameters=parameters)
    parameters.τAer_ref = [x[1]];
    #@show parameters.p₀
    parameters.p₀ = [x[2]];
    parameters.nᵣ = [x[3]];
    @show parameters.p₀
    model = model_from_parameters(parameters);

    model.params.architecture = RadiativeTransfer.Architectures.CPU()
    J = vSmartMOM.rt_run(model);
    #@show J
    return J[1,1,:]#; R_SFI[1,1,:]
end

#x = [0.1,90001.0]
x = [0.1,90001.0,1.3]
# Run FW model:
@time runner(x);
@time dfdx = ForwardDiff.jacobian(runner, x);


text:
include( "debug.jl" )
include( "bint-type-const.jl" )
include( "bint-bint.jl" )
include( "bint-range.jl" )
include( "bint-arithmetic.jl" )
include( "bint-conversions.jl" )
include( "bint-clamp.jl" )


text:
# MadNLP.jl
# Created by Sungho Shin (sungho.shin@wisc.edu)

module MadNLP

import Pkg.TOML: parsefile
import MathOptInterface
import Libdl: dlopen, dlext, RTLD_DEEPBIND, RTLD_GLOBAL
import Printf: @sprintf
import LinearAlgebra: BLAS, Adjoint, Symmetric, mul!, ldiv!, norm, dot
import LinearAlgebra.BLAS: axpy!, libblas, liblapack, BlasInt, @blasfunc
import SparseArrays: AbstractSparseMatrix, SparseMatrixCSC, sparse, getcolptr, rowvals, nnz
import Logging: @debug, @info,  @warn, @error
import Base: string, show, print, size, getindex, copyto!, @kwdef
import SuiteSparse: UMFPACK
import NLPModels: finalize, AbstractNLPModel, obj, grad!, cons!, jac_coord!, hess_coord!, hess_structure!, jac_structure!, NLPModelMeta, get_nvar, get_ncon, get_minimize, get_x0, get_y0, get_nnzj, get_nnzh, get_lvar, get_uvar, get_lcon, get_ucon, Counters as _Counters # get_zl,get_zu
import SolverCore: AbstractExecutionStats, getStatus

const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities
const NLPModelsCounters = _Counters

export MadNLPUmfpack, MadNLPLapackCPU, MadNLPPardisoMKL, madnlp

# Version info
version() = parsefile(joinpath(@__DIR__,"..","Project.toml"))["version"]
introduce() = "MadNLP version v$(version())"

# Linear solver dependencies
include("enums.jl")
include("utils.jl")
include("options.jl")
include("matrixtools.jl")
include(joinpath("LinearSolvers","linearsolvers.jl"))
include(joinpath("KKT", "KKTsystem.jl"))
include("interiorpointsolver.jl")
include(joinpath("Interfaces","interfaces.jl"))

# Initialize
function __init__()
    # check_deps()
    try
        @isdefined(libpardiso) && dlopen(libpardiso,RTLD_DEEPBIND)
    catch e
        println("Pardiso shared library cannot be loaded")
    end
    set_blas_num_threads(Threads.nthreads(); permanent=true)
end

end # end module