text:
# ------------------------------------------------------------------------------
# BDDC
# ------------------------------------------------------------------------------

# ------------------------------------------------------------------------------
# Lumped
# ------------------------------------------------------------------------------

"""
```julia
LumpedBDDC(operator)
```

Lumped BDDC preconditioner for finite element operators

# Arguments:
- `operator`:  finite element operator to precondition

# Returns:
- Lumped BDDC preconditioner object

# Example:
```jldoctest
# setup
mesh = Mesh2D(1.0, 1.0);

# operator
finediffusion = GalleryOperator("diffusion", 5, 5, mesh);

# preconditioner
bddc = LumpedBDDC(finediffusion);

# verify
println(bddc)
println(bddc.operator)

# output
lumped BDDC preconditioner
finite element operator:
2d mesh:
    dx: 1.0
    dy: 1.0

2 inputs:
operator field:
  tensor product basis:
    numbernodes1d: 5
    numberquadraturepoints1d: 5
    numbercomponents: 1
    dimension: 2
  evaluation mode:
    gradient
operator field:
  tensor product basis:
    numbernodes1d: 5
    numberquadraturepoints1d: 5
    numbercomponents: 1
    dimension: 2
  evaluation mode:
    quadratureweights

1 output:
operator field:
  tensor product basis:
    numbernodes1d: 5
    numberquadraturepoints1d: 5
    numbercomponents: 1
    dimension: 2
  evaluation mode:
    gradient
```
"""
function LumpedBDDC(operator::Operator)
    # common constructor
    return BDDC(operator, BDDCInjectionType.scaled)
end

# ------------------------------------------------------------------------------
# Dirichlet
# ------------------------------------------------------------------------------

"""
```julia
DirichletBDDC(operator)
```

Dirichlet BDDC preconditioner for finite element operators

# Arguments:
- `operator`:  finite element operator to precondition

# Returns:
- Dirichlet BDDC preconditioner object

# Example:
```jldoctest
# setup
mesh = Mesh2D(1.0, 1.0);

# operator
finediffusion = GalleryOperator("diffusion", 5, 5, mesh);

# preconditioner
bddc = DirichletBDDC(finediffusion);

# verify
println(bddc)
println(bddc.operator)

# output
Dirichlet BDDC preconditioner
finite element operator:
2d mesh:
    dx: 1.0
    dy: 1.0

2 inputs:
operator field:
  tensor product basis:
    numbernodes1d: 5
    numberquadraturepoints1d: 5
    numbercomponents: 1
    dimension: 2
  evaluation mode:
    gradient
operator field:
  tensor product basis:
    numbernodes1d: 5
    numberquadraturepoints1d: 5
    numbercomponents: 1
    dimension: 2
  evaluation mode:
    quadratureweights

1 output:
operator field:
  tensor product basis:
    numbernodes1d: 5
    numberquadraturepoints1d: 5
    numbercomponents: 1
    dimension: 2
  evaluation mode:
    gradient
```
"""
function DirichletBDDC(operator::Operator)
    # common constructor
    return BDDC(operator, BDDCInjectionType.harmonic)
end

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


text:
using ParameterizedFunctions, DiffEqBase, SciMLBase
using Test, InteractiveUtils, Latexify

using SpecialFunctions

### ODE Macros

println("Build some examples")
f_t = @ode_def SymCheck begin # Checks for error due to symbol on 1
  dx = x
  dy = -c*y + d*x*y*t^2
end a b c d

@test SciMLBase.__has_syms(f_t)

f_t2 = @ode_def SymCheck2 begin # Checks for error due to symbol on 1
  dx = 1
  dy = -c*y + d*x*y*t^2
end a b c d

f_t3 = @ode_def ExprCheck begin # Checks for error due to symbol on 1
  dx = a*x - b*x*y
  dy = -c*y + d*x*y
end a b c d # Change to π after unicode fix

f = @ode_def_all LotkaVolterra begin
  dx = a*x - b*x*y
  dy = -c*y + d*x*y
end a b c d

f_2 = @ode_def LotkaVolterra3 begin
  dx = a*x - b^2*x*y
  dy = -c*y + d*x*y
end a b c d

println("Test Values")
t = 1.0
u = [2.0,3.0]
p = [1.5,1,3,1]
du = zeros(2)
grad = similar(du)
J = zeros(2,2)
iJ= zeros(2,2)
iW= zeros(2,2)
f(du,u,p,t)
@test du == [-3.0,-3.0]
@test du == f(u,p,t)
f_t(du,u,p,t)
@test du == [2.0,-3.0]
f_t2(du,u,p,t)
@test du == [1.0,-3.0]

println("Test t-gradient")
f.tgrad(grad,u,p,t)
@test grad == zeros(2)
f_t.tgrad(grad,u,p,t)
@test grad == [0.0;12.0]

println("Test Jacobians")
f.jac(J,u,p,t)
@test J  == [-1.5 -2.0
             3.0 -1.0]
@test f.jac(u, p, t) == [-1.5 -2.0; 3.0 -1.0]

@code_llvm SciMLBase.has_jac(f)

println("Test booleans")
@test SciMLBase.has_jac(f) == true

println("Test difficult differentiable")
NJ = @ode_def DiffDiff begin
  dx = a*x - b*x*y
  dy = -c*y + erf(x*y/d)
end a b c d
NJ(du,u,[1.5,1,3,4],t)
@test du == [-3.0;-3*3.0 + erf(2.0*3.0/4)]
@test du == NJ(u, [1.5,1,3,4], t)

test_fail(x,y,d) = erf(x*y/d)
println("Test non-differentiable")
NJ = @ode_def NoJacTest begin
  dx = a*x - μ*x*y
  dy = -c*y + test_fail(x,y,d)
end a μ c d
NJ(du,u,[1.5,1,3,4],t)
@test du == [-3.0;-3*3.0 + erf(2.0*3.0/4)]
@test du == NJ(u,[1.5,1,3,4],t)
@test_throws MethodError NJ.jac(iJ,u,p,t)
# NJ(Val{:jac},t,u,J) # Currently gives E not defined, will be fixed by the next SymEgine

println("Test anonymous definition")

f_t_noname = @ode_def begin # Checks for error due to symbol on 1
  dx = x
  dy = -c*y + d*x*y*t^2
end a b c d

@test SciMLBase.__has_syms(f_t_noname)

f = @ode_def begin
  dx = a*x - b*x*y
  dy = -c*y + d*x*y
end a b c d
@test_nowarn f([0.1,0.2], [1,2], [1,2,3,4], 1)

# Test failures of derivatives should not have #undef
sir_ode = @ode_def SIRModel begin
  dS = -b*S*I
  dI =  b*S*I - g*I
  dR =  g*I
end b g

ext: jl


text:
using Sockets
using OpenSSL
using Nghttp2
using Test

include("testhelpers.jl")

"""
    Items:
[ ] send invalid headers
[ ] add faulty IO stream
[ ] add multiple requests

"""

macro catch_exception_object(code)
    quote
        err = try
            $(esc(code))
            nothing
        catch e
            e
        end
        if err == nothing
            error("Expected exception, got $err.")
        end
        err
    end
end

# Verifies calling into Nghttp library.
@testset "Nghttp2 " begin
    info = nghttp2_version()
    @test unsafe_string(info.proto_str) == "h2"
end

# Verifies creating unsecured Http2 connection
@testset "Http2 connection" begin
    tcp_connection = connect("www.nghttp2.org", 80)

    client_session = Nghttp2.open(tcp_connection)

    # TODO empty request
    io = IOBuffer()
    stream1 = submit_request(
        client_session,
        io,
        [
            ":method" => "GET",
            ":path" => "/",
            ":scheme" => "http",
            ":authority" => "www.nghttp2.org",
            "accept" => "*/*",
        ])

    stream2 = recv(client_session.session)
    stream3 = try_recv(client_session.session)

    @test isnothing(stream3)

    lengths = (length(read_all(stream2)), length(read_all(stream1)))
    @test minimum(lengths) == 6616
    @test maximum(lengths) == 39082

    header_lengths = (length(stream1.headers), length(stream2.headers))
    @test minimum(header_lengths) == 15
    @test maximum(header_lengths) == 18
end

@testset "Http2 greedy client connection" begin
    tcp_connection = connect("www.nghttp2.org", 80)

    client_session = Nghttp2.open(tcp_connection)

    io = IOBuffer()
    stream1 = submit_request(
        client_session,
        io,
        [
            ":method" => "GET",
            ":path" => "/",
            ":scheme" => "http",
            ":authority" => "www.nghttp2.org",
            "accept" => "*/*",
        ])

    err = @catch_exception_object read(stream1, 128 * 1024)
    @test typeof(err) == EOFError
end

@testset "Https2 connection" begin
    tcp_stream = connect("nghttp2.org", 443)

    ssl_ctx = OpenSSL.SSLContext(OpenSSL.TLSv12ClientMethod())
    result = OpenSSL.ssl_set_options(ssl_ctx, OpenSSL.SSL_OP_NO_COMPRESSION | OpenSSL.SSL_OP_NO_TLSv1_2)
    result = OpenSSL.ssl_set_alpn(ssl_ctx, OpenSSL.UPDATE_HTTP2_ALPN)

    ssl_stream = SSLStream(ssl_ctx, tcp_stream, tcp_stream)

    # TODO expose connect
    result = connect(ssl_stream)

    client_session = Nghttp2.open(ssl_stream)

    iob = IOBuffer()
    stream1 = submit_request(
        client_session,
        iob,
        [
            ":method" => "GET",
            ":path" => "/",
            ":scheme" => "https",
            ":authority" => "www.nghttp2.org",
            "accept" => "*/*",
        ])

    stream2 = recv(client_session.session)

    lengths = (length(read_all(stream1)), length(read_all_by_byte(stream2)))
    @test minimum(lengths) == 6616
    @test maximum(lengths) == 39082

    header_lengths = (length(stream1.headers), length(stream2.headers))
    @test minimum(header_lengths) == 16
    @test maximum(header_lengths) == 19
    close(ssl_stream)
end

@testset "Large request" begin
    f1 = @async test_server()
    f2 = @async test_client(IOBuffer(repeat('a', 9 * 65536)))

    fetch(f1)
    @test fetch(f2) == true
end

@testset "Submit request on server session" begin
    # Create an server session and send a request.
    tcp_connection = connect("www.nghttp2.org", 80)

    server_session = Nghttp2.from_accepted(tcp_connection)

    headers = [
        ":method" => "GET",
        ":path" => "/",
        ":scheme" => "http",
        ":authority" => "www.nghttp2.org",
        "accept" => "*/*",
        "user-agent" => "curl/7.75.0"]
    err = @catch_exception_object Nghttp2.send(server_session.session, IOBuffer(), headers)
    @test typeof(err) == Http2ProtocolError
    @test err.lib_error_code == Nghttp2.NGHTTP2_ERR_PROTO
    @test err.msg == "Protocol error"
end

@testset "Invaid request" begin
    f1 = @async test_server()
    f2 = @async test_client(IOBuffer(), INVALID_REQUEST_HEADERS)

    err = @catch_exception_object fetch(f2)
    @test err.task.exception isa Base.IOError

    err = @catch_exception_object fetch(f1)
    @test err.task.exception isa Http2ProtocolError
    @test err.task.exception.lib_error_code == Nghttp2.NGHTTP2_ERR_HTTP_HEADER
end

max_stars_repo_name: grlap/nghttpjl

max_issues_repo_name: grlap/nghttpjl

max_forks_repo_name: grlap/nghttpjl


text:
include("../aimajulia.jl");

using Test;

using Main.aimajulia;

#The following planning tests are from the aima-python doctests

precondition = (map(expr, ["P(x)", "Q(y, z)"]), [expr("Q(x)")]);

effect = ([expr("Q(x)")], [expr("P(x)")]);

plan_action = PlanningAction(expr("A(x, y, z)"), precondition, effect);

arguments = map(expr, ["A", "B", "C"]);

@test (substitute(plan_action, expr("P(x, z, y)"), (arguments...,)) == expr("P(A, C, B)"));

test_planning_kb = FirstOrderLogicKnowledgeBase(map(expr, ["P(A)", "Q(B, C)", "R(D)"]));

@test check_precondition(plan_action, test_planning_kb, (arguments...,));

execute_action(plan_action, test_planning_kb, (arguments...,));

# Found no valid substitutions!
@test (length(ask(test_planning_kb, expr("P(A)"))) == 0);

# Found valid substitutions!
@test (length(ask(test_planning_kb, expr("Q(A)"))) != 0);

# Found valid substitutions!
@test (length(ask(test_planning_kb, expr("Q(A)"))) != 0);

@test (!check_precondition(plan_action, test_planning_kb, (arguments...,)));

air_cargo = air_cargo_pddl();

@test (goal_test(air_cargo) == false);

for action in map(expr, ("Load(C1, P1, SFO)", "Fly(P1, SFO, JFK)", "Unload(C1, P1, JFK)",
                        "Load(C2, P2, JFK)", "Fly(P2, JFK, SFO)", "Unload(C2, P2, SFO)"))
    execute_action(air_cargo, action);
end

@test goal_test(air_cargo);

air_cargo = air_cargo_pddl();

@test (goal_test(air_cargo) == false);

for action in map(expr, ("Load(C2, P2, JFK)", "Fly(P2, JFK, SFO)", "Unload(C2, P2, SFO)",
                        "Load(C1, P1, SFO)", "Fly(P1, SFO, JFK)", "Unload(C1, P1, JFK)"))
	execute_action(air_cargo, action);
end

@test goal_test(air_cargo);

spare_tire = spare_tire_pddl();

@test (goal_test(spare_tire) == false);

for action in map(expr, ("Remove(Flat, Axle)", "Remove(Spare, Trunk)", "PutOn(Spare, Axle)"))
    execute_action(spare_tire, action);
end

@test goal_test(spare_tire);

three_block_tower = three_block_tower_pddl();

@test (goal_test(three_block_tower) == false);

for action in map(expr, ("MoveToTable(C, A)", "Move(B, Table, C)", "Move(A, Table, B)"))
    execute_action(three_block_tower, action);
end

@test goal_test(three_block_tower);

have_cake_and_eat_cake_too = have_cake_and_eat_cake_too_pddl();

@test (goal_test(have_cake_and_eat_cake_too) == false);

for action in map(expr, ("Eat(Cake)", "Bake(Cake)"))
    execute_action(have_cake_and_eat_cake_too, action);
end

@test goal_test(have_cake_and_eat_cake_too);

spare_tire = spare_tire_pddl();

negated_kb = FirstOrderLogicKnowledgeBase([expr("At(Flat, Trunk)")]);

spare_tire_graph = PlanningGraph(spare_tire, negated_kb);

untouched_graph_levels_count = length(spare_tire_graph.levels);

expand_graph(spare_tire_graph);

@test (untouched_graph_levels_count == (length(spare_tire_graph.levels) - 1));

# Apply graphplan() to spare tire planning problem.

spare_tire = spare_tire_pddl();

negated_kb = FirstOrderLogicKnowledgeBase([expr("At(Flat, Trunk)")]);

spare_tire_gp = GraphPlanProblem(spare_tire, negated_kb);

@test (!(typeof(graphplan(spare_tire_gp, ([expr("At(Spare, Axle)"), expr("At(Flat, Ground)")], []))) <: Nothing));

doubles_tennis = doubles_tennis_pddl();

@test (goal_test(doubles_tennis) == false);

for action in map(expr, ["Go(A, LeftBaseLine, RightBaseLine)", "Hit(A, RightBaseLine, Ball)", "Go(A, RightBaseLine, LeftNet)"])
    execute_action(doubles_tennis, action);
end

@test goal_test(doubles_tennis);

# Create dictionary representation of possible refinements for "going to San Francisco airport HLA" (Fig. 11.4).
go_to_sfo_refinements_dict = Dict([Pair("HLA", ["Go(Home,SFO)", "Go(Home,SFO)", "Drive(Home, SFOLongTermParking)", "Shuttle(SFOLongTermParking, SFO)", "Taxi(Home, SFO)"]),
                                    Pair("steps", [["Drive(Home, SFOLongTermParking)", "Shuttle(SFOLongTermParking, SFO)"], ["Taxi(Home, SFO)"], [], [], []]),
                                    Pair("precondition_positive", [["At(Home), Have(Car)"], ["At(Home)"], ["At(Home)", "Have(Car)"], ["At(SFOLongTermParking)"], ["At(Home)"]]),
                                    Pair("precondition_negated", [[], [], [], [], []]),
                                    Pair("effect_add_list", [["At(SFO)"], ["At(SFO)"], ["At(SFOLongTermParking)"], ["At(SFO)"], ["At(SFO)"]]),
                                    Pair("effect_delete_list", [["At(Home)"], ["At(Home)"], ["At(Home)"], ["At(SFOLongTermParking)"], ["At(Home)"]])
                                    ]);

# Base.Test tests for refinements().
function test_refinement_goal_test(kb::FirstOrderLogicKnowledgeBase)
    return ask(kb, expr("At(SFO)"));
end

refinement_lib = Dict([Pair("HLA", ["Go(Home, SFO)", "Taxi(Home, SFO)"]),
                        Pair("steps", [["Taxi(Home, SFO)"], []]),
                        Pair("precondition_positive", [["At(Home)"], ["At(Home)"]]),
                        Pair("precondition_negated", [[],[]]),
                        Pair("effect_add_list", [["At(SFO)"],["At(SFO)"]]),
                        Pair("effect_delete_list", [["At(Home)"], ["At(Home)"]])]);
# Go to San Francisco airport high-level action schema
precondition_positive = Array{Expression, 1}([expr("At(Home)")]);
precondition_negated = Array{Expression, 1}([]);
effect_add_list = Array{Expression, 1}([expr("At(SFO)")]);
effect_delete_list = Array{Expression, 1}([expr("At(Home)")]);
go_sfo = PlanningHighLevelAction(expr("Go(Home, SFO)"),
                                (precondition_positive, precondition_negated),
                                (effect_add_list, effect_delete_list));
# Take Taxi to San Francisco airport high-level action schema
precondition_positive = Array{Expression, 1}([expr("At(Home)")]);
precondition_negated = Array{Expression, 1}([]);
effect_add_list = Array{Expression, 1}([expr("At(SFO)")]);
effect_delete_list = Array{Expression, 1}([expr("At(Home)")]);
taxi_sfo = PlanningHighLevelAction(expr("Go(Home, SFO)"),
                                    (precondition_positive, precondition_negated),
                                    (effect_add_list, effect_delete_list));
go_sfo_pddl = HighLevelPDDL(Array{Expression, 1}([expr("At(Home)")]), [go_sfo, taxi_sfo], test_refinement_goal_test);
result = refinements(go_sfo, go_sfo_pddl, refinement_lib);
@test (length(result) == 1);
@test (result[1].name == "Taxi");
@test (result[1].arguments == (expr("Home"), expr("SFO")));

job_shop_scheduling = job_shop_scheduling_pddl();

@test (goal_test(job_shop_scheduling) == false);

for i in reverse(1:2)
    for j in 1:3
        execute_action(job_shop_scheduling, job_shop_scheduling.jobs[i][j]);
    end
end

@test goal_test(job_shop_scheduling);

max_stars_repo_name: mikhail-j/aima-julia-aimacode-fork

max_issues_repo_name: mikhail-j/aima-julia-aimacode-fork

max_forks_repo_name: mikhail-j/aima-julia-aimacode-fork


text:
workspace()
using Gallium
using HoTRG_lemon.SimulatorModule
using HoTRG_lemon.SpinModule
using HoTRG_lemon.LatticeModule
using TensorOperations
using TensorMatrices_lemon

Gallium.breakpoint(renormalizeSpace!)

dimM = 5
trotterparameter = 0.01
trotteriteration = 10
fieldrange = 0.1

isinginfo = SpinInfo("quantum_ising_2",1.0e-13,0.0)
trotterinfo = TrotterInfo(trotterparameter, trotteriteration)
ising = QuantumIsingModel(isinginfo, trotterinfo)
lattice = buildLattice("quantum_2d_fractal", ising)
simulator = buildSimulator(lattice,dimM)

simulator()


text:
"""
    `CalcVSDiffusionCoefs(boundary, alphas, dt, nu)`

Calculates the RBF coefficients for the diffusion step.

# ARGUMENTS
* `boundary::Boundary`     : Boundary object
* `etas::Array{T}`         : RBF coefficients for rho
* `dt::Real`               : Time step
* `nu::Real`               : Kinematic viscosity
* `sigma = 0.2`            : Gaussian width
where {T<:Real}

# RETURNS
* `betas::Array{Float64}`  : Array of RBF coefficients
"""
function CalcVSDiffusionCoefs(boundary::Boundary,
                              alphas::Array{T},
                              dt::Real,
                              nu::Real;
                              sigma=0.2) where {T<:Real}

    # === Extract geometry from boundary ===
    X = [[point[1], point[2], 0.0] for point in boundary.bodyPTS]
    n_hats = [[n_hat[1], n_hat[2], 0.0] for n_hat in boundary.nHats]
    t_hats = [[t_hat[1], t_hat[2], 0.0] for t_hat in boundary.tHats]
    NPTS = boundary.NPTS_BODY

    # === Constants ===
    CONST1 = -1 / (nu * dt)

    # === Build matrix ===
    A = zeros(NPTS, NPTS)
    for k in 1:NPTS
        for i in 1:NPTS
            if k != i
                R = X[k] - X[i]
                r = norm(R)
                r_hat = R ./ r
                n_hat = n_hats[k]
                A[k,i] = dot(RBF_gauss(r, deriv = 1, sigma=sigma) .* r_hat, n_hat)
            else
                R = [0, 0, 0]
                r = norm(R)
                r_hat = [0, 0, 0]
                n_hat = n_hats[k]
                A[k,i] = dot(RBF_gauss(r, deriv = 1, sigma=sigma) .* r_hat, n_hat)
            end
        end
    end

    # === Build RHS ===
    b = zeros(NPTS)
    for k in 1:NPTS
        b[k] = CalcVS(boundary, alphas, X[k]) * CONST1
    end

    # === Solve linear system ===
    betas = A\b
    return betas

end

"""
    `CalcDiffusion(boundary, beta, X_eval)`

Evaluates the RBF interpolation function: omega(X) = sum(beta_i * phi(X - X_i)).

# ARGUMENTS
* `boundary::Boundary` : Boundary object
* `beta::Array{T}`     : Array containing RBF coefficients
* `X_eval::Array{T}`   : Point to evaluate function at
* `sigma = 0.2`        : Gaussian width
where {T<:Real}

# RETURNS
* `omega::Float64`     : Evaluated function value
"""
function CalcDiffusion(boundary::Boundary,
                       beta::Array{T},
                       X_eval::Array{T};
                       sigma=0.2) where {T<:Real}

    # === Unpack geometry ===
    XP = [[point[1], point[2], 0.0] for point in boundary.bodyPTS]
    NPTS = boundary.NPTS_BODY

    # === Summation over all points for gamma ===
    omega = sum( [beta[i] * RBF_gauss(norm(X_eval - XP[i]), sigma=sigma) for i in 1:NPTS ] )
    return omega
end


text:
@doc raw"""
    sumSigmoidSolverSB(f,z,F,T,Pi)

This function implements Srivastava and Bullo's algorithm to maximize a sum of
general sigoidal functions and returns the maximizer and maximum. `f` is an array
of sigmoidal functions, `df` its derivatives and `z` is the array of inflection points.
"""
function sumSigmoidSolverSB(f,df,z,F,T,Pi,df_inv=nothing)
    ub = maximum([-optimize(x->-df[i](x),F[i],T[i]).minimum for i in 1:length(df)])
    lambdas = collect(range(0.,ub,length=1000))
    FLPs = [FLP_sigmoid(lambdas[i],f,df,df_inv,z,F,T,Pi)[2] for i in 1:length(lambdas)]
    lambdastar = lambdas[partialsortperm(FLPs,1:50,rev=true)]
    lambdas .= range(minimum(lambdastar),maximum(lambdastar),length=1000)
    FLPs .= [FLP_sigmoid(lambdas[i],f,df,df_inv,z,F,T,Pi)[2] for i in 1:length(lambdas)]
    return FLP_sigmoid(lambdas[argmax(FLPs)],f,df,df_inv,z,F,T,Pi)
end

@doc raw"""
    FLP_sigmoid(lambda,f,df,df_inv,z,F,T,Pi)

This function evaluates the linear program relaxation of the maximization of the
sum of sigmoidal functions.
"""
function FLP_sigmoid(lambda,f,df,df_inv,z,F,T,Pi::Real)
    n = length(f)
    Pi_eff = Pi - sum(F)

    df_invs = Vector{Float64}(undef,n)
    if isnothing(df_inv)
        for i in 1:n
            try
                df_invs[i] = find_zero(x-> df[i](x)-lambda,(z[i],T[i]))
            catch
                df_invs[i] = 0.
            end
        end
    else
        df_invs = [df_inv[i](lambda) for i in 1:n]
        df_invs .= [F[i] < df_invs[i] < T[i] ? df_invs[i] : 0.0 for i in 1:n]
    end

    LPindex = Vector{Float64}(undef,n)
    xLP = Vector{Float64}(undef,n)
    gain_rate = Vector{Float64}(undef,3)
    for i in 1:n
        gain_rate .= [df_invs[i] > 0 ? f[i](df_invs[i])/df_invs[i] : 0, F[i]>0 ? f[i](F[i])/F[i] : 0, T[i] > 0 ? f[i](T[i])/T[i] : 0]
        xLPind = argmax(gain_rate)
        xLP[i] = [df_invs[i],F[i],T[i]][xLPind]
        LPindex[i] = gain_rate[xLPind]
    end
    LPperm = sortperm(LPindex,rev = true)
    cum_capacity = 0.0
    k = n
    for i in 1:n
        cum_capacity += xLP[LPperm[i]] - F[LPperm[i]]
        if cum_capacity >= Pi_eff
            k = i
            break
        end
    end
    if k < n
        xLP[LPperm[(k+1):n]] .= F[(k+1):n]
    end
    xLP[LPperm[k]] -= max(cum_capacity - Pi_eff,0.0) # if the capacity is not reached
    return xLP, sum(f[i](xLP[i]) for i in 1:n)
end

max_stars_repo_name: andres-alban/res-alloc-bass-demands

max_issues_repo_name: andres-alban/res-alloc-bass-demands

max_forks_repo_name: andres-alban/res-alloc-bass-demands


text:
# This file was generated, do not modify it. # hide
tm = TunedModel(model=xgb, tuning=Grid(resolution=7),
                resampling=CV(rng=11), ranges=[r1,r2],
                measure=rms)
mtm = machine(tm, X, y)
fit!(mtm, rows=train)

rms(y[test], predict(mtm, rows=test))


text:
module CImGuiExtensions

using CImGui
using CImGui.CSyntax
using CImGui.CSyntax.CStatic
using CImGui.GLFWBackend
using CImGui.OpenGLBackend
using CImGui.GLFWBackend.GLFW
using CImGui.OpenGLBackend.ModernGL
using CImGui: ImVec2, ImVec4, IM_COL32, ImU32
using CSV
using DataFrames
using Printf

include("abstract_types.jl")
include("util.jl")
include("status.jl")
include("path.jl")
include("layout.jl")
include("displayproperties.jl")
include("schemas.jl")
include("image_importer.jl")
include("filedialog.jl")
include("control.jl")
include("plotlines.jl")
include("nestedinterval.jl")
include("labelintervals.jl")
include("truncatedplot.jl")
include("import.jl")
include("export.jl")

        # Abstract types
export  AbstractModel,
        AbstractDialogModel,
        AbstractStatus,
        AbstractDisplayProperties,
        AbstractControl,
        AbstractOperation,
        AbstractImporter,
        AbstractExporter,
        AbstractVendor,
        AbstractProduct,
        AbstractDialogControl,
        AbstractData,
        AbstractLayout,
        AbstractContext,
        AbstractPlotContext,
        AbstractSchema,
        AbstractModelViewControl,
        # Concrete types
        Path,
        FileDialogModel,
        EnabledStatus,
        DisabledStatus,
        ConfirmedStatus,
        UnconfirmedStatus,
        FileDialogDisplayProperties,
        NestedIntervalDisplayProperties,
        FileDialogControl,
        GenericVendor,
        GenericProduct,
        CSVSchema,
        CSVImporter,
        ImageImporter,
        CSVExporter,
        RectangularLayout,
        NestedIntervalControl,
        NestedIntervalContext,
        PlotlinesControl,
        PlotlinesModel,
        PlotlinesDisplayProperties,
        ImportContext,
        ExportContext,
        PlotContext,
        NestedInterval,
        Tickmark,
        Axis,
        TruncatedPlotContext,
        LabelledIntervals,
        LabelledInterval,
        LabelIntervalContext,
        LabelIntervalControl,
        LabelIntervalDisplayProperties,
        IntervalLabels,
        # Functions
        isenabled,
        enable!,
        disable!,
        isrunning,
        get_labelled_intervals,
        load_dataframe,
        get_start,
        get_stop,
        set_start!,
        set_stop!,
        get_width,
        get_height,
        get_padding,
        get_data,
        get_layout,
        stretch_linearly
end # module


text:
using MonteCarloMeasurements: Particles

function test_namedtuple_data(
    idata,
    group,
    names,
    nchains,
    ndraws;
    library = "MyLib",
    coords = Dict(),
    dims = Dict(),
)
    @test idata isa InferenceData
    @test group in ArviZ.groupnames(idata)
    ds = getproperty(idata, group)
    sizes = dimsizes(ds)
    @test length(sizes) == 2 + length(coords)
    vars = vardict(ds)
    for name in string.(names)
        @test name in keys(vars)
        dim = get(dims, name, [])
        s = (x -> length(get(coords, x, []))).(dim)
        @test size(vars[name]) == (nchains, ndraws, s...)
    end
    @test "inference_library" in keys(attributes(ds))
    @test attributes(ds)["inference_library"] == library
    return nothing
end

@testset "from_namedtuple" begin
    rng = MersenneTwister(42)

    nchains, ndraws = 4, 10
    sizes = (x = (), y = (2,), z = (3, 4))
    dims = Dict("y" => ["yx"], "z" => ["zx", "zy"])
    coords = Dict("yx" => ["y1", "y2"], "zx" => 1:3, "zy" => 1:4)

    nts = [
        "NamedTuple" =>
            (; (k => randn(rng, nchains, ndraws, v...) for (k, v) in pairs(sizes))...),
        "Vector{NamedTuple}" => [
            (; (k => randn(rng, ndraws, v...) for (k, v) in pairs(sizes))...)
            for _ in 1:nchains
        ],
        "Matrix{NamedTuple}" => [
            (; (k => randn(rng, v...) for (k, v) in pairs(sizes))...)
            for _ in 1:nchains, _ in 1:ndraws
        ],
        "Vector{Vector{NamedTuple}}" => [
            [(; (k => randn(rng, v...) for (k, v) in pairs(sizes))...) for _ in 1:ndraws] for _ in 1:nchains
        ],
        "Vector{NamedTuple} particles" => [
            (; (k => Particles(randn(rng, ndraws, v...)) for (k, v) in pairs(sizes))...) for _ in 1:nchains
        ],
    ]

    @testset "posterior::$(type)" for (type, nt) in nts
        idata1 = from_namedtuple(nt; dims = dims, coords = coords, library = "MyLib")
        test_namedtuple_data(
            idata1,
            :posterior,
            keys(sizes),
            nchains,
            ndraws;
            library = "MyLib",
            coords = coords,
            dims = dims,
        )

        idata2 =
            convert_to_inference_data(nt; dims = dims, coords = coords, library = "MyLib")
        test_namedtuple_data(
            idata2,
            :posterior,
            keys(sizes),
            nchains,
            ndraws;
            library = "MyLib",
            coords = coords,
            dims = dims,
        )
    end

    @testset "$(group)" for group in [
        :posterior_predictive,
        :sample_stats,
        :predictions,
        :log_likelihood,
    ]
        @testset "::$(type)" for (type, nt) in nts
            idata1 = convert_to_inference_data(
                nt;
                (group => nt,)...,
                dims = dims,
                coords = coords,
                library = "MyLib",
            )
            test_namedtuple_data(
                idata1,
                :posterior,
                keys(sizes),
                nchains,
                ndraws;
                library = "MyLib",
                coords = coords,
                dims = dims,
            )
            test_namedtuple_data(
                idata1,
                group,
                keys(sizes),
                nchains,
                ndraws;
                library = "MyLib",
                coords = coords,
                dims = dims,
            )

            idata2 = convert_to_inference_data(
                nt;
                (group => keys(sizes),)...,
                dims = dims,
                coords = coords,
                library = "MyLib",
            )
            test_namedtuple_data(
                idata2,
                group,
                keys(sizes),
                nchains,
                ndraws;
                library = "MyLib",
                coords = coords,
                dims = dims,
            )
        end
    end

    @testset "$(group)" for group in [:prior_predictive, :sample_stats_prior]
        @testset "::$(type)" for (type, nt) in nts
            idata1 = from_namedtuple(;
                prior = nt,
                (group => nt,)...,
                dims = dims,
                coords = coords,
                library = "MyLib",
            )
            test_namedtuple_data(
                idata1,
                :prior,
                keys(sizes),
                nchains,
                ndraws;
                library = "MyLib",
                coords = coords,
                dims = dims,
            )
            test_namedtuple_data(
                idata1,
                group,
                keys(sizes),
                nchains,
                ndraws;
                library = "MyLib",
                coords = coords,
                dims = dims,
            )

            idata2 = from_namedtuple(;
                prior = nt,
                (group => keys(sizes),)...,
                dims = dims,
                coords = coords,
                library = "MyLib",
            )
            test_namedtuple_data(
                idata2,
                group,
                keys(sizes),
                nchains,
                ndraws;
                library = "MyLib",
                coords = coords,
                dims = dims,
            )
        end
    end

    @testset "$(group)" for group in
                            [:observed_data, :constant_data, :predictions_constant_data]
        _, nt = nts[1]
        idata = from_namedtuple(
            nt;
            (group => Dict("w" => [1.0, 2.0]),)...,
            dims = Dict("w" => ["wx"]),
            coords = Dict("wx" => 1:2),
            library = "MyLib",
        )
        @test idata isa InferenceData
        @test group in ArviZ.groupnames(idata)
        ds = getproperty(idata, group)
        sizes = dimsizes(ds)
        @test length(sizes) == 1
        @test "w" in keys(vardict(ds))
        @test "inference_library" in keys(attributes(ds))
        @test attributes(ds)["inference_library"] == "MyLib"
    end
end


text:
""" Categories of C-sets and attributed C-sets.
"""
module CSets
export ACSetTransformation, CSetTransformation, components, force, is_natural

using Compat: isnothing

using AutoHashEquals
using Reexport
using StaticArrays: SVector

@reexport using ...CSetDataStructures
using ...GAT, ..FreeDiagrams, ..Limits, ..FinSets
import ..Limits: limit, colimit, universal
import ..FinSets: force
using ...Theories: Category, CatDesc, AttrDesc
import ...Theories: dom, codom, compose, ⋅, id

# C-set transformations
#######################

""" Transformation between attributed C-sets.

A morphism of C-sets is a natural transformation: a transformation between
functors C -> Set satisfying the naturality axiom for all morphisms in C. This
struct records the data of a transformation; it does not enforce naturality.

The transformation has a component for every object in C. When C-sets have
attributes, the data types are assumed to be fixed. Thus, the naturality axiom
for data attributes is a commutative triangle, rather than a commutative square.
"""
@auto_hash_equals struct ACSetTransformation{CD <: CatDesc, AD <: AttrDesc{CD},
    Comp <: NamedTuple, Dom <: AbstractACSet{CD,AD}}
  components::Comp
  dom::Dom
  codom::Dom

  function ACSetTransformation{CD,AD}(components::NamedTuple, X::Dom, Y::Dom) where
      {Ob, CD <: CatDesc{Ob}, AD <: AttrDesc{CD}, Dom <: AbstractACSet{CD,AD}}
    @assert keys(components) ⊆ Ob
    coerced_components = NamedTuple{Ob}(
      coerce_component(ob, get(components, ob) do; Int[] end, X, Y)
      for ob in Ob)
    new{CD,AD,typeof(coerced_components),Dom}(coerced_components, X, Y)
  end
end

function coerce_component(ob::Symbol, f::FinFunction{Int}, X, Y)
  @assert length(dom(f)) == nparts(X,ob) "Domain error in component $ob"
  @assert length(codom(f)) == nparts(Y,ob) "Codomain error in component $ob"
  return f
end
function coerce_component(ob::Symbol, f, X, Y)::FinFunction{Int}
  FinFunction(f, nparts(X,ob), nparts(Y,ob))
end

ACSetTransformation(components::NamedTuple, X::Dom, Y::Dom) where
    {CD <: CatDesc, AD <: AttrDesc{CD}, Dom <: AbstractACSet{CD,AD}} =
  ACSetTransformation{CD,AD}(components, X, Y)
ACSetTransformation(X::Dom, Y::Dom; components...) where
    {CD <: CatDesc, AD <: AttrDesc{CD}, Dom <: AbstractACSet{CD,AD}} =
  ACSetTransformation{CD,AD}((; components...), X, Y)

const CSetTransformation{CD, Comp, Dom <: AbstractCSet{CD}} =
  ACSetTransformation{CD,AttrDesc{CD,(),(),(),()},Comp,Dom}

CSetTransformation(components, X::Dom, Y::Dom) where Dom <: AbstractCSet =
  ACSetTransformation(components, X, Y)
CSetTransformation(X::Dom, Y::Dom; components...) where Dom <: AbstractCSet =
  ACSetTransformation(X, Y; components...)

components(α::ACSetTransformation) = α.components
Base.getindex(α::ACSetTransformation, ob) = α.components[ob]

""" Is the transformation between C-sets a natural transformation?

Uses the fact that to check whether a transformation is natural, it suffices to
check the naturality equation on a generating set of morphisms.
"""
function is_natural(α::ACSetTransformation{CD,AD}) where {CD,AD}
  X, Y = dom(α), codom(α)
  for (f, c, d) in zip(CD.hom, CD.dom, CD.codom)
    Xf, Yf, α_c, α_d = subpart(X,f), subpart(Y,f), α[c], α[d]
    all(Yf[α_c(i)] == α_d(Xf[i]) for i in eachindex(Xf)) || return false
  end
  for (f, c) in zip(AD.attr, AD.adom)
    Xf, Yf, α_c = subpart(X,f), subpart(Y,f), α[c]
    all(Yf[α_c(i)] == Xf[i] for i in eachindex(Xf)) || return false
  end
  return true
end

force(α::ACSetTransformation) =
  ACSetTransformation(map(force, components(α)), dom(α), codom(α))

# Category of C-sets
####################

@instance Category{ACSet, ACSetTransformation} begin
  dom(α::ACSetTransformation) = α.dom
  codom(α::ACSetTransformation) = α.codom

  id(X::ACSet) = ACSetTransformation(map(id, finsets(X)), X, X)

  function compose(α::ACSetTransformation, β::ACSetTransformation)
    # Question: Should we incur cost of checking that codom(β) == dom(α)?
    ACSetTransformation(map(compose, components(α), components(β)),
                        dom(α), codom(β))
  end
end

finsets(X::ACSet) = map(table -> FinSet(length(table)), tables(X))

# Limits and colimits
#####################

""" Limit of C-sets that stores the pointwise limits in FinSet.
"""
struct CSetLimit{Ob <: AbstractCSet, Diagram, Cone <: Multispan{Ob},
                 Limits <: NamedTuple} <: AbstractLimit{Ob,Diagram}
  diagram::Diagram
  cone::Cone
  limits::Limits
end

""" Colimit of attributed C-sets that stores the pointwise colimits in FinSet.
"""
struct ACSetColimit{Ob <: AbstractACSet, Diagram, Cocone <: Multicospan{Ob},
                    Colimits <: NamedTuple} <: AbstractColimit{Ob,Diagram}
  diagram::Diagram
  cocone::Cocone
  colimits::Colimits
end

# Compute limits and colimits of C-sets by reducing to those in FinSet using the
# "pointwise" formula for (co)limits in functor categories.

function limit(diagram::AbstractFreeDiagram{ACS}) where
    {CD <: CatDesc, ACS <: AbstractCSet{CD}}
  limits = map(limit, unpack_diagram(diagram))
  Xs = cone_objects(diagram)
  Y = ACS()
  for (c, lim) in pairs(limits)
    add_parts!(Y, c, length(ob(lim)))
  end
  for (f, c, d) in zip(CD.hom, CD.dom, CD.codom)
    Yfs = map(legs(limits[c]), Xs) do π, X
      compose(π, FinFunction(subpart(X, f), nparts(X, d)))
    end
    Yf = universal(limits[d], Multispan(ob(limits[c]), Yfs))
    set_subpart!(Y, f, collect(Yf))
  end
  πs = pack_components(map(legs, limits), map(X -> Y, Xs), Xs)
  CSetLimit(diagram, Multispan(Y, πs), limits)
end

function universal(lim::CSetLimit, cone::Multispan)
  components = map(universal, lim.limits, unpack_diagram(cone))
  CSetTransformation(components, apex(cone), ob(lim))
end

function colimit(diagram::AbstractFreeDiagram{ACS}) where
    {CD <: CatDesc, AD <: AttrDesc{CD}, Ts, ACS <: AbstractACSet{CD,AD,Ts}}
  # Colimit of C-set without attributes.
  colimits = map(colimit, unpack_diagram(diagram))
  Xs = cocone_objects(diagram)
  Y = ACS()
  for (c, colim) in pairs(colimits)
    add_parts!(Y, c, length(ob(colim)))
  end
  for (f, c, d) in zip(CD.hom, CD.dom, CD.codom)
    Yfs = map(legs(colimits[d]), Xs) do ι, X
      compose(FinFunction(subpart(X, f), nparts(X, d)), ι)
    end
    Yf = universal(colimits[c], Multicospan(ob(colimits[d]), Yfs))
    set_subpart!(Y, f, collect(Yf))
  end
  ιs = pack_components(map(legs, colimits), Xs, map(X -> Y, Xs))

  # Set data attributes by canonical inclusion from attributes in diagram.
  for (attr, c, d) in zip(AD.attr, AD.adom, AD.acodom)
    T = Ts.parameters[d]
    data = Vector{Union{Some{T},Nothing}}(nothing, nparts(Y, c))
    for (ι, X) in zip(ιs, Xs)
      for i in parts(X, c)
        j = ι[c](i)
        if isnothing(data[j])
          data[j] = Some(subpart(X, i, attr))
        else
          val1, val2 = subpart(X, i, attr), something(data[j])
          val1 == val2 || error(
            "ACSet colimit does not exist: $attr attributes $val1 != $val2")
        end
      end
    end
    set_subpart!(Y, attr, map(something, data))
  end

  ACSetColimit(diagram, Multicospan(Y, ιs), colimits)
end

function universal(colim::ACSetColimit, cocone::Multicospan)
  components = map(universal, colim.colimits, unpack_diagram(cocone))
  ACSetTransformation(components, ob(colim), apex(cocone))
end

""" Diagram in C-Set → named tuple of diagrams in FinSet
"""
unpack_diagram(diagram::DiscreteDiagram{<:AbstractACSet}) =
  map(DiscreteDiagram, unpack_finsets(ob(diagram)))
unpack_diagram(span::Multispan{<:AbstractACSet}) =
  map(Multispan, finsets(apex(span)), unpack_components(legs(span)))
unpack_diagram(cospan::Multicospan{<:AbstractACSet}) =
  map(Multicospan, finsets(apex(cospan)), unpack_components(legs(cospan)))
unpack_diagram(para::ParallelMorphisms{<:AbstractACSet}) =
  map(ParallelMorphisms, unpack_components(hom(para)))

""" Vector of C-sets → named tuple of vectors of FinSets
"""
unpack_finsets(Xs::AbstractVector{<:AbstractACSet{CD}}) where
    {Ob, CD <: CatDesc{Ob}} =
  NamedTuple{Ob}([ map(X -> FinSet{Int,Int}(nparts(X, ob)), Xs) for ob in Ob ])

""" Vector of C-set transformations → named tuple of vectors of FinFunctions
"""
unpack_components(αs::AbstractVector{<:ACSetTransformation{CD}}) where
    {Ob, CD <: CatDesc{Ob}} =
  NamedTuple{Ob}([ map(α -> α[ob], αs) for ob in Ob ])

""" Named tuple of vectors of FinFunctions → vector of C-set transformations
"""
function pack_components(fs::NamedTuple{Ob}, doms, codoms) where Ob
  components = map((x...) -> NamedTuple{Ob}(x), fs...) # XXX: Is there a better way?
  map(ACSetTransformation, components, doms, codoms)
end

""" Objects in diagram that will have explicit legs in limit cone.

Encodes common conventions such as, when taking a pullback of a cospan, not
explicitly including a cone leg for the cospan apex since it can be computed
from the other legs.
"""
cone_objects(diagram) = ob(diagram)
cone_objects(cospan::Multicospan) = feet(cospan)
cone_objects(para::ParallelMorphisms) = SVector(dom(para))

""" Objects in diagram that will have explicit legs in colimit cocone.
"""
cocone_objects(diagram) = ob(diagram)
cocone_objects(span::Multispan) = feet(span)
cocone_objects(para::ParallelMorphisms) = SVector(codom(para))

end


text:
# # Эксперимент №1 Не зашумленные

# ## Набор параметров
α       = 0.004;        # Параметр регуляризации
w       = 0.0001;       # Эмпирический параметр регуляризации
S       =  5000;        # Количество итераций
β       = 0.001;        # Шаг минимизации
#
x_tp    = 0.05;         # Стартовое местоположение фронта
T_end   = 0.47;         # Регулируем конечное местоположение фронта
ε       = 10^(-1.5);    # Крутизна фронта
Nx      = 15000;          # Число интервалов по ``X``
Mt      = 30000;         # Число интервалов по ``T``

using Dates;
timestamp           = String(Dates.format(now(), "yyyy-mm-dd_HH-MM"));
wdir                = String("Expirement1___" * timestamp);
prefix              = String("Exp1_" * timestamp * "_");
n_direct_draft      = prefix * "direct_draft.png";
n_direct_gif        = prefix * "direct.gif";
n_data              = prefix * "data.png";
n_velocity          = prefix * "velocity.png";
n_initial_guess     = prefix * "initial_guess.png";
n_Js_min_data_guess = prefix * "Js_withguess.jld";
n_Qs_min_data_guess = prefix * "Qs_withguess.jld";
n_Js_min_data       = prefix * "Js_noguess.jld";
n_Qs_min_data       = prefix * "Qs_noguess.jld";
n_guess_draft       = prefix * "minimization_draft_withguess.png"
n_noguess_draft     = prefix * "minimization_draft_noguess.png"
n_minim_gif_guess   = prefix * "minimization_withguess.gif";
n_minim_gif_noguess = prefix * "minimization_noguess.gif";
n_minim_avi_guess   = prefix * "minimization_withguess.avi";
n_minim_avi_noguess = prefix * "minimization_noguess.avi";
n_minim_mp4_guess   = prefix * "minimization_withguess.mp4";
n_minim_mp4_noguess = prefix * "minimization_noguess.mp4";

using Markdown
info = """Вы зпустили эксперимент №1. Сценарий скрипта:
  - Создадим каталог `$(wdir)`.
  - Перейдем в него и все результаты поместим там.
  - Все файлы будем сохранять с префиксом `$(prefix)`
  - Сгенерируем экспериментальные данные на сетке `$(Nx+1) × $(Mt+1)`
    * Сохраним эскиз прямого решения `$(n_direct_draft)`
    * Сохраним анимацию решения прямой задачи, с помощью которой генерировали
        экспериментальные даные `$(n_direct_gif)`.
  - Найдем начальное приближение с параметром регуляризации `$(α)`
    * Сохраним начальное приближение `$(n_initial_guess)`
  - Начнем процесс минимизации используя найденное приближение
    * Сериализуем данные и сохраним `$(n_Js_min_data_guess)`, `$(n_Qs_min_data_guess)`
    * Сохраним Черновик процесса минимизации `$(n_guess_draft)`
    * Сохраним Анимацию процесса минимизации `$(n_minim_gif_guess)`
    * Сохраним mp4 Анимации процесса минимизации `$(n_minim_mp4_guess)`
    * Сохраним mp4 Анимации процесса минимизации `$(n_minim_avi_guess)`
  - Начнем процесс минимизации с нулевого начального приближения
    * Сериализуем данные и сохраним `$(n_Js_min_data)`, `$(n_Qs_min_data)`
    * Сохраним Черновик процесса минимизации `$(n_noguess_draft)`
    * Сохраним Анимацию процесса минимизации `$(n_minim_gif_noguess)`
    * Сохраним mp4 Анимации процесса минимизации `$(n_minim_mp4_noguess)`
    * Сохраним mp4 Анимации процесса минимизации `$(n_minim_avi_noguess)`
"""
display(Markdown.parse(info))

@info "Переходим в $(wdir)"
run(`mkdir $(wdir)`)
run(`cd $(wdir)`)

# -----------------------------------------------------------------------------
using NonLinearReactionAdvectionDiffusionWithFrontData
using NonLinearReactionAdvectionDiffusionWithFrontData: heterogeneity_map;
using Serialization;
using Plots; pyplot();
using Dierckx;
# -----------------------------------------------------------------------------

# ## Решение на точных данных

# ### Решение прямой задачи для генерирования экспериментальной информации
@info "Решение прямой задачи для генерирования экспериментальной информации"
# -----------------------------------------------------------------------------
a, b, t₀, T, N, M, ε, Xₙ, Tₘ, qₙ, ulₘ, urₘ, u₀ = dparams(x_tp = x_tp,
                                                         ε = ε,
                                                         Nx = Nx,
                                                         Mt = Mt,
                                                         T_end = T_end);
u, XX, TP = solve(u₀, Xₙ, N, Tₘ, M, ε, ulₘ, urₘ, qₙ);
ϕl, ϕr, ϕ, f1_data, f2_data = generate_obs_data(u, Xₙ, N, Tₘ, M, qₙ, ulₘ, urₘ);
directP = draft(u, Xₙ, N, Tₘ, M, title = "Эскиз прямого решения")
savefig(directP, n_direct_draft);
# -----------------------------------------------------------------------------

# ### Найдем начальное приближение
@info "Найдем начальное приближение"
# -----------------------------------------------------------------------------
q_guess = initial_guess(f1_data, Xₙ, N, Tₘ, M, ulₘ, urₘ, 0.005);
initial_guessP = plot(Xₙ, qₙ, xlabel = "X", ylabel = "q(x)", label="Истинное");
initial_guessP = plot!(Xₙ, q_guess, label="Найденное")
savefig(n_initial_guess);
# -----------------------------------------------------------------------------

# -----------------------------------------------------------------------------
# В следующей строке, мы записываем все параметры в latex формате,
# не забываем экранировать все спецсимволы.
# Дальше, мы отобразим все эти параметры на графике.
using Printf;
params = "\$w = $(w), \\varepsilon = $(@sprintf("%.4f", ε)), "*
"\\beta = $(β), f_1 \\in $(@sprintf("[%.2f, %.2f]", extrema(f1_data)...)), "*
"\\alpha = $(α)\$"
nothing #hide
# -----------------------------------------------------------------------------

# -----------------------------------------------------------------------------
spl = Spline1D(Xₙ, q_guess);
Nx      = 500;          # Число интервалов по ``X``
Mt      = 1000;         # Число интервалов по ``T``
a, b, t₀, T, N, M, ε, Xₙ, Tₘ, qₙ, ulₘ, urₘ, u₀ = dparams(x_tp = x_tp,
                                                         ε = ε,
                                                         Nx = Nx,
                                                         Mt = Mt,
                                                         T_end = T_end);
u, XX, TP = solve(u₀, Xₙ, N, Tₘ, M, ε, ulₘ, urₘ, qₙ);
ϕl, ϕr, ϕ, f1_data, f2_data = generate_obs_data(u, Xₙ, N, Tₘ, M, qₙ, ulₘ, urₘ);
q₀ = spl(Xₙ);

# ### Старт с найденного приближения на точных данных
# -----------------------------------------------------------------------------
@time qs, guess_Js, guess_Qs = minimize(q₀, u₀, ulₘ, urₘ, Xₙ, N, Tₘ, M, ε, f1_data, f2_data,
                            S = S, β = β, w = w, showProgress = true)
serialize(n_Js_min_data_guess, guess_Js);
serialize(n_Qs_min_data_guess, guess_Qs);
# -----------------------------------------------------------------------------

# ### Отрисовка результатов минимизации
# -----------------------------------------------------------------------------
a, b, c = minimization_draft(qₙ, guess_Qs, Xₙ, N, guess_Js, zoom = true,
                              annotate_string = params, zoom_chunk = 17//18)
plot(a, size = (800, 800))
plot(b)
plot(c)
withguessP = plot(a, b, c, layout = (3,1), size = (1000, 2400));
savefig(withguessP, n_guess_draft)
nothing; #hide
# -----------------------------------------------------------------------------
frames = collect(1:div(S, 50):S);
make_minimzation_gif(
    guess_Js, guess_Qs, qₙ, Xₙ,
    name = n_minim_gif_guess,
    frames_to_write = frames,
    β = β
)
run(`gif2mp4 $(n_minim_gif_guess) $(n_minim_mp4_guess)`)
run(`ffmpeg -i  $(n_minim_gif_guess) $(n_minim_avi_guess) -vcodec copy -acodec copy`)
# -----------------------------------------------------------------------------

# -----------------------------------------------------------------------------
# ------------------ Старт с нулевого приближения -----------------------------
# -----------------------------------------------------------------------------
q₀ = zero(Xₙ)
# -----------------------------------------------------------------------------
@time qs, noguess_Js, noguess_Qs = minimize(q₀, u₀, ulₘ, urₘ, Xₙ, N, Tₘ, M, ε, f1_data, f2_data,
                            S = S, β = β, w = w, showProgress = true)

serialize(n_Js_min_data, noguess_Js);
serialize(n_Qs_min_data, noguess_Qs);
# -----------------------------------------------------------------------------

# ### Эскиз процесса минимизации
# -----------------------------------------------------------------------------
a, b, c = minimization_draft(qₙ, noguess_Qs, Xₙ, N, noguess_Js, zoom = true,
                              annotate_string = params, zoom_chunk = 17//18)
plot(a, size = (800, 800))
plot(b)
plot(c)
withguessP = plot(a, b, c, layout = (3,1), size = (1000, 2400));
savefig(withguessP, n_noguess_draft)
nothing; #hide
frames = collect(1:div(S, 50):S);
make_minimzation_gif(
    noguess_Js, noguess_Qs, qₙ, Xₙ,
    name = n_minim_gif_noguess,
    frames_to_write = frames,
    β = β
)
run(`gif2mp4 $(n_minim_gif_noguess) $(n_minim_mp4_noguess)`)
run(`ffmpeg -i  $(n_minim_gif_noguess) $(n_minim_avi_noguess) -vcodec copy -acodec copy`)
# -----------------------------------------------------------------------------


text:
@inline function generate_compat_entry(v::VersionNumber)::String
    if v.major == 0 && v.minor == 0
        return "0.0.$(v.patch)"
    else
        return "$(v.major).$(v.minor)"
    end
end

# The `_remove_trailing_zeros` function is based on
# the `compat_version` function from:
# https://github.com/invenia/PkgTemplates.jl/blob/master/src/plugins/project_file.jl

@inline function _remove_trailing_zeros(v::VersionNumber)::String
    if v.patch == 0
        if v.minor == 0
            if v.major == 0 # v.major is 0, v.minor is 0, v.patch is 0
                throw(DomainError("0.0.0 is not a valid input"))
            else # v.major is nonzero and v.minor is 0 and v.patch is 0
                return "$(v.major)"
            end
        else # v.minor is nonzero, v.patch is 0
            return "$(v.major).$(v.minor)"
        end
    else # v.patch is nonzero
        return "$(v.major).$(v.minor).$(v.patch)"
    end
end

# @inline function MajorMinorVersion(version::VersionNumber)
#     major = version.major::Base.VInt
#     minor = version.minor::Base.VInt
#     result = MajorMinorVersion(major, minor)
#     return result
# end

# @inline MajorMinorVersion(::Nothing) = nothing

# @inline function VersionNumber(x::MajorMinorVersion)
#     major = x.major::Base.VInt
#     minor = x.minor::Base.VInt
#     result = VersionNumber(x, y)::VersionNumber
#     return result
# end

# @inline function Base.isless(a::MajorMinorVersion, b::MajorMinorVersion)
#     a_versionnumber = VersionNumber(a)::VersionNumber
#     b_versionnumber = VersionNumber(b)::VersionNumber
#     result = Base.isless(a_versionnumber, b_versionnumber)
#     return result
# end


text:
"""
    WebMSISE00

Load this module to access the functions below...

"""
module WebMSISE00

export ReadMSISE
export WriteMSISEData
export ReadMSISEData

using HTTP, Gumbo, DataStructures
using DataFrames
using CSV
using JLD2
using AbstractTrees

"""
    ReadMSISE(case_inputs::Dict, casename::String = "")

    Executes an MSISE00 run on NASA's Model Web site.

    Uses HTTP.jl and Gumbo.jl to access the appropriate PHP form on the NASA
    Model Web site, submits the form, then scrapes the resulting data. The data
    is returned in a DataFrame.

    INPUT:
        `case_inputs` is a Dict that gets merged with the default settings
            contained in the payload Dict below.

        `casename` is used as the basis of the filename for storing output into
            a JLD2 file (no string means "don't write a file")

    OUTPUT: casename.jld2 (optional)

    RETURNS: DataFrame containing the requested atmospheric data
"""
function ReadMSISE(case_inputs::Dict, casename::String = "")

    url = "https://ccmc.gsfc.nasa.gov/modelweb/models/nrlmsise00.php"

    resp = HTTP.get(url) # grab the entire page

# The hierarchical structure of the HTML body is broken into a many
#   dimensional array. The trick is to page thru until you find what you want.

# The following commented code would be needed in case we wish to autodetect
#   the form data and do something with it

# parsed = parsehtml(String(resp.body))  # parse it into a ginormous array
# body = parsed.root[2] # the body (as opposed to the header)
# form = body[4][1][6]

# TIP: Chrome developer mode is useful for looking at the form to determine
#   the variables and values to POST.
#
# the form items to POST start with a dict of the responses of the format:
#   key(name) => value

    url2 = "https://ccmc.gsfc.nasa.gov/cgi-bin/modelweb/models/vitmo_model.cgi"

# NOTE: everything is submitted to the form as a string
# NOTE2: using OrderedDict in case order actually matters
# NOTE3: multiple selected checkboxes in the form result in needing to
#   use multiple values with a single key... obviously a problem when using
#   anything related to a Dict. **SOLVED** by storing checkbox items in an
#   array, then parsing that array using an inner loop

# vars Array values
#   Independent variables:
#   00 = Year; 01 = Month; 02 = Day of Month; 03 = Day of Year;
#   04 = Hour of Day; 05 = Height (km); 06 = Latitude; 07 = Longitude
#   Calculated Parameters:
#   08 = Atomic Oxygen number density (cm^-3); 09 = N2 number density;
#   10 = O2 number density; 11 = total mass density (g/cm^3);
#   12 = Neutral Temp (K); 13 = Exospheric Temp (K); 14 = He number density;
#   15 = Ar number density; 16 = H number density; 17 = N number density;
#   18 = Anomalous Oxygen number density
#   Indices taken from user input vs. data base
#   19 = F10.7 (daily); 20 = F10.7 (3 month avg); 21 = ap (daily);
#       (remainder won't generally be used)

    payload = OrderedDict{String, Any}(
        "model" => "nrlmsise",
        "year" => "2020",
        "month" => "01",
        "day" => "01",
        "time_flag" => "0",
        "hour" => "1.5",
        "geo_flag" => "0.",
        "latitude" => "55.",
        "longitude" => "45.",
        "height" => "100.",
        "profile" => "1",
        "start" => "0.",
        "stop" => "1000.",
        "step" => "50.",
        "f10_7" => "",
        "f10_7_3" => "",
        "ap" => "",
        "format" => "2",
        "vars" => ["05", "08", "09"],
        "linestyle" => "solid",
        "charsize" => "1.0",
        "symbol" => "2",
        "symsize" => "1.0",
        "yscale" => "Lin",
        "xscale" => "Lin",
        "imagex" => "640",
        "imagey" => "480"
    )

# Merge the input with the payload Dict

    merge!(payload, case_inputs)

# Convert the Dict into the format expected...
#   That format was discovered using Chrome developer tool, looking at the outgoing
#   form submission

    form_data = ""

    for (k, v) in payload
        if isa(v, Array) || isa(v, Tuple)
            for w in v
                form_data = form_data * k * "=" * w * "&"
            end
        else
            form_data = form_data * k * "=" * v * "&"
        end
    end
    form_data = rstrip(form_data, ['&']) # get rid of extraneous trailing '&'

    resp2 = HTTP.request("POST", url2,
        ["Content-Type" => "application/x-www-form-urlencoded"], form_data)

    body2 = parsehtml(String(resp2.body)).root[2]

# find all HREFs in body2 (there are many)
# body2 is a Tree, which is why we needed AbstractTrees

    links = String[]
    for element in PreOrderDFS(body2)
        if element isa HTMLElement
           if "href" in keys(attrs(element))
               push!(links, attrs(element)["href"])
           end
        end
    end

# extract the one that matches what we need
    link = ""
    url3 = ""
    for link in links
        if startswith(link, "http") && endswith(link,".lst")
            url3 = link
            break
        end
    end

    resp3 = HTTP.get(url3)

    result = parsehtml(String(resp3.body)).root[2][1].text

# Next step: parse the actual data from 'result'
# result is pure ASCII. It contains an echo of the important inputs,
#   a 'decoder ring' for the columns, and then the actual data (in columns
#   labeled 1, 2, ..., hence the need for the 'decoder ring')
# Optional step: Write result out to a file for safe keeping?

    model_input_echo = String[]
    cols_key = String[]
    data = String[]
    i = 0
    df1 = DataFrame()

# parse line by line, breaking result into 3 pieces

    for line in eachline(IOBuffer(result))
        if !occursin("Selected output parameters", line) && i==0
            push!(model_input_echo, line)
        elseif occursin("Selected output parameters", line)
            i = i + 1
            continue
        elseif i == 1 && occursin(",", line)
            push!(cols_key, line)
        else
            push!(data, line)
        end
    end

# Parse column key to a DataFrame so we can pull out the column names

    for l in cols_key
        thisline = CSV.read(IOBuffer(l), DataFrame; header = 0, delim = "  ")
        df1 = vcat(df1, thisline)
    end

    cols_names = df1.Column2

# need to reassemble data into a long string so we can read into DataFrame
#   without trying to parse line by line (as in the above loop)
#   NOTE: start on line 3

    data2 = ""
    for l in data[3:length(data)]
        data2 = data2 * l * "\n"
    end

# Seems inefficient to read the data twice (probably can replace last
#   push! in the line-by-line parsing with what's immediately above)

    df2 = CSV.read(IOBuffer(data2), DataFrame; header = 0, delim = " ",
        ignorerepeated = true)

# Now we need to insert the column names...

    rename!(df2,cols_names)

    if length(casename) > 0
# Convert model_input_echo from a vector of one line strings
        input_echo = ""
        for l in model_input_echo
            input_echo = input_echo * l * "\n"
        end

        WriteMSISEData(casename, input_echo, df2)
    end

# the DataFrame is now ready for use... so return it
    return df2

end # ReadMSISE

# RunMSISEUI: a simple user interface to set up a case to grab from ModelWeb

"""
    WriteMSISEData(casename::String, run_echo::String, rundata::DataFrame)

Writes the MSISE00 results to a JLD2 file.

Used for storing chunks of the MSISE00 output to a file for later use.

INPUT:
    `casename` = basis for the filename ('.jld2' extension is added)
    `run_echo` = the ModelWeb site echos back the inputs at the top of its
        output, this string contains just that portion of the output
    `rundata` = the DataFrame containing the parsed MSISE00 output

OUTPUT: 'casename.jld2' file containing stored data

RETURNS: none
"""
function WriteMSISEData(casename::String, run_echo::String, rundata::DataFrame)
    filename = casename * ".jld2"

    @save filename run_echo rundata
end

"""
    ReadMSISEData(casename::String)

Reads MSISE00 results from a JLD2 file.

INPUT:
    `casename` = basis for the filename ('.jld2' extension is added)

OUTPUT: none

RETURNS:
    `run_echo` = the ModelWeb site echos back the inputs at the top of its
        output, this string contains just that portion of the output
    `rundata` = the DataFrame containing the parsed MSISE00 output
"""
function ReadMSISEData(casename::String)
    filename = casename * ".jld2"

    @load filename run_echo rundata

    return run_echo, rundata
end

end # module


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "MathGL"
version = v"2.4.4"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://sourceforge.net/projects/mathgl/files/mathgl/mathgl%202.4.4/mathgl-2.4.4.tar.gz", "0e5977196635962903eaff9b2f759e5b89108339b6e71427036c92bfaf3149e9"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
apk add g++
cd $WORKSPACE/srcdir
for f in ${WORKSPACE}/srcdir/patches/*.patch; do
    atomic_patch -p1 ${f}
done
cd mathgl-2.4.4/
mkdir build
cd build
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release -Denable-all-docs=OFF -Denable-qt5=ON -Denable-qt4=OFF -Denable-mpi=OFF  ../.
make -j
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "windows"),
    Platform("x86_64", "linux"; libc="glibc"),
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libmgl-qt5", :mgl_qt5),
    LibraryProduct("libmgl", :mgl),
    LibraryProduct("libmgl-qt", :mgl_qt)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Qt_jll", uuid="ede63266-ebff-546c-83e0-1c6fb6d0efc8"))
    Dependency(PackageSpec(name="libpng_jll", uuid="b53b4c65-9356-5827-b1ea-8c7a1a84506f"))
    Dependency(PackageSpec(name="Zlib_jll", uuid="83775a58-1f1d-513f-b197-d71354ab007a"))
    Dependency(PackageSpec(name="GSL_jll", uuid="1b77fbbe-d8ee-58f0-85f9-836ddc23a7a4"))
    Dependency(PackageSpec(name="Libglvnd_jll", uuid="7e76a0d4-f3c7-5321-8279-8d96eeed0f29"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"8.1.0")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
import Term: RenderableText, Spacer, vLine, hLine, cleantext, textlen, chars, Panel
import Term.layout: pad

@testset "Layout - pad" begin
    @test pad("aaa", 20, :left) == "aaa                 "
    @test pad("aaa", 20, :right) == "                 aaa"
    @test pad("aaa", 20, :center) == "        aaa         "
    @test pad("aaa", 10, 20) == "          aaa                    "
    p = Panel(; width=20, height=10)
    padded = pad(p.segments, 10, 10)
    @test padded[1].measure.w == 40
end

@testset "\e[34mlayout - spacer" begin
    sizes = [(22, 1), (44, 123), (21, 1), (4334, 232)]
    for (w, h) in sizes
        spacer = Spacer(w, h)
        @test spacer.measure.w == w
        @test spacer.measure.h == h
    end
end

@testset "\e[34mlayout - vLine " begin
    for h in [1, 22, 55, 11]
        line = vLine(h)
        @test length(line.segments) == h
        @test line.measure.h == h
    end

    lines = [(22, "bold"), (55, "red on_green")]
    for (h, style) in lines
        line = vLine(h; style = style)
        @test length(line.segments) == h
        @test line.measure.h == h
    end
    line = vLine(5; style="red")
    @test line.segments[1].text == "\e[31m│\e[39m"

    for box in (:MINIMAL_DOUBLE_HEAD, :DOUBLE, :ASCII, :DOUBLE_EDGE)
        @test vLine(22; box = box).measure.h == 22
    end

    panel = Panel(; width=20, height=5)
    @test length(vLine(panel).segments) == 5
    @test vLine().measure.h == displaysize(stdout)[1]
end

@testset "\e[34mlayout - hLine " begin
    for w in [1, 342, 433, 11, 22]
        line = hLine(w)
        @test length(line.segments) == 1
        @test textlen(line.segments[1].text) == w
        @test line.measure.w == w
    end

    for box in (:MINIMAL_DOUBLE_HEAD, :DOUBLE, :ASCII, :DOUBLE_EDGE)
        @test hLine(22; box = box).measure.w == 22
        @test hLine(22, "title"; box = box).measure.w == 22
    end

    for style in ("bold", "red on_green", "blue")
        @test textlen(hLine(11; style = style).segments[1].text) == 11
        @test textlen(hLine(11, "ttl"; style = style).segments[1].text) == 11
    end

    panel = Panel(; width=20, height=5)
    @test hLine().measure.w == displaysize(stdout)[2]
    @test textlen(hLine(panel).segments[1].text) == 20
end

@testset "\e[34mlayout - stack strings" begin
    s1 = "."^50
    s2 = ".\n"^5 * "."
    @test s1 / s2 isa String
end

function testlayout(p, w, h)
    _p = string(p)
    widths = textwidth.(cleantext.(split(_p, "\n")))
    @test length(unique(widths)) == 1

    @test p.measure.w == w
    @test textlen(cleantext(p.segments[1].text)) == w
    @test length(chars(cleantext(p.segments[1].text))) == w

    @test p.measure.h == h
    @test length(p.segments) == h
end

@testset "\e[34mlayout - renderable" begin
    r1 = RenderableText("."^100; width = 25)
    r2 = RenderableText("."^100; width = 50)

    r = r1 / r2
    @test r.measure.w == 50
    @test r.measure.h == 6

    h1 = hLine(22)
    h2 = hLine(33)
    @test (h1 / h2).measure.w == 33
    @test (h1 / h2).measure.h == 2

    r1 = RenderableText("."^100; width = 25)
    r2 = RenderableText("."^100; width = 50)

    r = r1 * r2
    @test r.measure.w == 75
    @test r.measure.h == 4

    # stack other renderables
    h1 = vLine(22)
    h2 = vLine(33)
    @test (h1 * h2).measure.w == 2
    @test (h1 * h2).measure.h == 33

end

@testset "\e[34mlayout - panels" begin
    p1 = Panel()
    p2 = Panel(; width=24, height=3)
    p3 = Panel("this [red]panel[/red]"^5, width=12)

    testlayout(p1 * p2, 112, 3)
    @test string(p1 * p2) == "\e[22m╭──────────────────────────────────────────────────────────────────────────────────────╮\e[22m\e[22m╭──────────────────────╮\e[22m\n\e[22m╰──────────────────────────────────────────────────────────────────────────────────────╯\e[22m\e[22m│\e[22m                      \e[22m│\e[22m\n                                                                                \e[22m╰──────────────────────╯\e[22m"

    testlayout(p1 / p2, 88, 5)
    @test string(p1 / p2) == "\e[22m╭──────────────────────────────────────────────────────────────────────────────────────╮\e[22m\n\e[22m╰──────────────────────────────────────────────────────────────────────────────────────╯\e[22m\n\e[22m╭──────────────────────╮\e[22m                                                                \n\e[22m│\e[22m                      \e[22m│\e[22m                                                                \n\e[22m╰──────────────────────╯\e[22m                                                                "

    testlayout(p2 * p1, 112, 3)
    testlayout(p2 / p1, 88, 5)

    testlayout(p1 * p2 * p3, 124, 12)
    testlayout(p1 / p2 / p3, 88, 17)
    testlayout(p3 * p1 * p2, 124, 12)
    testlayout(p3 / p1 / p2, 88, 17)
end


text:
#=
This script runs step-by-step through md-RQA analysis in Julia. It uses DynamicalSystems.jl for delay embedding and analysis.
=#

# ----- this script relies on the following packages -----
using Pkg
Pkg.add("DataFrames");
Pkg.add("CSV");
Pkg.add("DynamicalSystems"); # RQA stuffs
Pkg.add("StatsBase");
Pkg.add("Plots"); # plotting
Pkg.add("ElasticArrays"); # efficent matrix iteration
Pkg.add("LinearAlgebra") # diag caculation
Pkg.add("Distances") # matrix pairwise distances

using DataFrames, Plots, CSV, DynamicalSystems, StatsBase, ElasticArrays, LinearAlgebra, Distances;

# you should also include the follwing functions from this repo
include("mdFnn.jl")
include("mdEmbed.jl")

# Step 1 ----- load in an example dataset
data = DataFrame(CSV.File("exampleData.csv"));
#data = DataFrame(CSV.File("lorenz.csv", header=false));

# turn to matrix (and in this case only use first 3 columns)
ts_names = names(data[:,1:3]) |> permutedims; # use this for ami plotting
data = Matrix(data)[:,1:3];

# get the dimensions (number of columns) of the data
ncol = size(data,2);

# Step 2 ----- find the delay (tau)

# parameters:
nbins = 10;
maxlag = 200;

# get delay for each dimension
lags = zeros(ncol);
plotMI = zeros(maxlag+1, ncol);
for dimension in 1:ncol
    plotMI[:,dimension] = selfmutualinfo(data[:,dimension],0:1:maxlag; nbins=nbins);
    lags[dimension] = estimate_delay(data[:,dimension], "mi_min", 0:1:maxlag; nbins=nbins);
end

τ = maximum(lags) |> round |> Int;

# Step 3 ------ find the embedding dimension
# parameters
tau = τ;
maxEmb = 10;
numSamples = 500;
Rtol = 10;
Atol = 2;
abs_threshold = 1;

fnnPerc = mdFnn(data[:,:], tau, maxEmb, numSamples, Rtol, Atol);

# embedding dimension

# simple version
D = findall(fnnPerc .< abs_threshold) |> minimum;

#= a more comple version would be to consider if the rate of FNN function is slowing. This would be useful for cases where the FNN function never goes below the threshold but does "level off"

ex:

# the minimum difference between FNN of two consecutive D's
diff_threshold = 3

D_below_absolute_threshold = findall(fnnPerc .< abs_threshold)

if size(D_below_absolute_threshold)[1]>0
    D = (findall(fnnPerc .< abs_threshold) |> minimum);
else
    # if no values are below the absolute threshold,
    # then use the difference threshold

    D = (findall(abs.(diff(fnnPerc)) .< diff_threshold) |> minimum) + 1
end
=#

# Step 4 ------ embed timeseries
embed_TS = mdEmbed(data, D, τ);

# Step 5 ------ construct Recurrence Matrix
R = RecurrenceMatrix(embed_TS, .05; fixedrate = true);

# Step 6 ------ Recurrence quantification
rqaOUT = rqa(R; lmin = 50, thieller = 1);

# Step 7 ------- Construct output data frame

# get the obtained parameters
parameters = (delay = tau, embed = D);

# join parameters and rqaOUT
df_out = [merge(parameters, rqaOUT)] # be careful not to duplicate column names

#optionally save to csv
CSV.write("rqa_output.csv", df_out)

# Step 8 ------ Diagnostic Plots

## DO NOT ACCEPT ANY VALUES WITHOUT LOOKING AT THE PLOTS!!!!

# plot of the ami function
ami_plot = plot(plotMI[:,:], label = ts_names);
plot!(ami_plot, [τ], seriestype="vline", linestyle=:dash,label = "delay selected")

# false nearest neighbors plot
fnn_plot = plot(fnnPerc, label = "%FNN");
plot!(fnn_plot, [D], seriestype="vline", linestyle=:dash,label = "embed selected")

# recurrence plot
xs, ys = coordinates(R);
rec_plot = scatter(xs, ys, markersize = 0.1, markercolor = :black, legend = false);

# combine plots
l = @layout [a{0.5w} [b; c]]
diag_plots = plot(rec_plot, ami_plot,fnn_plot, layout = l);
plot!(diag_plots, size=(800,400))
png(diag_plots, "simple_diagnostic_plots_filename.png")


text:
insts = OrderedDict()
insts["minimal"] = [
    # tr neglog
    ((3, false, MatLogdetCone()),),
    ((3, true, MatLogdetCone()),),
    ((3, false, MatNegLog()),),
    ((3, true, MatNegLog()),),
    ((3, false, MatNegLogDirect()),),
    ((3, true, MatNegLogDirect()),),
    ((3, false, MatNegLogEigOrd()),),
    # negrtdet
    ((3, false, MatNegRtdet()),),
    ((3, true, MatNegRtdet()),),
    ((3, false, MatNegRtdetEFExp()),),
    ((3, true, MatNegRtdetEFExp()),),
    ((3, false, MatNegRtdetEFPow()),),
    ((3, true, MatNegRtdetEFPow()),),
    # tr negsqrtconj
    ((3, false, MatNegSqrtConj()),),
    ((3, false, MatNegSqrtConjDirect()),),
    ((3, false, MatNegSqrtConjEigOrd()),),
    ((3, true, MatNegSqrtConj()),),
    ((3, true, MatNegSqrtConjDirect()),),
    ((3, true, MatNegSqrtConjEigOrd()),),
    # tr negentropy
    ((3, false, MatNegEntropy()),),
    ((3, false, MatNegEntropyEigOrd()),),
    # tr negsqrt
    ((3, false, MatNegSqrt()),),
    ((3, false, MatNegSqrtEigOrd()),),
    # tr negpower01
    ((3, false, MatNegPower01(0.7)),),
    ((3, false, MatNegPower01(0.7)),),
    # tr power12
    ((3, false, MatPower12(1.3)),),
    ((3, false, MatPower12EigOrd(1.3)),),
    ]
insts["fast"] = [
    ((30, true, MatNegRtdet()),),
    ((15, false, MatNegRtdetEFExp()),),
    ((15, false, MatNegRtdetEFPow()),),
    ((40, false, MatNegSqrtConj()),),
    ((8, true, MatNegSqrtConjEigOrd()),),
    ((20, false, MatNegSqrtConjDirect()),),
    ((15, true, MatNegSqrtConjDirect()),),
    ((50, false, MatLogdetCone()),),
    ((30, false, MatNegLog()),),
    ((20, true, MatNegLog()),),
    ((6, false, MatNegLogEigOrd()),),
    ((15, false, MatNegLogDirect()),),
    ((10, true, MatNegLogDirect()),),
    ((30, false, MatNegEntropy()),),
    ((8, true, MatNegEntropyEigOrd()),),
    ((30, false, MatPower12(1.5)),),
    ((7, false, MatPower12EigOrd(1.5)),),
    ]
insts["various"] = [
    ((100, false, MatNegRtdet()),),
    ((60, true, MatNegRtdet()),),
    ((50, false, MatNegRtdetEFExp()),),
    ((30, true, MatNegRtdetEFPow()),),
    ((100, false, MatNegSqrtConj()),),
    ((12, false, MatNegSqrtConjEigOrd()),),
    ((50, false, MatNegSqrtConjDirect()),),
    ((30, true, MatNegSqrtConjDirect()),),
    ((150, false, MatLogdetCone()),),
    ((80, true, MatLogdetCone()),),
    ((100, false, MatNegLog()),),
    ((50, true, MatNegLog()),),
    ((18, false, MatNegLogEigOrd()), nothing, (default_tol_relax = 1000,)),
    ((80, false, MatNegLogDirect()),),
    ((40, true, MatNegLogDirect()),),
    ((75, false, MatNegEntropy()),),
    ((14, false, MatNegEntropyEigOrd()),),
    ((30, false, MatPower12(1.5)),),
    ((6, true, MatPower12EigOrd(1.5)),),
    ]
return (CovarianceEstJuMP, insts)


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL

function randomshapes()
	V,EV = Lar.randomcuboids(50, .4)
	V = GL.normalize2(V,flag=true)
	model2d = V,EV
	Sigma =  Lar.spaceindex(model2d);
	GL.VIEW([ GL.GLLines(V,EV, GL.COLORS[1]) ]);

	W,EW = Lar.fraglines(1.5,1.5,1.5)((V,EV))
	GL.VIEW([ GL.GLLines(W,EW, GL.COLORS[1]) ]);

	W = convert(Lar.Points, V')
	cop_EV = Lar.coboundary_0(EV::Lar.Cells)
	cop_EW = convert(Lar.ChainOp, cop_EV)
	V, copEV, copFE = Lar.Arrangement.planar_arrangement(W::Lar.Points, cop_EW::Lar.ChainOp)

	triangulated_faces = Lar.triangulate2D(V, [copEV, copFE])
	FVs = convert(Array{Lar.Cells}, triangulated_faces)
	W = convert(Lar.Points, V')
	GL.VIEW(GL.GLExplode(W,FVs,1.2,1.2,1.2,99,1));
	GL.VIEW(GL.GLExplode(W,FVs,1,1,1,99,1));
end

randomshapes()


text:
module JW3126PkgDev
export gen
using PkgTemplates
const PT = PkgTemplates

const gen = Template(;
       plugins=[
            PT.Git(; ssh=true, manifest=false),
            PT.Develop(),
            PT.TagBot(),
            # PT.TravisCI(windows=false, osx=false, linux=true),
            PT.Documenter{GitHubActions}(),
            PT.GitHubActions(),
            PT.Codecov(),
       ],
   )

end # module


text:
module problem108

description = """
Diophantine reciprocals I
Problem 108

In the following equation x, y, and n are positive integers.
1/x + 1/y = 1/n

For n = 4 there are exactly three distinct solutions:
1/5 + 1/20 = 1/4
1/6 + 1/12 = 1/4
1/8 + 1/8 = 1/4

What is the least value of n for which the number of distinct solutions exceeds one-thousand?

NOTE: This problem is an easier version of problem 110; it is strongly advised that you solve this one first.
"""

# Brute force check for 1/y + 1/x = 1/n for x in the range n+1 to 2n
# Integer y exists when (x-n) divides x*n, so no need to actually calulate the fractions.
function solutionCount(n :: Int)
  c = 0
  xn = n*(n+1)
  xmn = 1
  for x = n+1 : n * 2
    if xn % xmn == 0
      c += 1
    end
    xn += n
    xmn += 1
  end
  return c
end

using Base.Test

@test solutionCount(4) == 3

# all solutions take the form:
# a/(n*k) + (k-a)/(n*k) = 1/n,
# where a is a factor of n*k, and (k-a) is a factor of n*k
# I arrived at this solution by trial and error playing with the solutionCount above.
# Originally, I had only considered the case where a=1, reducing the problem to N having >1000 divisors.
# The solution to that problem requires a large value for Sigma0(n) which implies any prime factors.
@test solutionCount(2*2*3*3*5*7*11*13) > 1000

end

max_stars_repo_name: mbuhot/mbuhot-euler-solutions

max_issues_repo_name: mbuhot/mbuhot-euler-solutions

max_forks_repo_name: mbuhot/mbuhot-euler-solutions


text:
#------------------------------------------------------------------------------#
#                                    Module                                    #
#------------------------------------------------------------------------------#

"""
# Description

    module iso

Isotropic tensor module.

# Usage

This module is automatically `include`d by the `IsoOrthoTensor` package, and  is
not used in isolation.

```julia-repl
julia> using IsoOrthoTensor

julia> typeof(IOT.iso)
Module

julia> Δ(1) # Δ: U+394
2×2 Array{Int64,2}:
 1  0
 0  1

julia> Δ(2) # Δ: U+394
2×2×2×2 Array{Int64,4}:
[:, :, 1, 1] =
 3  0
 0  1

[:, :, 2, 1] =
 0  1
 1  0

[:, :, 1, 2] =
 0  1
 1  0

[:, :, 2, 2] =
 1  0
 0  3

```
"""
module iso

#------------------------------------------------------------------------------#
#                                   Imports                                    #
#------------------------------------------------------------------------------#

include("basicTensor.jl")
using .basicTensor

import Combinatorics

#------------------------------------------------------------------------------#
#                                   Exports                                    #
#------------------------------------------------------------------------------#

export i𝕔, 𝕔, Δ

#------------------------------------------------------------------------------#
#                             Auxiliary Functions                              #
#------------------------------------------------------------------------------#

"""
# Description

    i𝕔(
        OPD::Tuple{Vararg{Int64,N}} where N,
        FID::Tuple{Vararg{Int64,N}} where N;
        PRE = ([()], 1)
    )::Array{Tuple{Vararg{Int64,N}} where N,1} # 𝕔: U+1d554

Returns a one-dimensional array with tensor index combinations (as NTuples) that
encode summation terms for a combinatorial type of nonstandard  tensor  product,
referred to by the '⊛': Unicode U+229b symbol in some Lattice  Boltzmann  theory
literature, as [1].

`OPD` is an NTuple containing the operand dimensions, like (2, 2),  and  has  as
many elements as there are operands in the  combinatorial  type  of  nonstandard
tensor product.

`FID` is an NTuple containing the fixed indices, i.e., indices that take a fixed
position in all terms, like (1,), indicating that the first index is kept  fixed
in all summation terms.

`PRE` is **not** a user input and should be left alone. It  serves  to  pass  on
information across recursive calls.

# Usage

Nonstandard combinatorial tensor product between two rank-2 tensors, keeping the
first index fixed: arguments `(2, 2)` and  `(1,)`  to  the  `i𝕔`  function  —  a
summation of `3 choose 1 = 3` products between the operands:

```julia-repl
julia> using IsoOrthoTensor

julia> i𝕔((2, 2), (1,)) # 𝕔: U+1d554
3-element Array{Tuple{Vararg{Int64,N}} where N,1}:
 (1, 2, 3, 4)
 (1, 3, 2, 4)
 (1, 4, 2, 3)

```

Nonstandard combinatorial tensor product between three rank-2  tensors,  keeping
the first and sixth indices fixed: arguments `(2, 2, 2)` and  `(1,  6)`  to  the
`i𝕔` function — a summation of `(4 choose 1)*(3 choose 2) = 4*3 =  12`  products
between the operands:

```julia-repl
julia> i𝕔((2, 2, 2), (1, 6)) # 𝕔: U+1d554
12-element Array{Tuple{Vararg{Int64,N}} where N,1}:
 (1, 2, 3, 4, 5, 6)
 (1, 2, 3, 5, 4, 6)
 (1, 2, 4, 5, 3, 6)
 (1, 3, 2, 4, 5, 6)
 (1, 3, 2, 5, 4, 6)
 (1, 3, 4, 5, 2, 6)
 (1, 4, 2, 3, 5, 6)
 (1, 4, 2, 5, 3, 6)
 (1, 4, 3, 5, 2, 6)
 (1, 5, 2, 3, 4, 6)
 (1, 5, 2, 4, 3, 6)
 (1, 5, 3, 4, 2, 6)

```

# References

[1]: K. K.  Mattila,  L.  A.  Hegele  Júnior,  P.  C.  Philippi,  “High-Accuracy
Approximation of High-Rank Derivatives: Isotropic Finite  Differences  Based  of
Lattice-Boltzmann Stencils,” The Scientific World Journal, vol. 2014, article ID
142907, 16 pages, 2014.
"""
function i𝕔(
    OPD::Tuple{Vararg{Int64,N}} where N,
    FID::Tuple{Vararg{Int64,N}} where N;
    PRE = ([()], 1)
)::Array{Tuple{Vararg{Int64,N}} where N,1} # 𝕔: U+1d554
    RED = sum(OPD) + length(PRE[1][1]) # REsult Dimensions
    VID = Tuple(i for i in 1:RED if !(i in FID)) # Variable InDices
    ARR = Array{NTuple{N,Int} where N,1}()
    BEG = PRE[2]
    RNG = BEG:BEG+OPD[1]-1
    SLO = Tuple(i for i in RNG if !(i in FID)) # free SLOts in current operand
    for i in PRE[1]
        ALW = Tuple(k for k in VID if !(k in i)) # ALoWed free indices
        for j in Combinatorics.multiset_combinations(ALW, length(SLO))
            ITM = Int[]
            l = 1
            for k in RNG
                if k in FID
                    push!(ITM, k)
                else
                    push!(ITM, j[l])
                    l += 1
                end
            end
            push!(ARR, Tuple(vcat(i..., ITM...)))
        end
    end
    if length(OPD) <= 1
        return ARR
    else
        return i𝕔(OPD[2:end], FID, PRE = (ARR, BEG+OPD[1]))
    end
end

"""
# Description

    function 𝕔(
        OPS::Tuple{Vararg{Tensor{Int64},N}} where N
        FID::Tuple{Vararg{Int64,N}} where N;
        D::Int64 = 2
    )::Tensor{Int64} # 𝕔: U+1d554

Performs  a  combinatorial  type  of  nonstandard  tensor  product  between  the
operands, which are elements of the `OPS` `Tuple`,  keeping  the  indices  `FID`
fixed in a `D`-dimensional Euclidean space, and returns  the  resulting  tensor.
The combinatorial nonstandard tensor product performed is  referred  to  by  the
'⊛': Unicode U+229b symbol in some Lattice Boltzmann theory literature, as [1].

`OPS` is an NTuple containing the operand tensors, like `(δ, δ)`  —  a  pair  of
Kronecker-delta tensors stored in the temporary `δ` identifier.

`FID` is an NTuple containing the fixed indices, i.e., indices that take a fixed
position in all terms, like (1,), indicating that the first index is kept  fixed
in all summation terms.

`D` is the dimensionality of the Euclidean space in which  the  tensor  operands
apply.

# Usage

Calculate the `𝝙⁽ⁿ⁾` tensor of order `2n` that is isotropic with respect to  all
of its `2n` indices [1] with `n=2` through the nonstandard product

    Δ⁽²⁾αβγε = δαβ*δγε + δαγ*δβε + δαε*δβγ,

in which terms combine the three free indices  `βγε`  while  keeping  the  first
index `α` fixed.

```julia-repl
julia> 𝕔((K(2), K(2)), ()) # 𝕔: U+1d554
2×2×2×2 Array{Int64,4}:
[:, :, 1, 1] =
 6  0
 0  2

[:, :, 2, 1] =
 0  2
 2  0

[:, :, 1, 2] =
 0  2
 2  0

[:, :, 2, 2] =
 2  0
 0  6

```

# References

[1]: K. K.  Mattila,  L.  A.  Hegele  Júnior,  P.  C.  Philippi,  “High-Accuracy
Approximation of High-Rank Derivatives: Isotropic Finite  Differences  Based  of
Lattice-Boltzmann Stencils,” The Scientific World Journal, vol. 2014, article ID
142907, 16 pages, 2014.
"""
function 𝕔(
    OPS::Tuple{Vararg{Array{Int64},N}} where N,
    FID::Tuple{Vararg{Int64,N}} where N;
    D::Int64 = 2
)::Tensor{Int64} # 𝕔: U+1d554
    OPD = Tuple(ndims(A) for A in OPS) # OPerand Dimensions
    RED = sum(OPD) # REsult Dimensions
    RET = fill(zero(OPS[1][1]), Tuple(D for i in 1:RED)) # RETurn tensor
    ONE = one(OPS[1][1])
    for idx in Base.product(fill(1:D, RED)...) # Return tensor index
        for IDP in i𝕔(OPD, FID) # InDex Permutations
            PID = [idx[perm] for perm in IDP] # Permuted InDices
            prod, BEG = ONE, 1
            for fact in OPS
                SIZ = ndims(fact)
                prod *= fact[PID[BEG:BEG-1+SIZ]...]
                BEG += SIZ
            end
            RET[idx...] += prod
        end
    end
    return RET
end

#------------------------------------------------------------------------------#
#                                User Functions                                #
#------------------------------------------------------------------------------#

"""
# Description

    Δ(n::Int64; D::Int64 = 2)::Tensor{Int64} # Δ: U+394

Computes and returns an `n`-th order Isotropic Tensor (of `Int64` components) in
a `D`-dimensional Euclidean space, checking bounds on `D`, i.e.,  whether  `D  ∈
[1, 3]`, and on `n`, i.e., whether `n ∈ [0, ∞)`.

```julia-repl
julia> Δ(1) # Δ: U+394
2×2 Array{Int64,2}:
 1  0
 0  1

julia> Δ(2) # Δ: U+394
2×2×2×2 Array{Int64,4}:
[:, :, 1, 1] =
 3  0
 0  1

[:, :, 2, 1] =
 0  1
 1  0

[:, :, 1, 2] =
 0  1
 1  0

[:, :, 2, 2] =
 1  0
 0  3

```
"""
function Δ(n::Int64; D::Int64 = 2)::Tensor{Int64} # Δ: U+394
    # Validation
    if D <= 0 || D >= 4
        throw(DomainError("dimension D = $D outside the valid domain [1, 3]"))
    elseif n < 0
        throw(DomainError("order n = $n outside the valid domain [0, ∞)"))
        # ∞: U+221e
    end
    # Execution
    if D == 1 || n == 0
        return 1
    else
        return n == 1 ? K(D) : 𝕔((K(D), Δ(n-1, D = D)), (1,), D = D)
    end
end

end


text:
struct SubProblem{T <: AbstractFloat, A <: AbstractVector, S <: LQSolver, PT <: PenaltyTerm}
    id::Int
    π::T
    solver::S
    penalty::PT
    c::A
    x::A
    y::A
    ρ::A
    optimvector::A

    function (::Type{SubProblem})(model::JuMP.Model,
                                  id::Integer,
                                  π::AbstractFloat,
                                  xdim::Integer,
                                  subsolver::MPB.AbstractMathProgSolver,
                                  penalty::PenaltyTerm)
        solver = LQSolver(model, subsolver)
        solver()
        optimvector = getsolution(solver)
        x₀ = optimvector[1:xdim]
        y₀ = optimvector[xdim+1:end]
        T = promote_type(eltype(optimvector),Float32)
        c_ = convert(AbstractVector{T},JuMP.prepAffObjective(model))
        c_ *= model.objSense == :Min ? 1 : -1
        x₀_ = convert(AbstractVector{T},x₀)
        y₀_ = convert(AbstractVector{T},y₀)
        optimvector_ = convert(AbstractVector{T},optimvector)
        A = typeof(x₀_)
        S = typeof(solver)
        PT = typeof(penalty)
        subproblem = new{T,A,S,PT}(id,
                                   π,
                                   solver,
                                   penalty,
                                   c_,
                                   x₀_,
                                   y₀_,
                                   zero(x₀_),
                                   optimvector_)
        initialize_penaltyterm!(penalty, solver, x₀_)
        return subproblem
    end

    function (::Type{SubProblem})(model::JuMP.Model,
                                  id::Integer,
                                  π::AbstractFloat,
                                  x₀::AbstractVector,
                                  y₀::AbstractVector,
                                  subsolver::MPB.AbstractMathProgSolver,
                                  penalty::PenaltyTerm)
        T = promote_type(eltype(x₀),eltype(y₀),Float32)
        c_ = convert(AbstractVector{T},JuMP.prepAffObjective(model))
        c_ *= model.objSense == :Min ? 1 : -1
        x₀_ = convert(AbstractVector{T},x₀)
        y₀_ = convert(AbstractVector{T},y₀)
        A = typeof(x₀_)
        solver = LQSolver(model, subsolver)
        S = typeof(solver)
        PT = typeof(penalty)
        subproblem = new{T,A,S,PT}(id,
                                   π,
                                   solver,
                                   penalty,
                                   c_,
                                   x₀_,
                                   y₀_,
                                   zero(x₀_),
                                   [x₀_...,y₀_...])
        initialize_penaltyterm!(penalty, solver, x₀_)
        return subproblem
    end
end

function update_subproblem!(subproblem::SubProblem, ξ::AbstractVector, r::AbstractFloat)
    subproblem.ρ .= subproblem.ρ + r*(subproblem.x - ξ)
end
update_subproblems!(subproblems::Vector{<:SubProblem}, ξ::AbstractVector, r::AbstractFloat) = map(prob -> update_subproblem!(prob,ξ,r), subproblems)

function reformulate_subproblem!(subproblem::SubProblem, ξ::AbstractVector, r::AbstractFloat)
    model = subproblem.solver.lqmodel
    # Cache initial cost
    c = copy(subproblem.c)
    # Linear part
    c[1:length(ξ)] += subproblem.ρ
    # Update penalty
    update_penaltyterm!(subproblem.penalty, subproblem.solver, c, r, ξ)
    return nothing
end
reformulate_subproblems!(subproblems::Vector{<:SubProblem}, ξ::AbstractVector, r::AbstractFloat) = map(prob -> reformulate_subproblem!(prob,ξ,r), subproblems)

function get_objective_value(subproblem::SubProblem)
    return subproblem.π*subproblem.c⋅subproblem.optimvector
end

function get_solution(subproblem::SubProblem)
    return copy(subproblem.y), getredcosts(subproblem.solver), getduals(subproblem.solver)
end

function fill_submodel!(submodel::JuMP.Model, subproblem::SubProblem, nrows::Integer, ncols::Integer)
    fill_submodel!(submodel, get_solution(subproblem)..., nrows, ncols)
end

function fill_submodel!(submodel::JuMP.Model, x::AbstractVector, μ::AbstractVector, λ::AbstractVector, nrows::Integer, ncols::Integer)
    submodel.colVal = x[1:length(submodel.colVal)]
    submodel.redCosts = μ[ncols+1:ncols+length(submodel.colVal)]
    submodel.linconstrDuals = λ[nrows+1:nrows+length(submodel.linconstr)]
    submodel.objVal = JuMP.prepAffObjective(submodel)⋅x[1:length(submodel.colVal)]
end

function (subproblem::SubProblem)(ξ::AbstractVector)
    solve_penalized!(subproblem.penalty, subproblem.solver, subproblem.optimvector, subproblem.x, ξ)
    solvestatus = status(subproblem.solver)
    if solvestatus == :Optimal
        ncols = length(subproblem.optimvector)
        xdim = length(subproblem.x)
        subproblem.optimvector[1:ncols] = getsolution(subproblem.solver)[1:ncols]
        subproblem.x .= subproblem.optimvector[1:xdim]
        subproblem.y .= subproblem.optimvector[xdim+1:ncols]
        return get_objective_value(subproblem)
    elseif solvestatus == :Infeasible
        return Inf
    elseif solvestatus == :Unbounded
        return -Inf
    else
        error(@sprintf("Subproblem %d was not solved properly, returned status code: %s", subproblem.id, string(solvestatus)))
    end
end