text:
## Convert_a functions for the mdp
# Converts from vector to an action
function POMDPs.convert_a(::Type{Vector{Disturbance}}, avec::AbstractArray{Float64}, mdp::AdversarialDrivingMDP)
    a = Disturbance[]
    index = 1
    for agent in adversaries(mdp)
        d = agent.vec_to_disturbance(avec[index : index + agent.disturbance_dim - 1])
        index += agent.disturbance_dim
        push!(a, d)
    end
    a
end

# Convert from action to a vector
function POMDPs.convert_a(::Type{AbstractArray}, a::Vector{Disturbance}, mdp::AdversarialDrivingMDP)
    vec = Vector{Float64}(undef, sum([agent.disturbance_dim for agent in adversaries(mdp)]))
    index = 1
    for i = 1:length(adversaries(mdp))
        agent = adversaries(mdp)[i]
        vec[index : index + agent.disturbance_dim - 1] .= agent.disturbance_to_vec(a[i])
        index += agent.disturbance_dim
    end
    vec
end

## PedestrianControl Disturbances
const PEDESTRIAN_DISTURBANCE_DIM = 5

# Converts PedestrianControl disturbance to a vector
function PedestrianControl_to_vec(pc::PedestrianControl)
    Float64[pc.da..., pc.noise.pos..., pc.noise.vel]
end

# Converts from a vector to a PedestrianConrol Disturbance
function vec_to_PedestrianControl(arr::AbstractArray)
    @assert length(arr) == PEDESTRIAN_DISTURBANCE_DIM
    da = VecE2(arr[1], arr[2])
    noise_pos = VecE2(arr[3], arr[4])
    noise_v = arr[5]
    PedestrianControl(VecE2(0., 0.), da, Noise(noise_pos, noise_v))
end

## BlinkerVehicle Disturbances
const BLINKERVEHICLE_DISTURBANCE_DIM = 5

# Converts PedestrianControl disturbance to a vector
function BlinkerVehicleControl_to_vec(bv::BlinkerVehicleControl)
    g = Float64(bv.toggle_goal)
    b = Float64(bv.toggle_blinker)
    Float64[bv.da, g, b, bv.noise.pos[1], bv.noise.vel]
end

# Converts from a vector to a PedestrianConrol Disturbance
function vec_to_BlinkerVehicleControl(arr::AbstractArray)
    @assert length(arr) == BLINKERVEHICLE_DISTURBANCE_DIM
    da = arr[1]
    toggle_goal = arr[2] > 0
    toggle_blinker = arr[3] > 0
    noise_pos = VecE2(arr[4], 0.)
    noise_v = arr[5]
    BlinkerVehicleControl(0., da, toggle_goal, toggle_blinker, Noise(noise_pos, noise_v))
end

## Continuous action spaces
get_actions(m::Vector{Sampleable}) = error("Error! `actions` called on continuous action space")
get_actionindex(m::Vector{Sampleable}, a::Vector{Disturbance}) = error("Error! `actionindex` called on continuous action space")
function Distributions.logpdf(m::Vector{Sampleable}, a::Vector{Disturbance}, mdp::AdversarialDrivingMDP)
    avec = convert_a(AbstractArray, a, mdp)
    sum([logpdf(m[i], avec[i]) for i=1:length(avec)])
end
Base.rand(rng::AbstractRNG, m::Vector{Sampleable}, mdp::AdversarialDrivingMDP) = convert_a(Vector{Disturbance}, [rand(rng, d) for d in m], mdp)
Base.rand(rng::AbstractRNG, m::Vector{Sampleable}, convert_fn::Function) = convert_fn([rand(rng, d) for d in m])

## Discrete BlinkerVehicle Actions
struct DiscreteActionModel
    actions::Vector{Vector{Disturbance}}
    actionindex::Dict{Vector{Disturbance}, Int}
    probs::Vector{Float64}
end

function DiscreteActionModel(actions::Vector{Vector{Disturbance}}, probs::Vector{Float64})
    aind = Dict{Vector{Disturbance}, Int}(a => i for (a, i) in zip(actions, 1:length(actions)))
    DiscreteActionModel(actions, aind, probs)
end

get_actions(m::DiscreteActionModel) = m.actions
get_actionindex(m::DiscreteActionModel, a::Vector{Disturbance}) = m.actionindex[a]
Distributions.logpdf(m::DiscreteActionModel, a::Vector{Disturbance}, mdp::AdversarialDrivingMDP) = log(m.probs[findfirst(m.actions .== [a])])
Base.rand(rng::AbstractRNG, m::DiscreteActionModel, mdp::AdversarialDrivingMDP) = m.actions[rand(rng, Categorical(m.probs))]

## Combine disturbance_models
# Automatically determine the type of disturbance (continuous or discrete)
combine_disturbance_models(adversaries::Vector{Agent}) = adversaries[1].disturbance_model isa DiscreteActionModel ? combine_discrete(adversaries) : combine_continuous(adversaries)

# Combines continuous disturbance models
function combine_continuous(adversaries::Vector{Agent})
    m = Vector{Sampleable}()
    for a in adversaries
        push!(m, a.disturbance_model...)
    end
    m
end

# Combines discrete disturbance models
function combine_discrete(adversaries::Vector{Agent})
    actions = Vector{Disturbance}[]
    action_probabilities = Float64[]

    # Add the baseline action where all agents do nothing
    base_action = Disturbance[a.disturbance_model.actions[1][1] for a in adversaries]
    aprob = sum([a.disturbance_model.probs[1] for a in adversaries])
    push!(actions, base_action)
    push!(action_probabilities, aprob)

    # Loops through the rest of the actions
    for adv_i=1:length(adversaries)
        adv = adversaries[adv_i]
        for (act, prob) in zip(adv.disturbance_model.actions[2:end], adv.disturbance_model.probs[2:end])
            new_action = deepcopy(base_action)
            new_action[adv_i] = act[1]
            push!(actions, new_action)
            push!(action_probabilities, prob)
        end
    end
    action_probabilities ./= sum(action_probabilities)
    DiscreteActionModel(actions, action_probabilities)
end

## default actions
# Function to setup default actions for the blinker verhicle
get_bv_actions(med_accel = 1.5, large_accel = 3.0, med_prob = 1e-2, large_prob = 1e-3) = DiscreteActionModel(
                      Vector{Vector{Disturbance}}([ [BlinkerVehicleControl(0, 0., false, false, Noise())],
                        [BlinkerVehicleControl(0, -large_accel, false, false, Noise())],
                        [BlinkerVehicleControl(0, -med_accel, false, false, Noise())],
                        [BlinkerVehicleControl(0, med_accel, false, false, Noise())],
                        [BlinkerVehicleControl(0, large_accel, false, false, Noise())],
                        [BlinkerVehicleControl(0, 0., true, false, Noise())], # toggle goal
                        [BlinkerVehicleControl(0, 0., false, true, Noise())] # toggle blinker
                       ]),
                        [1 - (4*large_prob + 2*med_prob), large_prob, med_prob, med_prob, large_prob, large_prob, large_prob]
                        )

function get_rand_bv_actions(rng; med_accel = 1.5, large_accel = 3.0)
    probs = 10 .^ [rand(rng, Uniform(-4, -2)),
                   rand(rng, Uniform(-3, -1)),
                   rand(rng, Uniform(-3, -1)),
                   rand(rng, Uniform(-4, -2)),
                   rand(rng, Uniform(-4, -2)),
                   rand(rng, Uniform(-4, -2))]
    leftover = 1.0 - sum(probs)
    DiscreteActionModel( Vector{Vector{Disturbance}}([ [BlinkerVehicleControl(0, 0., false, false, Noise())],
                            [BlinkerVehicleControl(0, -large_accel, false, false, Noise())],
                            [BlinkerVehicleControl(0, -med_accel, false, false, Noise())],
                            [BlinkerVehicleControl(0, med_accel, false, false, Noise())],
                            [BlinkerVehicleControl(0, large_accel, false, false, Noise())],
                            [BlinkerVehicleControl(0, 0., true, false, Noise())], # toggle goal
                            [BlinkerVehicleControl(0, 0., false, true, Noise())] # toggle blinker
                           ]),
                            [leftover, probs...]
                            )
end

# TODO Add default distributions as needed


text:
@testset "fsm_retx_win.jl" begin
    DEFAULT_CAPACITY = 65000
    TIMEOUT_DFLT = 1000
    NREPS = 64
    MAX_PAYLOAD_SIZE = 1460
    MAX_RETX_ATTEMPTS = 8
    @testset "multiple segments with intervening ack" begin
        tx_ackno = WrappingInt32(rand(UInt32))
        conn = fsm_in_established(tx_ackno - 1, tx_ackno - 1; cap=DEFAULT_CAPACITY)
        d1 = "asdf"
        d2 = "qwer"

        write!(conn, d1)
        tick!(conn, 1)
        tick!(conn, 20)
        write!(conn, d2)
        tick!(conn, 1)
        @test conn.segments_out |> !isempty
        expect_seg(conn; data=d1, ack=true, payload_size=length(d1))
        expect_one_seg(conn; data=d2, ack=true, payload_size=length(d1))
        tick!(conn, TIMEOUT_DFLT - 23)
        expect_no_seg(conn)

        tick!(conn, 4)
        expect_one_seg(conn; data=d1, ack=true, payload_size=length(d1))
        tick!(conn, 2TIMEOUT_DFLT - 2)
        expect_no_seg(conn)
        send_ack!(conn, tx_ackno, tx_ackno + 4)
        tick!(conn, TIMEOUT_DFLT - 2)
        expect_no_seg(conn)
        tick!(conn, 3)
        expect_one_seg(conn; data=d2, ack=true, payload_size=length(d1))
    end

    @testset "multiple segments without intervening ack" begin
        tx_ackno = WrappingInt32(rand(UInt32))
        conn = fsm_in_established(tx_ackno - 1, tx_ackno - 1; cap=DEFAULT_CAPACITY)
        d1 = "asdf"
        d2 = "qwer"

        write!(conn, d1)
        tick!(conn, 1)
        tick!(conn, 20)
        write!(conn, d2)
        tick!(conn, 1)

        @test conn.segments_out |> !isempty
        expect_seg(conn; data=d1, ack=true, payload_size=length(d1))
        expect_one_seg(conn; data=d2, ack=true, payload_size=length(d1))
        tick!(conn, TIMEOUT_DFLT - 23)
        expect_no_seg(conn)

        tick!(conn, 4)
        expect_one_seg(conn; data=d1, ack=true, payload_size=length(d1))
        tick!(conn, 2TIMEOUT_DFLT - 2)
        expect_no_seg(conn)
        tick!(conn, 3)
        expect_one_seg(conn; data=d1, ack=true, payload_size=length(d1))
    end

    @testset "check that ACK of new data resets exponential backoff and restarts timer" begin
        backoff_test = function (num_backoffs::Int)
            tx_ackno = WrappingInt32(rand(UInt32))
            conn = fsm_in_established(tx_ackno - 1, tx_ackno - 1; cap=DEFAULT_CAPACITY)
            d1 = "asdf"
            d2 = "qwer"

            write!(conn, d1)
            tick!(conn, 1)
            tick!(conn, 20)
            write!(conn, d2)
            tick!(conn, 1)

            @test conn.segments_out |> !isempty
            expect_seg(conn; data=d1, ack=true, payload_size=length(d1))
            expect_one_seg(conn; data=d2, ack=true, payload_size=length(d1))
            tick!(conn, TIMEOUT_DFLT - 23)
            expect_no_seg(conn)

            tick!(conn, 4)
            expect_one_seg(conn; data=d1, ack=true, payload_size=length(d1))
            for i in 1:num_backoffs-1
                tick!(conn, (TIMEOUT_DFLT << i) - i)
                expect_no_seg(conn)
                tick!(conn, i)
                expect_one_seg(conn; ack=true, data=d1, payload_size=length(d1))
            end
            send_ack!(conn, tx_ackno, tx_ackno + 4)
            tick!(conn, TIMEOUT_DFLT - 2)
            expect_no_seg(conn)
            tick!(conn, 3)
            expect_one_seg(conn; ack=true, data=d2, payload_size=length(d2))
        end

        for i in 1:MAX_RETX_ATTEMPTS
            backoff_test(i - 1)
        end
    end
end

ext: jl

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
using Xpress
using MathOptInterface
using Test

const MOI = MathOptInterface
const MOIT = MathOptInterface.Test

const MOIT = MOI.Test
const MOIU = MOI.Utilities

const OPTIMIZER = Xpress.Optimizer(OUTPUTLOG = 0)
const OPTIMIZER_2 = Xpress.Optimizer(OUTPUTLOG = 0)

# Xpress can only obtain primal and dual rays without presolve. Check more on
# https://www.fico.com/fico-xpress-optimization/docs/latest/solver/optimizer/HTML/XPRSgetprimalray.html

const BRIDGED_OPTIMIZER = MOI.Bridges.full_bridge_optimizer(
    OPTIMIZER, Float64)
const BRIDGED_CERTIFICATE_OPTIMIZER =
    MOI.Bridges.full_bridge_optimizer(
        Xpress.Optimizer(OUTPUTLOG = 0, PRESOLVE = 0), Float64)

const CONFIG = MOIT.TestConfig()
const CONFIG_LOW_TOL = MOIT.TestConfig(atol = 1e-3, rtol = 1e-3)

@testset "MOI interface" begin
    optimizer = Xpress.Optimizer(OUTPUTLOG = 0)
    @test MOI.get(optimizer, MOI.RawParameter("logfile")) == ""
    optimizer = Xpress.Optimizer(OUTPUTLOG = 0, logfile = "output.log")
    @test MOI.get(optimizer, MOI.RawParameter("logfile")) == "output.log"
    @test MOI.set(optimizer, MOI.TimeLimitSec(), 100) === nothing
    @test MOI.set(optimizer, MOI.TimeLimitSec(), 3600.0) === nothing

    @test MOI.get(optimizer, MOI.RawParameter("MIPRELSTOP")) == 0.0001
    MOI.set(optimizer, MOI.RawParameter("MIPRELSTOP"), 0.123)
    @test MOI.get(optimizer, MOI.RawParameter("MIPRELSTOP")) == 0.123

    @test MOI.get(optimizer, MOI.RawParameter("MPSOBJNAME")) == ""
    MOI.set(optimizer, MOI.RawParameter("MPSOBJNAME"), "qwerty")
    @test MOI.get(optimizer, MOI.RawParameter("MPSOBJNAME")) == "qwerty"

    @test MOI.get(optimizer, MOI.RawParameter("PRESOLVE")) == 1
    MOI.set(optimizer, MOI.RawParameter("PRESOLVE"), 3)
    @test MOI.get(optimizer, MOI.RawParameter("PRESOLVE")) == 3

    @show MOI.get(optimizer, MOI.RawParameter("XPRESSVERSION"))
    @test MOI.get(optimizer, MOI.RawParameter("MATRIXNAME")) == ""
    @test MOI.get(optimizer, MOI.RawParameter("SUMPRIMALINF")) == 0.0
    @test MOI.get(optimizer, MOI.RawParameter("NAMELENGTH")) == 8
end

@testset "SolverName" begin
    @test MOI.get(OPTIMIZER, MOI.SolverName()) == "Xpress"
end

@testset "supports_default_copy_to" begin
    @test MOIU.supports_default_copy_to(OPTIMIZER, true)
end

@testset "Unit Tests Constraints" begin
    MOIT.basic_constraint_tests(OPTIMIZER, CONFIG)
end

@testset "Unit Tests" begin
    MOIT.unittest(BRIDGED_OPTIMIZER, CONFIG, #= excludes =# String[
        # TODO: fix errors
        # does not work with caching optimizer
        "delete_soc_variables",

        # These tests fail due to tolerance issues; tested below.
        "solve_qcp_edge_cases",
        "solve_qp_edge_cases",
        "solve_qp_zero_offdiag",

        # These tests require extra parameters to obtain certificates.
        "solve_farkas_equalto_upper",
        "solve_farkas_equalto_lower",
        "solve_farkas_lessthan",
        "solve_farkas_greaterthan",
        "solve_farkas_interval_upper",
        "solve_farkas_interval_lower",
        "solve_farkas_variable_lessthan",
        "solve_farkas_variable_lessthan_max",
       ],
    )
    # certificates
    MOIT.solve_farkas_equalto_upper(BRIDGED_CERTIFICATE_OPTIMIZER,CONFIG)
    MOIT.solve_farkas_equalto_lower(BRIDGED_CERTIFICATE_OPTIMIZER,CONFIG)
    MOIT.solve_farkas_lessthan(BRIDGED_CERTIFICATE_OPTIMIZER,CONFIG)
    MOIT.solve_farkas_greaterthan(BRIDGED_CERTIFICATE_OPTIMIZER,CONFIG)
    MOIT.solve_farkas_interval_upper(BRIDGED_CERTIFICATE_OPTIMIZER,CONFIG)
    MOIT.solve_farkas_interval_lower(BRIDGED_CERTIFICATE_OPTIMIZER,CONFIG)
    MOIT.solve_farkas_variable_lessthan(BRIDGED_CERTIFICATE_OPTIMIZER,CONFIG)
    MOIT.solve_farkas_variable_lessthan_max(BRIDGED_CERTIFICATE_OPTIMIZER,CONFIG)
    # tolerance
    MOIT.solve_qcp_edge_cases(BRIDGED_OPTIMIZER, CONFIG_LOW_TOL)
    MOIT.solve_qp_edge_cases(BRIDGED_OPTIMIZER, CONFIG_LOW_TOL)
    MOIT.solve_qp_zero_offdiag(BRIDGED_OPTIMIZER, CONFIG_LOW_TOL)
    # MOIT.delete_soc_variables(OPTIMIZER, CONFIG_LOW_TOL)
    MOIT.modificationtest(BRIDGED_OPTIMIZER, CONFIG)
end

@testset "Linear tests" begin
    MOIT.contlineartest(
        BRIDGED_OPTIMIZER,
        MOIT.TestConfig(
            basis = true,
            infeas_certificates = false
        ), [
            # These tests require extra parameters to obtain certificates.
            "linear8a", "linear8b", "linear8c","linear12",
        ]
    )
    MOIT.linear8atest(BRIDGED_CERTIFICATE_OPTIMIZER, CONFIG)
    MOIT.linear8btest(BRIDGED_CERTIFICATE_OPTIMIZER, CONFIG)
    MOIT.linear8ctest(BRIDGED_CERTIFICATE_OPTIMIZER, CONFIG)
    MOIT.linear12test(BRIDGED_CERTIFICATE_OPTIMIZER, CONFIG)
end

@testset "Quadratic tests" begin
    MOIT.contquadratictest(
        BRIDGED_OPTIMIZER,
        MOIT.TestConfig(atol = 1e-3, rtol = 1e-3),
        [
            # Xpress does NOT support ncqcp.
            "ncqcp",
        ]
    )
end

@testset "Conic tests" begin
    MOIT.lintest(BRIDGED_OPTIMIZER, CONFIG, [
        # These tests require extra parameters to obtain certificates.
        "lin3", "lin4"
    ])
    MOIT.lin3test(BRIDGED_CERTIFICATE_OPTIMIZER, CONFIG)
    MOIT.lin4test(BRIDGED_CERTIFICATE_OPTIMIZER, CONFIG)
    MOIT.soctest(BRIDGED_OPTIMIZER, MOIT.TestConfig(duals = false, atol=1e-3), ["soc3"])
    MOIT.soc3test(
        BRIDGED_OPTIMIZER,
        MOIT.TestConfig(duals = false, infeas_certificates = false, atol = 1e-3)
    )
    MOIT.rsoctest(BRIDGED_OPTIMIZER, MOIT.TestConfig(duals = false, atol=1e-3), ["rotatedsoc3"])
    MOIT.rotatedsoc3test(BRIDGED_OPTIMIZER, MOIT.TestConfig(duals = false, atol=1e-2))
    MOIT.geomeantest(BRIDGED_OPTIMIZER, MOIT.TestConfig(duals = false, atol=1e-3))
end

@testset "Integer Linear tests" begin
    MOIT.intlineartest(BRIDGED_OPTIMIZER, CONFIG)
end

@testset "ModelLike tests" begin
    @testset "default_objective_test" begin
        MOIT.default_objective_test(OPTIMIZER)
    end

    @testset "default_status_test" begin
        MOIT.default_status_test(OPTIMIZER)
    end

    @testset "nametest" begin
        MOIT.nametest(BRIDGED_OPTIMIZER)
    end

    @testset "validtest" begin
        MOIT.validtest(OPTIMIZER)
    end

    @testset "emptytest" begin
        MOIT.emptytest(BRIDGED_OPTIMIZER)
    end

    @testset "orderedindicestest" begin
        MOIT.orderedindicestest(OPTIMIZER)
    end

    @testset "copytest" begin
        BRIDGED_OPTIMIZER_2 = MOI.Bridges.full_bridge_optimizer(
            OPTIMIZER_2, Float64
        )
        MOIT.copytest(BRIDGED_OPTIMIZER, BRIDGED_OPTIMIZER_2)
    end
end

@testset "IIS tests" begin
    for warning in [true, false]
        @testset "Variable bounds" begin
            model = Xpress.Optimizer(OUTPUTLOG = 0, DEFAULTALG = 3, PRESOLVE = 0)
            MOI.set(model, MOI.RawParameter("MOI_WARNINGS"), warning)
            x = MOI.add_variable(model)
            c1 = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.GreaterThan(2.0))
            c2 = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.LessThan(1.0))

            # Getting the results before the conflict refiner has been called must return an error.
            @test MOI.get(model, MOI.ConflictStatus()) == MOI.COMPUTE_CONFLICT_NOT_CALLED
            @test_throws ErrorException MOI.get(model, MOI.ConstraintConflictStatus(), c1)

            # Once it's called, no problem.
            MOI.compute_conflict!(model)
            @test MOI.get(model, MOI.ConflictStatus()) == MOI.CONFLICT_FOUND
            @test MOI.get(model, MOI.ConstraintConflictStatus(), c1) == MOI.IN_CONFLICT
            @test MOI.get(model, MOI.ConstraintConflictStatus(), c2) == MOI.IN_CONFLICT
        end
    end

    @testset "Two conflicting constraints (GreaterThan, LessThan)" begin
        model = Xpress.Optimizer(OUTPUTLOG = 0)
        x = MOI.add_variable(model)
        y = MOI.add_variable(model)
        b1 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 0.0), MOI.GreaterThan(0.0))
        b2 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [y]), 0.0), MOI.GreaterThan(0.0))
        cf1 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0)
        c1 = MOI.add_constraint(model, cf1, MOI.LessThan(-1.0))
        cf2 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, -1.0], [x, y]), 0.0)
        c2 = MOI.add_constraint(model, cf2, MOI.GreaterThan(1.0))

        # Getting the results before the conflict refiner has been called must return an error.
        @test MOI.get(model, MOI.ConflictStatus()) == MOI.COMPUTE_CONFLICT_NOT_CALLED
        @test_throws ErrorException MOI.get(model, MOI.ConstraintConflictStatus(), c1)

        # Once it's called, no problem.
        # Two possible IISes: b1, b2, c1 OR b2, c1, c2
        MOI.compute_conflict!(model)
        @test MOI.get(model, MOI.ConflictStatus()) == MOI.CONFLICT_FOUND
        @test MOI.get(model, MOI.ConstraintConflictStatus(), b1) in [MOI.IN_CONFLICT, MOI.NOT_IN_CONFLICT]
        @test MOI.get(model, MOI.ConstraintConflictStatus(), b2) == MOI.IN_CONFLICT
        @test MOI.get(model, MOI.ConstraintConflictStatus(), c1) == MOI.IN_CONFLICT
        @test MOI.get(model, MOI.ConstraintConflictStatus(), c2) in [MOI.IN_CONFLICT, MOI.NOT_IN_CONFLICT]
    end

    @testset "Two conflicting constraints (EqualTo)" begin
        model = Xpress.Optimizer(OUTPUTLOG = 0)
        x = MOI.add_variable(model)
        y = MOI.add_variable(model)
        b1 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 0.0), MOI.GreaterThan(0.0))
        b2 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [y]), 0.0), MOI.GreaterThan(0.0))
        cf1 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0)
        c1 = MOI.add_constraint(model, cf1, MOI.EqualTo(-1.0))
        cf2 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, -1.0], [x, y]), 0.0)
        c2 = MOI.add_constraint(model, cf2, MOI.GreaterThan(1.0))

        # Getting the results before the conflict refiner has been called must return an error.
        @test MOI.get(model, MOI.ConflictStatus()) == MOI.COMPUTE_CONFLICT_NOT_CALLED
        @test_throws ErrorException MOI.get(model, MOI.ConstraintConflictStatus(), c1)

        # Once it's called, no problem.
        # Two possible IISes: b1, b2, c1 OR b2, c1, c2
        MOI.compute_conflict!(model)
        @test MOI.get(model, MOI.ConflictStatus()) == MOI.CONFLICT_FOUND
        @test MOI.get(model, MOI.ConstraintConflictStatus(), b1) in [MOI.IN_CONFLICT, MOI.NOT_IN_CONFLICT]
        @test MOI.get(model, MOI.ConstraintConflictStatus(), b2) == MOI.IN_CONFLICT
        @test MOI.get(model, MOI.ConstraintConflictStatus(), c1) == MOI.IN_CONFLICT
        @test MOI.get(model, MOI.ConstraintConflictStatus(), c2) in [MOI.IN_CONFLICT, MOI.NOT_IN_CONFLICT]
    end

    @testset "Variables outside conflict" begin
        model = Xpress.Optimizer(OUTPUTLOG = 0)
        x = MOI.add_variable(model)
        y = MOI.add_variable(model)
        z = MOI.add_variable(model)
        b1 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 0.0), MOI.GreaterThan(0.0))
        b2 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [y]), 0.0), MOI.GreaterThan(0.0))
        b3 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [z]), 0.0), MOI.GreaterThan(0.0))
        cf1 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0)
        c1 = MOI.add_constraint(model, cf1, MOI.LessThan(-1.0))
        cf2 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, -1.0, 1.0], [x, y, z]), 0.0)
        c2 = MOI.add_constraint(model, cf2, MOI.GreaterThan(1.0))

        # Getting the results before the conflict refiner has been called must return an error.
        @test MOI.get(model, MOI.ConflictStatus()) == MOI.COMPUTE_CONFLICT_NOT_CALLED
        @test_throws ErrorException MOI.get(model, MOI.ConstraintConflictStatus(), c1)

        # Once it's called, no problem.
        MOI.compute_conflict!(model)
        @test MOI.get(model, MOI.ConflictStatus()) == MOI.CONFLICT_FOUND
        @test MOI.get(model, MOI.ConstraintConflictStatus(), b1) == MOI.IN_CONFLICT
        @test MOI.get(model, MOI.ConstraintConflictStatus(), b2) == MOI.IN_CONFLICT
        @test MOI.get(model, MOI.ConstraintConflictStatus(), b3) == MOI.NOT_IN_CONFLICT
        @test MOI.get(model, MOI.ConstraintConflictStatus(), c1) == MOI.IN_CONFLICT
        @test MOI.get(model, MOI.ConstraintConflictStatus(), c2) == MOI.NOT_IN_CONFLICT
    end

    @testset "No conflict" begin
        model = Xpress.Optimizer(OUTPUTLOG = 0)
        x = MOI.add_variable(model)
        c1 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 0.0), MOI.GreaterThan(1.0))
        c2 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 0.0), MOI.LessThan(2.0))

        # Getting the results before the conflict refiner has been called must return an error.
        @test MOI.get(model, MOI.ConflictStatus()) == MOI.COMPUTE_CONFLICT_NOT_CALLED
        @test_throws ErrorException MOI.get(model, MOI.ConstraintConflictStatus(), c1)

        # Once it's called, no problem.
        MOI.compute_conflict!(model)
        @test MOI.get(model, MOI.ConflictStatus()) == MOI.NO_CONFLICT_FOUND
        @test MOI.get(model, MOI.ConstraintConflictStatus(), c1) == MOI.NOT_IN_CONFLICT
        @test MOI.get(model, MOI.ConstraintConflictStatus(), c2) == MOI.NOT_IN_CONFLICT
    end
end

@testset "test_farkas_dual_min" begin
    model = Xpress.Optimizer(OUTPUTLOG = 0, PRESOLVE = 0)
    x = MOI.add_variables(model, 2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.SingleVariable}(),
        MOI.SingleVariable(x[1]),
    )
    clb = MOI.add_constraint.(
        model, MOI.SingleVariable.(x), MOI.GreaterThan(0.0)
    )
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 1.0], x), 0.0),
        MOI.LessThan(-1.0),
    )
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
    @test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
    clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb)
    c_dual = MOI.get(model, MOI.ConstraintDual(), c)
    @test clb_dual[1] > 1e-6
    @test clb_dual[2] > 1e-6
    @test c_dual[1] < -1e-6
    @test clb_dual[1] ≈ -2 * c_dual atol = 1e-6
    @test clb_dual[2] ≈ -c_dual atol = 1e-6
end

@testset "test_farkas_dual_min_interval" begin
    model = Xpress.Optimizer(OUTPUTLOG = 0, PRESOLVE = 0)
    x = MOI.add_variables(model, 2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.SingleVariable}(),
        MOI.SingleVariable(x[1]),
    )
    clb = MOI.add_constraint.(
        model, MOI.SingleVariable.(x), MOI.Interval(0.0, 10.0)
    )
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 1.0], x), 0.0),
        MOI.LessThan(-1.0),
    )
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
    @test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
    clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb)
    c_dual = MOI.get(model, MOI.ConstraintDual(), c)
    @test clb_dual[1] > 1e-6
    @test clb_dual[2] > 1e-6
    @test c_dual[1] < -1e-6
    @test clb_dual[1] ≈ -2 * c_dual atol = 1e-6
    @test clb_dual[2] ≈ -c_dual atol = 1e-6
end

@testset "test_farkas_dual_min_equalto" begin
    model = Xpress.Optimizer(OUTPUTLOG = 0, PRESOLVE = 0)
    x = MOI.add_variables(model, 2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.SingleVariable}(),
        MOI.SingleVariable(x[1]),
    )
    clb = MOI.add_constraint.(model, MOI.SingleVariable.(x), MOI.EqualTo(0.0))
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 1.0], x), 0.0),
        MOI.LessThan(-1.0),
    )
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
    @test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
    clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb)
    c_dual = MOI.get(model, MOI.ConstraintDual(), c)
    @test clb_dual[1] > 1e-6
    @test clb_dual[2] > 1e-6
    @test c_dual[1] < -1e-6
    @test clb_dual[1] ≈ -2 * c_dual atol = 1e-6
    @test clb_dual[2] ≈ -c_dual atol = 1e-6
end

@testset "test_farkas_dual_min_ii" begin
    model = Xpress.Optimizer(OUTPUTLOG = 0, PRESOLVE = 0)
    x = MOI.add_variables(model, 2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(-1.0, x[1])], 0.0),
    )
    clb = MOI.add_constraint.(
        model, MOI.SingleVariable.(x), MOI.LessThan(0.0)
    )
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-2.0, -1.0], x), 0.0),
        MOI.LessThan(-1.0),
    )
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
    @test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
    clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb)
    c_dual = MOI.get(model, MOI.ConstraintDual(), c)
    @test clb_dual[1] < -1e-6
    @test clb_dual[2] < -1e-6
    @test c_dual[1] < -1e-6
    @test clb_dual[1] ≈ 2 * c_dual atol = 1e-6
    @test clb_dual[2] ≈ c_dual atol = 1e-6
end

@testset "test_farkas_dual_max" begin
    model = Xpress.Optimizer(OUTPUTLOG = 0, PRESOLVE = 0)
    x = MOI.add_variables(model, 2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.SingleVariable}(),
        MOI.SingleVariable(x[1]),
    )
    clb = MOI.add_constraint.(
        model, MOI.SingleVariable.(x), MOI.GreaterThan(0.0)
    )
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 1.0], x), 0.0),
        MOI.LessThan(-1.0),
    )
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
    @test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
    clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb)
    c_dual = MOI.get(model, MOI.ConstraintDual(), c)
    @test clb_dual[1] > 1e-6
    @test clb_dual[2] > 1e-6
    @test c_dual[1] < -1e-6
    @test clb_dual[1] ≈ -2 * c_dual atol = 1e-6
    @test clb_dual[2] ≈ -c_dual atol = 1e-6
end

@testset "test_farkas_dual_max_ii" begin
    model = Xpress.Optimizer(OUTPUTLOG = 0, PRESOLVE = 0)
    x = MOI.add_variables(model, 2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(-1.0, x[1])], 0.0),
    )
    clb = MOI.add_constraint.(
        model, MOI.SingleVariable.(x), MOI.LessThan(0.0)
    )
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-2.0, -1.0], x), 0.0),
        MOI.LessThan(-1.0),
    )
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
    @test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
    clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb)
    c_dual = MOI.get(model, MOI.ConstraintDual(), c)
    @test clb_dual[1] < -1e-6
    @test clb_dual[2] < -1e-6
    @test c_dual[1] < -1e-6
    @test clb_dual[1] ≈ 2 * c_dual atol = 1e-6
    @test clb_dual[2] ≈ c_dual atol = 1e-6
end

@testset "Delete equality constraint in binary variable" begin
    atol = rtol = 1e-6
    model = Xpress.Optimizer(OUTPUTLOG = 0)
    # an example on mixed integer programming
    #
    #   maximize 1.1x + 2 y + 5 z
    #
    #   s.t.  x + y + z <= 10
    #         x + 2 y + z <= 15
    #
    #         x is continuous: 0 <= x <= 5
    #         y is integer: 0 <= y <= 10
    #         z is binary
    v = MOI.add_variables(model, 3)

    cf = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0, 1.0], v), 0.0)
    c = MOI.add_constraint(model, cf, MOI.LessThan(10.0))
    cf2 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 2.0, 1.0], v), 0.0)
    c2 = MOI.add_constraint(model, cf2, MOI.LessThan(15.0))

    vc1 = MOI.add_constraint(
        model,
        MOI.SingleVariable(v[1]),
        MOI.Interval(0.0, 5.0),
    )
    vc2 = MOI.add_constraint(
        model,
        MOI.SingleVariable(v[2]),
        MOI.Interval(0.0, 10.0),
    )
    vc3 = MOI.add_constraint(model, MOI.SingleVariable(v[2]), MOI.Integer())
    vc4 = MOI.add_constraint(model, MOI.SingleVariable(v[3]), MOI.ZeroOne())

    objf =
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.1, 2.0, 5.0], v), 0.0)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        objf,
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    MOI.optimize!(model)

    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(model, MOI.ResultCount()) >= 1
    @test MOI.get(model, MOI.PrimalStatus()) in
            [MOI.FEASIBLE_POINT, MOI.NEARLY_FEASIBLE_POINT]
    @test MOI.get(model, MOI.ObjectiveValue()) ≈ 19.4 atol = atol rtol =
        rtol
    @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [4, 5, 1] atol = atol rtol =
        rtol
    @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 10 atol = atol rtol =
        rtol
    @test MOI.get(model, MOI.ConstraintPrimal(), c2) ≈ 15 atol = atol rtol =
        rtol
    @test MOI.get(model, MOI.ObjectiveBound()) >= 19.4 - atol

    z_value = MOI.get(model, MOI.VariablePrimal(), v)[3]

    # Relax binary bounds of z variable
    MOI.delete(model, vc4)

    # Fix z to its optimal value
    vc5 = MOI.add_constraint(model, MOI.SingleVariable(v[3]), MOI.EqualTo(z_value))

    MOI.optimize!(model)

    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(model, MOI.ResultCount()) >= 1
    @test MOI.get(model, MOI.PrimalStatus()) in
            [MOI.FEASIBLE_POINT, MOI.NEARLY_FEASIBLE_POINT]
    @test MOI.get(model, MOI.ObjectiveValue()) ≈ 19.4 atol = atol rtol =
        rtol
    @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [4, 5, 1] atol = atol rtol =
        rtol
    @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 10 atol = atol rtol =
        rtol
    @test MOI.get(model, MOI.ConstraintPrimal(), c2) ≈ 15 atol = atol rtol =
        rtol
    @test MOI.get(model, MOI.ObjectiveBound()) >= 19.4 - atol

    # Add back binary bounds to z
    vc4 = MOI.add_constraint(model, MOI.SingleVariable(v[3]), MOI.ZeroOne())

    # Remove equality constraint
    MOI.delete(model, vc5)

    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(model, MOI.ResultCount()) >= 1
    @test MOI.get(model, MOI.PrimalStatus()) in
            [MOI.FEASIBLE_POINT, MOI.NEARLY_FEASIBLE_POINT]
    @test MOI.get(model, MOI.ObjectiveValue()) ≈ 19.4 atol = atol rtol =
        rtol
    @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [4, 5, 1] atol = atol rtol =
        rtol
    @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 10 atol = atol rtol =
        rtol
    @test MOI.get(model, MOI.ConstraintPrimal(), c2) ≈ 15 atol = atol rtol =
        rtol
    @test MOI.get(model, MOI.ObjectiveBound()) >= 19.4 - atol
end

@testset "Binary Variables Infeasibility" begin
    atol = rtol = 1e-6
    model = Xpress.Optimizer(OUTPUTLOG = 0)
    v = MOI.add_variable(model)

    infeas_err = ErrorException("The problem is infeasible")
    vc1 = MOI.add_constraint(model, MOI.SingleVariable(v), MOI.ZeroOne())
    @test_throws infeas_err vc2 = MOI.add_constraint(
        model,
        MOI.SingleVariable(v),
        MOI.GreaterThan(2.0),
    )
    @test_throws infeas_err vc3 = MOI.add_constraint(
        model,
        MOI.SingleVariable(v),
        MOI.LessThan(-1.0),
        )
    @test_throws infeas_err vc4 = MOI.add_constraint(
        model,
        MOI.SingleVariable(v),
        MOI.Interval(-1.0,-0.5),
    )
    @test_throws infeas_err vc5 = MOI.add_constraint(
        model,
        MOI.SingleVariable(v),
        MOI.EqualTo(2.0),
    )
end

@testset "MIP Start testing" begin
    # The idea of the test is the following:
    # * Give Xpress a problem it can solve quickly but not simple enough to be
    #   solved at the root node when heuristics and cuts are disabled.
    # * Assert the solution was not found at the root node (query attribute
    #   MIPSOLNODE, it should be higher than one).
    # * Limit the model to searching the root node.
    # * Run again without the MIP-start, so the solution found is worse.
    #   (This is mostly done to be sure no status since the optimal solve
    #   was kept.)
    # * MIP-start the model with the optimal solution.
    # * If the solution found is the MIP-started solution then XPRESS has
    #   used the MIP-start (we checked XPRESS could not solve it using
    #   a single node).
    atol = rtol = 1e-6
    weight = Float64[
        677, 460, 752, 852, 580, 116, 457, 121, 454, 870,
        443, 196, 411, 539, 348, 187, 771, 127, 338, 527,
        932, 961,  48, 135, 367, 998, 363, 153, 921, 578,
        311, 560, 293, 258, 474, 884, 162, 136, 479, 289,
        813, 139, 795, 825, 945, 750, 462, 659, 270,  34,
        758, 865, 238, 367, 444, 116,  69, 894, 584,  96,
         29, 199, 712, 703, 856, 692, 396, 409, 603, 632,
        479, 848, 822, 248, 424, 978, 738, 655, 210, 173,
        731, 100, 889, 195, 245, 329, 446,  47, 235,  25,
        254, 150, 520, 665, 391, 907, 123, 826, 959, 176
    ]
    profit = weight .- (5.0,)
    capacity = 10000.0
    model = Xpress.Optimizer(
        OUTPUTLOG    = 0,
        HEURSTRATEGY = 0,
        CUTSTRATEGY  = 0,
        PRESOLVE     = 0,
        MIPPRESOLVE  = 0,
        PRESOLVEOPS  = 0,
    )
    # The variables: x[1:100], Bin
    x, _ = MOI.add_constrained_variables(model, fill(MOI.ZeroOne(), 100))
    # The objective function: maximize sum(profit' * x)
    objf = MOI.ScalarAffineFunction(
        MOI.ScalarAffineTerm.(profit, x), 0.0
    )
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        objf,
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    # The capacity constraint: sum(weight' * x) <= capacity
    MOI.add_constraint(
        model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(
            weight, x), 0.0
        ),
        MOI.LessThan(capacity)
    )

    # FIRST RUN: get the optimal value. Check that the optimal
    # value was not discovered at the root node.
    MOI.optimize!(model)

    solution = MOI.get(model, MOI.VariablePrimal(), x)
    computed_obj_value = profit' * solution
    obtained_obj_value = MOI.get(model, MOI.ObjectiveValue()) :: Float64
    @test isapprox(
        obtained_obj_value, computed_obj_value; rtol = rtol, atol = atol
    )
    @test isapprox(9945.0, computed_obj_value; rtol = rtol, atol = atol)

    node_solution_was_found = MOI.get(model, MOI.RawParameter("MIPSOLNODE"))

    @test node_solution_was_found > 1

    # SECOND RUN: run without MIP-start and only searching the first node.
    # Should give a worse solution than the previous one.
    MOI.set(model, MOI.RawParameter("MAXNODE"), 1)

    MOI.optimize!(model)

    # One node may not be enough to even get any solution.
    if MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
        solution2 = MOI.get(model, MOI.VariablePrimal(), x)
        computed_obj_value2 = profit' * solution2
        obtained_obj_value2 = MOI.get(model, MOI.ObjectiveValue()) :: Float64
        @test isapprox(
            obtained_obj_value2, computed_obj_value2; rtol = rtol, atol = atol
        )
        # There should be at least one unit of difference.
        @test obtained_obj_value1 > obtained_obj_value2 + 0.5
    end

    # THIRD RUN: run with MIP-start and searching only the root node.
    # Should find the optimal solution impossible to get in one node.
    MOI.set(model, MOI.VariablePrimalStart(), x, solution)

    # This postsolve is necessary because of an unrelated bug. Apparently,
    # if Xpress is stopped because MAXNODE and has no feasible solution,
    # then MOI.optimize! will not call postsolve over it, however calling
    # postsolve is necessary otherwise a new call to MOI.optimize will
    # trigger error 707 ("707 Error: Function cannot be called during the
    # global search, except in callbacks.").
    Xpress.postsolve(model.inner)

    MOI.optimize!(model)

    solution3 = MOI.get(model, MOI.VariablePrimal(), x)
    computed_obj_value3 = profit' * solution3
    obtained_obj_value3 = MOI.get(model, MOI.ObjectiveValue()) :: Float64
    @test isapprox(
        obtained_obj_value3, computed_obj_value3; rtol = rtol, atol = atol
    )
    @test isapprox(9945.0, computed_obj_value3; rtol = rtol, atol = atol)
end


text:
"""
    second-order cone
        assumes: z ∈ R³
"""
function κ_so(z)
    z1 = view(z, 2:3)
    z0 = z[1]

    z_proj = zero(z)
    status = false

    if norm(z1) <= z0
        z_proj = copy(z)
        if norm(z1) < z0
            status = true
        end
        # status = true
    elseif norm(z1) <= -z0
        z_proj = zero(z)
    else
        a = 0.5 * (1.0 + z0 / norm(z1))
        z_proj[1:end-1] = a * z1
        z_proj[end] = a * norm(z1)
    end

    return z_proj, status
end

function cone_product(z, s)
    [z' * s; z[1] * view(s, 2:3) + s[1] * view(z, 2:3)]
end


text:
import Term.consoles: console_height, console_width, console, err_console

println("\nTesting tprint, stdout temporarily disabled")

@suppress_out begin
    @testset "\e[34mTPRINT" begin
        @test_nowarn tprint("string")
        @test_nowarn tprint("[red]adasd[/red]")
        @test_nowarn tprint("[blue on_green]adasda")
        @test_nowarn tprint("[red]dadas[green]insdai[/green]outssdrse[blue]fsfsf[/blue]")

        @test_nowarn tprint(Panel("test"))
        @test_nowarn tprint(TextBox("test"))

        @test_nowarn tprint(1)
        @test_nowarn tprint(:x)

        @test_nowarn tprint(1, Panel("test"), "test")
    end
end

@testset "\e[34mCONSOLE" begin
    @test console_height() == displaysize(stdout)[1]
    @test console_height(stdout) == displaysize(stdout)[1]

    @test console_width() == displaysize(stdout)[2]
    @test console_width(stdout) == displaysize(stdout)[2]
end


text:
#cd("/home/user/github/TOF-Tracer2/")

include("$(pwd())/startupAPiTOF.jl")

fp = "./ExampleFiles/APiTOFDATA/" # All files in this path will be processed
filefilterRegexp = r"\.h5$"
rf = "./ExampleFiles/APiTOFDATA/Data_09_46_20_04mm.h5"  # The mass scale from this file defines the mass scale of all

masslist = MasslistFunctions.loadMasslist("./ExampleFiles/MASSLISTS/exampleMasslist.csv")
cr = [37.028406 55.038971 282.118343] # pos
#cr = [62 125] # neg

# alternatively: use an auto generated masslist
# masslistMasses, masslistElements, masslistElementsMasses, masslistCompositions = createMassList(C=0:20, O=0:20, N=0:1, allowRadicals=false) #
(masslistMasses, masslistElements, masslistElementsMasses, masslistCompositions) = masslist
s = (masslistMasses.>17) .& ( masslistMasses.<900)
masslistMasses = masslistMasses[s]
masslistCompositions = masslistCompositions[s,:]

####################### END OF SETTINGS ###############################################################

####################### Processing sequence ###########################################################

correctMassScaleAndExtractSumSpec(
    fp,
    masslistMasses,
    masslistElements,
    masslistElementsMasses,
    masslistCompositions,
    rf,
    cr,
    filefilterRegexp=filefilterRegexp,
    onlyUseAverages = true,
    plotControlMass = true,
    testRangeStart = 54.5, # the mass shift of this region will be shown if plot control mass is set true. Should not be part of calibRegions
    testRangeEnd = 55.5,
    recalibInterval = 900,
    resolution = 5500,
    firstNFiles= 1,
    lastNFiles = 0,
    UMRpeaks = true, # create Unit Mass Resolution peaks
    UMRmasses = length(s),
    UMRmassLow = 0.2,
    UMRmassHigh = 0.4,
    umrBinWidth =100
    #binWidth = 10
    )
#
baselineAndPeakshape(
    fp,
    peakshapeRegions=4,
    peakshapeRegionStretch=1.0,
    peakshapeQuantileValue = 0.7,
    peakfindingNoiseThresholdValue = 5,
    peakWindowWidth = 100,
    peakfindingSignalLimit=0.5
    )
    #
#
mtrx = deconvolute(
    fp,
    calcTransposed = false
    )
#


text:
abstract type InitializationScheme end
struct TaylorModeInit <: InitializationScheme end
Base.@kwdef struct ClassicSolverInit{ALG} <: InitializationScheme
    alg::ALG = Tsit5()
    init_on_du::Bool = false
end

function initial_update!(integ, cache)
    return initial_update!(integ, cache, integ.alg.initialization)
end

########################################################################
# Some utilities below
"""Quick and dirty wrapper to make OOP functions IIP"""
function oop_to_iip(f)
    function f!(du, u, p, t)
        du .= f(u, p, t)
        return nothing
    end
    return f!
end

"""Basically an Kalman update"""
function condition_on!(
    x::SRGaussian,
    H::AbstractMatrix,
    data::AbstractVector,
    meascache,
    Kcache,
    covcache,
    Mcache,
)
    z, S = meascache

    _matmul!(z, H, x.μ)
    X_A_Xt!(S, x.Σ, H)
    @assert isdiag(S)
    S_diag = diag(S)
    if any(iszero.(S_diag)) # could happen with a singular mass-matrix
        S_diag .+= 1e-20
    end

    _matmul!(Kcache, x.Σ.mat, H')
    K = Kcache ./= S_diag'

    _matmul!(x.μ, K, data - z, 1.0, 1.0)
    # x.μ .+= K*(data - z)

    D = length(x.μ)
    mul!(Mcache, K, H, -1.0, 0.0)
    @inbounds @simd ivdep for i in 1:D
        Mcache[i, i] += 1
    end
    X_A_Xt!(covcache, x.Σ, Mcache)
    copy!(x.Σ, covcache)
    return nothing
end


text:
__precompile__()

module IOIndents

import Base: convert, show, pipe_reader, pipe_writer, lock, unlock, write,
             getindex, in, haskey, get, print

export IOIndent, Indent, Dedent, Align, Dealign, indent_string!, alignment_char!

struct Indent end
struct Dedent end

struct Align end
struct Dealign end

Base.show(io::IO, ::Union{Indent, Dedent, Align, Dealign}) = nothing

mutable struct IOIndent{IO_t <: IO} <: Base.AbstractPipe
    io::IO_t
    indent_level::Int
    aligns::Vector{Int}
    offset::Int
    indented_line::Bool
    indent_str::String
    align_char::Char

    function IOIndent{IO_t}(io::IO_t, indent_level::Int, aligns::Vector{Int},
                            offset::Int, indented_line::Bool, indent_str::String,
                            align_char::Char) where IO_t <: IO
        assert(!(IO_t <: IOIndent))
        return new(io, indent_level, aligns, offset, indented_line, indent_str, align_char)
    end
end

indent_string!(io::IOIndent, str::String) = (io.indent_str = str; io)
alignment_char!(io::IOIndent, chr::Char) = (io.align_char = chr; io)

IOIndent(io::IO) = IOIndent{typeof(io)}(io, 0, Int[], 0, false, "    ", ' ')

convert(::Type{IOIndent}, io::IOIndent) = io

in(key_value::Pair, io::IOIndent) = in(key_value, io.io, ===)
haskey(io::IOIndent, key) = haskey(io.io, key)
getindex(io::IOIndent, key) = getindex(io.io, key)
get(io::IOIndent, key, default) = get(io.io, key, default)

function show(_io::IO, io::IOIndent)
    ioi = IOIndent(_io)
    print(ioi, "IOIndent:", Indent())
    print(ioi, "\nIO: "); show(ioi, io.io)
    print(ioi, "\nIndent string: \"", io.indent_str, "\"")
    print(ioi, "\nAlign char: \"", io.align_char, "\"")
    print(ioi, "\nIndent: ", io.indent_level)
    print(ioi, "\nAligns: ", join(io.aligns, ","))
end

pipe_reader(io::IOIndent) = io.io
pipe_writer(io::IOIndent) = io.io
lock(io::IOIndent) = lock(io.io)
unlock(io::IOIndent) = unlock(io.io)

displaysize(io::IOContext) = displaysize(io.io)

write(io::IOIndent, ::Indent) = (io.indent_level += 1; 0)
print(io::IOIndent, ::Indent) = write(io, Indent())
write(io::IOIndent, ::Dedent) = (io.indent_level = max(0, io.indent_level - 1); 0)
print(io::IOIndent, ::Dedent) = write(io, Dedent())

_align_length(io) = length(io.aligns) == 0 ? 0 : io.aligns[end]
write(io::IOIndent, ::Align) = (push!(io.aligns, io.offset + _align_length(io)); 0)
print(io::IOIndent, ::Align) = write(io, Align())
write(io::IOIndent, ::Dealign) = (length(io.aligns) > 0 && pop!(io.aligns); 0)
print(io::IOIndent, ::Dealign) = write(io, Dealign())

write_indent(io::IOIndent) = write(io.io, io.indent_str^io.indent_level)
write_offset(io::IOIndent) = write(io.io, string(io.align_char)^_align_length(io))

write(io::IOIndent, chr::Char) = write(io, string(chr)) # Need to catch writing a '\n'
function write(io::IOIndent, str::String)
    written = 0
    for (i, line) in enumerate(split(str, "\n"))
        if i != 1
            written += write(io.io, "\n")
            io.indented_line = false
            io.offset = 0
        end
        if !io.indented_line && !isempty(line)
            written += write_indent(io)
            written += write_offset(io)
            io.indented_line = true
        end
        written += write(io.io, line)
        io.offset += strwidth(line)
    end
    return written
end

end # module


text:
function new_cdatanode(xdoc::XMLDocument, txt::ASCIIString)
        p = ccall(xmlNewCDataBlock, Xptr, (Xptr, Xstr, Cint), xdoc.ptr, txt, length(txt)+1)
        XMLNode(p)
end

add_cdata(xdoc::XMLDocument, x::XMLElement, txt::ASCIIString) = add_child(x, new_cdatanode(xdoc,txt))


text:
export RewardOverriddenEnv

struct RewardOverriddenEnv{F,E<:AbstractEnv} <: AbstractEnvWrapper
    env::E
    f::F
end

RewardOverriddenEnv(f) = env -> RewardOverriddenEnv(f, env)

RLBase.reward(env::RewardOverriddenEnv, args...; kwargs...) =
    env.f(reward(env.env, args...; kwargs...))


text:
struct Point{T}
  x::T
  y::T
end

struct Point3{T}
  x::T
  y::T
  z::T
end

struct Size{T}
  width::T
  height::T
end

struct Rect{T}
  x::T
  y::T
  width::T
  height::T
end

struct RotatedRect
    center::Point{Float32}
    size::Size{Float32}
    angle::Float32
end

struct Range
    start::Int32
    end_::Int32
end

struct TermCriteria
    type::Int32
    maxCount::Int32
    epsilon::Float64
end

struct cvComplex{T}
  re::T
  im::T
end


text:
using SparseArrays, DiffEqOperators, LinearAlgebra, Random,
      Test, BandedMatrices, FillArrays, LazyArrays, BlockBandedMatrices

# This test file tests for the correctness of higher dimensional concretization.
# The tests verify that multiplication in the concretized case agrees with the matrix-free
# multiplication

@testset "First Dimension" begin

      # Test that even when we have a vector, the concretizations using the higher dimension dispatch still function
      # correctly
      M = rand(22)

      L1 = CenteredDifference(1,2,1.0,20)
      L2 = CenteredDifference(1,2,1.0,20)
      L3 = CenteredDifference(3,3,1.0,20)

      @test L1*M ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*M ≈ BandedMatrix(L1, size(M))*M ≈ BandedBlockBandedMatrix(L1,size(M))*M
      @test L2*M ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*M ≈ BandedMatrix(L2, size(M))*M ≈ BandedBlockBandedMatrix(L2,size(M))*M
      @test L3*M ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*M ≈ BandedMatrix(L3, size(M))*M ≈ BandedBlockBandedMatrix(L3,size(M))*M

      M = rand(22,2,2,2,2)

      @test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L1,size(M))*vec(M)
      @test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L2,size(M))*vec(M)
      @test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L3,size(M))*vec(M)

end

@testset "Second Dimension" begin

      M = rand(2,22,2)

      L1 = CenteredDifference{2}(1,2,1.0,20)
      L2 = CenteredDifference{2}(1,2,1.0,20)
      L3 = CenteredDifference{2}(3,3,1.0,20)

      @test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L1,size(M))*vec(M)
      @test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L2,size(M))*vec(M)
      @test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L3,size(M))*vec(M)

      M = rand(2,22,2,2)

      @test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L1,size(M))*vec(M)
      @test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L2,size(M))*vec(M)
      @test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L3,size(M))*vec(M)

      M = rand(2,22,2,2,3)

      @test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L1,size(M))*vec(M)
      @test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L2,size(M))*vec(M)
      @test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L3,size(M))*vec(M)

end

@testset "Third Dimension" begin

      M = rand(3,2,22)

      L1 = CenteredDifference{3}(1,2,1.0,20)
      L2 = CenteredDifference{3}(1,2,1.0,20)
      L3 = CenteredDifference{3}(3,3,1.0,20)

      @test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L1,size(M))*vec(M)
      @test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L2,size(M))*vec(M)
      @test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L3,size(M))*vec(M)

      M = rand(3,2,22,2)

      @test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L1,size(M))*vec(M)
      @test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L2,size(M))*vec(M)
      @test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L3,size(M))*vec(M)

      M = rand(3,2,22,2,3)

      @test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L1,size(M))*vec(M)
      @test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L2,size(M))*vec(M)
      @test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L3,size(M))*vec(M)

end

@testset "Fifth Dimension" begin

      M = rand(3,2,3,2,22)

      L1 = CenteredDifference{5}(1,2,1.0,20)
      L2 = CenteredDifference{5}(1,2,1.0,20)
      L3 = CenteredDifference{5}(3,3,1.0,20)

      @test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L1,size(M))*vec(M)
      @test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L2,size(M))*vec(M)
      @test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L3,size(M))*vec(M)

      M = rand(3,2,3,2,22,3)

      @test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L1,size(M))*vec(M)
      @test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L2,size(M))*vec(M)
      @test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L3,size(M))*vec(M)

      M = rand(2,2,2,2,22,3,2)

      @test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L1,size(M))*vec(M)
      @test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L2,size(M))*vec(M)
      @test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L3,size(M))*vec(M)

end

@testset "Biased Upwind : Positive Wind with offside, regular grid" begin

      # Test that Biased Upwind with maximum offside concretizes as CenteredDifference of same order
      M = rand(22)
      P = PeriodicBC(Float64)*M[2:21]

      L1 = CenteredDifference(1,2,1.0,20)
      L2 = CenteredDifference(1,4,1.0,20)
      L3 = CenteredDifference(1,6,1.0,20)
      L4 = CenteredDifference(2,3,1.0,20)
      L5 = CenteredDifference(3,4,1.0,20)

      K1 = UpwindDifference(1,2,1.0,20,offside=1)
      K2 = UpwindDifference(1,4,1.0,20,offside=2)
      K3 = UpwindDifference(1,6,1.0,20,offside=3)
      K4 = UpwindDifference(2,3,1.0,20,offside=2)
      K5 = UpwindDifference(3,4,1.0,20,offside=3)

      @test L1*M ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*M ≈ BandedMatrix(L1, size(M))*M ≈ BandedBlockBandedMatrix(L1,size(M))*M ≈ K1*M ≈ Array(K1, size(M))*vec(M) ≈ sparse(K1,size(M))*M ≈ BandedMatrix(K1, size(M))*M ≈ BandedBlockBandedMatrix(K1,size(M))*M
      @test L2*M ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*M ≈ BandedMatrix(L2, size(M))*M ≈ BandedBlockBandedMatrix(L2,size(M))*M ≈ K2*M ≈ Array(K2, size(M))*vec(M) ≈ sparse(K2,size(M))*M ≈ BandedMatrix(K2, size(M))*M ≈ BandedBlockBandedMatrix(K2,size(M))*M
      @test L2*P ≈ Array(L2, size(P))*vec(P) ≈ sparse(L2,size(P))*P ≈ BandedMatrix(L2, size(P))*P ≈ BandedBlockBandedMatrix(L2,size(P))*P ≈ K2*P ≈ Array(K2, size(P))*vec(P) ≈ sparse(K2,size(P))*P ≈ BandedMatrix(K2, size(P))*P ≈ BandedBlockBandedMatrix(K2,size(P))*P
      @test L3*M ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*M ≈ BandedMatrix(L3, size(M))*M ≈ BandedBlockBandedMatrix(L3,size(M))*M ≈ K3*M ≈ Array(K3, size(M))*vec(M) ≈ sparse(K3,size(M))*M ≈ BandedMatrix(K3, size(M))*M ≈ BandedBlockBandedMatrix(K3,size(M))*M
      @test L4*M ≈ Array(L4, size(M))*vec(M) ≈ sparse(L4,size(M))*M ≈ BandedMatrix(L4, size(M))*M ≈ BandedBlockBandedMatrix(L4,size(M))*M ≈ K4*M ≈ Array(K4, size(M))*vec(M) ≈ sparse(K4,size(M))*M ≈ BandedMatrix(K4, size(M))*M ≈ BandedBlockBandedMatrix(K4,size(M))*M
      @test L5*M ≈ Array(L5, size(M))*vec(M) ≈ sparse(L5,size(M))*M ≈ BandedMatrix(L5, size(M))*M ≈ BandedBlockBandedMatrix(L5,size(M))*M ≈ K5*M ≈ Array(K5, size(M))*vec(M) ≈ sparse(K5,size(M))*M ≈ BandedMatrix(K5, size(M))*M ≈ BandedBlockBandedMatrix(K5,size(M))*M

      M = rand(22,2,2,2,2)

      @test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L1,size(M))*vec(M) ≈ vec(K1*M) ≈ Array(K1, size(M))*vec(M) ≈ sparse(K1,size(M))*vec(M) ≈ BandedMatrix(K1, size(M))*vec(M) ≈ BandedBlockBandedMatrix(K1,size(M))*vec(M)
      @test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L2,size(M))*vec(M) ≈ vec(K2*M) ≈ Array(K2, size(M))*vec(M) ≈ sparse(K2,size(M))*vec(M) ≈ BandedMatrix(K2, size(M))*vec(M) ≈ BandedBlockBandedMatrix(K2,size(M))*vec(M)
      @test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L3,size(M))*vec(M) ≈ vec(K3*M) ≈ Array(K3, size(M))*vec(M) ≈ sparse(K3,size(M))*vec(M) ≈ BandedMatrix(K3, size(M))*vec(M) ≈ BandedBlockBandedMatrix(K3,size(M))*vec(M)
      @test vec(L4*M) ≈ Array(L4, size(M))*vec(M) ≈ sparse(L4,size(M))*vec(M) ≈ BandedMatrix(L4, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L4,size(M))*vec(M) ≈ vec(K4*M) ≈ Array(K4, size(M))*vec(M) ≈ sparse(K4,size(M))*vec(M) ≈ BandedMatrix(K4, size(M))*vec(M) ≈ BandedBlockBandedMatrix(K4,size(M))*vec(M)
      @test vec(L5*M) ≈ Array(L5, size(M))*vec(M) ≈ sparse(L5,size(M))*vec(M) ≈ BandedMatrix(L5, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L5,size(M))*vec(M) ≈ vec(K5*M) ≈ Array(K5, size(M))*vec(M) ≈ sparse(K5,size(M))*vec(M) ≈ BandedMatrix(K5, size(M))*vec(M) ≈ BandedBlockBandedMatrix(K5,size(M))*vec(M)

end

@testset "Biased Upwind : Negative Wind with offside, regular grid" begin

      M = rand(22)

      L1 = CenteredDifference(1,2,1.0,20,-1)
      L2 = CenteredDifference(1,4,1.0,20,-1)
      L3 = CenteredDifference(1,6,1.0,20,-1)
      L4 = CenteredDifference(2,3,1.0,20,-1)
      L5 = CenteredDifference(3,4,1.0,20,-1)

      K1 = UpwindDifference(1,2,1.0,20,-1,offside=1)
      K2 = UpwindDifference(1,4,1.0,20,-1,offside=2)
      K3 = UpwindDifference(1,6,1.0,20,-1,offside=3)
      K4 = UpwindDifference(2,3,1.0,20,-1,offside=2)
      K5 = UpwindDifference(3,4,1.0,20,-1,offside=3)

      @test L1*M ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*M ≈ BandedMatrix(L1, size(M))*M ≈ BandedBlockBandedMatrix(L1,size(M))*M ≈ K1*M ≈ Array(K1, size(M))*vec(M) ≈ sparse(K1,size(M))*M ≈ BandedMatrix(K1, size(M))*M ≈ BandedBlockBandedMatrix(K1,size(M))*M
      @test L2*M ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*M ≈ BandedMatrix(L2, size(M))*M ≈ BandedBlockBandedMatrix(L2,size(M))*M ≈ K2*M ≈ Array(K2, size(M))*vec(M) ≈ sparse(K2,size(M))*M ≈ BandedMatrix(K2, size(M))*M ≈ BandedBlockBandedMatrix(K2,size(M))*M
      @test L3*M ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*M ≈ BandedMatrix(L3, size(M))*M ≈ BandedBlockBandedMatrix(L3,size(M))*M ≈ K3*M ≈ Array(K3, size(M))*vec(M) ≈ sparse(K3,size(M))*M ≈ BandedMatrix(K3, size(M))*M ≈ BandedBlockBandedMatrix(K3,size(M))*M
      @test L4*M ≈ Array(L4, size(M))*vec(M) ≈ sparse(L4,size(M))*M ≈ BandedMatrix(L4, size(M))*M ≈ BandedBlockBandedMatrix(L4,size(M))*M ≈ K4*M ≈ Array(K4, size(M))*vec(M) ≈ sparse(K4,size(M))*M ≈ BandedMatrix(K4, size(M))*M ≈ BandedBlockBandedMatrix(K4,size(M))*M
      @test L5*M ≈ Array(L5, size(M))*vec(M) ≈ sparse(L5,size(M))*M ≈ BandedMatrix(L5, size(M))*M ≈ BandedBlockBandedMatrix(L5,size(M))*M ≈ K5*M ≈ Array(K5, size(M))*vec(M) ≈ sparse(K5,size(M))*M ≈ BandedMatrix(K5, size(M))*M ≈ BandedBlockBandedMatrix(K5,size(M))*M

      M = rand(22,2,2,2,2)

      @test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L1,size(M))*vec(M) ≈ vec(K1*M) ≈ Array(K1, size(M))*vec(M) ≈ sparse(K1,size(M))*vec(M) ≈ BandedMatrix(K1, size(M))*vec(M) ≈ BandedBlockBandedMatrix(K1,size(M))*vec(M)
      @test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L2,size(M))*vec(M) ≈ vec(K2*M) ≈ Array(K2, size(M))*vec(M) ≈ sparse(K2,size(M))*vec(M) ≈ BandedMatrix(K2, size(M))*vec(M) ≈ BandedBlockBandedMatrix(K2,size(M))*vec(M)
      @test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L3,size(M))*vec(M) ≈ vec(K3*M) ≈ Array(K3, size(M))*vec(M) ≈ sparse(K3,size(M))*vec(M) ≈ BandedMatrix(K3, size(M))*vec(M) ≈ BandedBlockBandedMatrix(K3,size(M))*vec(M)
      @test vec(L4*M) ≈ Array(L4, size(M))*vec(M) ≈ sparse(L4,size(M))*vec(M) ≈ BandedMatrix(L4, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L4,size(M))*vec(M) ≈ vec(K4*M) ≈ Array(K4, size(M))*vec(M) ≈ sparse(K4,size(M))*vec(M) ≈ BandedMatrix(K4, size(M))*vec(M) ≈ BandedBlockBandedMatrix(K4,size(M))*vec(M)
      @test vec(L5*M) ≈ Array(L5, size(M))*vec(M) ≈ sparse(L5,size(M))*vec(M) ≈ BandedMatrix(L5, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L5,size(M))*vec(M) ≈ vec(K5*M) ≈ Array(K5, size(M))*vec(M) ≈ sparse(K5,size(M))*vec(M) ≈ BandedMatrix(K5, size(M))*vec(M) ≈ BandedBlockBandedMatrix(K5,size(M))*vec(M)

end

@testset "Biased Upwind : Positive Wind with offside, irregular grid" begin

      M = rand(7)
      P = PeriodicBC(Float64)*M[2:6]

      L1 = CenteredDifference(1,2,[0.08, 0.02, 0.05, 0.04, 0.07, 0.03],5)
      L2 = CenteredDifference(1,4,[0.08, 0.02, 0.05, 0.04, 0.07, 0.03],5)
      L3 = CenteredDifference(2,3,[0.08, 0.02, 0.05, 0.04, 0.07, 0.03],5)

      K1 = UpwindDifference(1,2,[0.08, 0.02, 0.05, 0.04, 0.07, 0.03],5,offside=1)
      K2 = UpwindDifference(1,4,[0.08, 0.02, 0.05, 0.04, 0.07, 0.03],5,offside=2)
      K3 = UpwindDifference(2,3,[0.08, 0.02, 0.05, 0.04, 0.07, 0.03],5,offside=2)

      @test L1*M ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*M ≈ BandedMatrix(L1, size(M))*M ≈ BandedBlockBandedMatrix(L1,size(M))*M ≈ K1*M ≈ Array(K1, size(M))*vec(M) ≈ sparse(K1,size(M))*M ≈ BandedMatrix(K1, size(M))*M ≈ BandedBlockBandedMatrix(K1,size(M))*M
      @test L2*M ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*M ≈ BandedMatrix(L2, size(M))*M ≈ BandedBlockBandedMatrix(L2,size(M))*M ≈ K2*M ≈ Array(K2, size(M))*vec(M) ≈ sparse(K2,size(M))*M ≈ BandedMatrix(K2, size(M))*M ≈ BandedBlockBandedMatrix(K2,size(M))*M
      @test L2*P ≈ Array(L2, size(P))*vec(P) ≈ sparse(L2,size(P))*P ≈ BandedMatrix(L2, size(P))*P ≈ BandedBlockBandedMatrix(L2,size(P))*P ≈ K2*P ≈ Array(K2, size(P))*vec(P) ≈ sparse(K2,size(P))*P ≈ BandedMatrix(K2, size(P))*P ≈ BandedBlockBandedMatrix(K2,size(P))*P
      @test L3*M ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*M ≈ BandedMatrix(L3, size(M))*M ≈ BandedBlockBandedMatrix(L3,size(M))*M ≈ K3*M ≈ Array(K3, size(M))*vec(M) ≈ sparse(K3,size(M))*M ≈ BandedMatrix(K3, size(M))*M ≈ BandedBlockBandedMatrix(K3,size(M))*M

      M = rand(7,2,2,2)

      @test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L1,size(M))*vec(M) ≈ vec(K1*M) ≈ Array(K1, size(M))*vec(M) ≈ sparse(K1,size(M))*vec(M) ≈ BandedMatrix(K1, size(M))*vec(M) ≈ BandedBlockBandedMatrix(K1,size(M))*vec(M)
      @test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L2,size(M))*vec(M) ≈ vec(K2*M) ≈ Array(K2, size(M))*vec(M) ≈ sparse(K2,size(M))*vec(M) ≈ BandedMatrix(K2, size(M))*vec(M) ≈ BandedBlockBandedMatrix(K2,size(M))*vec(M)
      @test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L3,size(M))*vec(M) ≈ vec(K3*M) ≈ Array(K3, size(M))*vec(M) ≈ sparse(K3,size(M))*vec(M) ≈ BandedMatrix(K3, size(M))*vec(M) ≈ BandedBlockBandedMatrix(K3,size(M))*vec(M)

end

@testset "Biased Upwind : Negative Wind with offside, irregular grid" begin

      M = rand(7)

      L1 = CenteredDifference(1,2,[0.08, 0.02, 0.05, 0.04, 0.07, 0.03],5,-1)
      L2 = CenteredDifference(1,4,[0.08, 0.02, 0.05, 0.04, 0.07, 0.03],5,-1)
      L3 = CenteredDifference(2,3,[0.08, 0.02, 0.05, 0.04, 0.07, 0.03],5,-1)

      K1 = UpwindDifference(1,2,[0.08, 0.02, 0.05, 0.04, 0.07, 0.03],5,-1,offside=1)
      K2 = UpwindDifference(1,4,[0.08, 0.02, 0.05, 0.04, 0.07, 0.03],5,-1,offside=2)
      K3 = UpwindDifference(2,3,[0.08, 0.02, 0.05, 0.04, 0.07, 0.03],5,-1,offside=2)

      @test L1*M ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*M ≈ BandedMatrix(L1, size(M))*M ≈ BandedBlockBandedMatrix(L1,size(M))*M ≈ K1*M ≈ Array(K1, size(M))*vec(M) ≈ sparse(K1,size(M))*M ≈ BandedMatrix(K1, size(M))*M ≈ BandedBlockBandedMatrix(K1,size(M))*M
      @test L2*M ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*M ≈ BandedMatrix(L2, size(M))*M ≈ BandedBlockBandedMatrix(L2,size(M))*M ≈ K2*M ≈ Array(K2, size(M))*vec(M) ≈ sparse(K2,size(M))*M ≈ BandedMatrix(K2, size(M))*M ≈ BandedBlockBandedMatrix(K2,size(M))*M
      @test L3*M ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*M ≈ BandedMatrix(L3, size(M))*M ≈ BandedBlockBandedMatrix(L3,size(M))*M ≈ K3*M ≈ Array(K3, size(M))*vec(M) ≈ sparse(K3,size(M))*M ≈ BandedMatrix(K3, size(M))*M ≈ BandedBlockBandedMatrix(K3,size(M))*M

      M = rand(7,2,2,2)

      @test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L1,size(M))*vec(M) ≈ vec(K1*M) ≈ Array(K1, size(M))*vec(M) ≈ sparse(K1,size(M))*vec(M) ≈ BandedMatrix(K1, size(M))*vec(M) ≈ BandedBlockBandedMatrix(K1,size(M))*vec(M)
      @test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L2,size(M))*vec(M) ≈ vec(K2*M) ≈ Array(K2, size(M))*vec(M) ≈ sparse(K2,size(M))*vec(M) ≈ BandedMatrix(K2, size(M))*vec(M) ≈ BandedBlockBandedMatrix(K2,size(M))*vec(M)
      @test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L3,size(M))*vec(M) ≈ vec(K3*M) ≈ Array(K3, size(M))*vec(M) ≈ sparse(K3,size(M))*vec(M) ≈ BandedMatrix(K3, size(M))*vec(M) ≈ BandedBlockBandedMatrix(K3,size(M))*vec(M)

end

@testset "BC concretizations" begin

      M2 = rand(11)
      M = M2[2:10]
      L1 = CenteredDifference(2,2,1.0,9)
      L2 = CenteredDifference(1,2,1.0,9)

      deriv_start, deriv_end = (L2*M2)[1], (L2*M2)[end]
      params = [1.0,0.5]
      left_RBC = params[1]*M[1] - params[2]*deriv_start
      right_RBC = params[1]*M[end] + params[2]*deriv_end
      bc1 = Dirichlet0BC(Float64)
      bc2 = DirichletBC(0.0,0.0)
      bc3 = NeumannBC((deriv_start,deriv_end),1.0,1)
      bc4 = RobinBC((params[1],-params[2],left_RBC), (params[1],params[2],right_RBC),1.0,1)
      bc5 = GeneralBC([-left_RBC,params[1],-params[2]],[-right_RBC,params[1],params[2]],1.0,1)
      bc6 = PeriodicBC(Float64)

      @test L1*bc1*M ≈ L1*(Array(bc1, length(M))[1]*vec(M) + Array(bc1, length(M))[2]) ≈ L1*(sparse(bc1, length(M))[1]*vec(M) + sparse(bc1, length(M))[2]) ≈ L1*(BandedMatrix(bc1, length(M))[1]*vec(M) + BandedMatrix(bc1, length(M))[2])
      @test L1*bc2*M ≈ L1*(Array(bc2, length(M))[1]*vec(M) + Array(bc2, length(M))[2]) ≈ L1*(sparse(bc2, length(M))[1]*vec(M) + sparse(bc2, length(M))[2]) ≈ L1*(BandedMatrix(bc2, length(M))[1]*vec(M) + BandedMatrix(bc2, length(M))[2])
      @test L1*bc3*M ≈ L1*(Array(bc3, length(M))[1]*vec(M) + Array(bc3, length(M))[2]) ≈ L1*(sparse(bc3, length(M))[1]*vec(M) + sparse(bc3, length(M))[2]) ≈ L1*(BandedMatrix(bc3, length(M))[1]*vec(M) + BandedMatrix(bc3, length(M))[2])
      @test L1*bc4*M ≈ L1*(Array(bc4, length(M))[1]*vec(M) + Array(bc4, length(M))[2]) ≈ L1*(sparse(bc4, length(M))[1]*vec(M) + sparse(bc4, length(M))[2]) ≈ L1*(BandedMatrix(bc4, length(M))[1]*vec(M) + BandedMatrix(bc4, length(M))[2])
      @test L1*bc5*M ≈ L1*(Array(bc5, length(M))[1]*vec(M) + Array(bc5, length(M))[2]) ≈ L1*(sparse(bc5, length(M))[1]*vec(M) + sparse(bc5, length(M))[2]) ≈ L1*(BandedMatrix(bc5, length(M))[1]*vec(M) + BandedMatrix(bc5, length(M))[2])
      @test L1*bc6*M ≈ L1*(Array(bc6, length(M))[1]*vec(M) + Array(bc6, length(M))[2]) ≈ L1*(sparse(bc6, length(M))[1]*vec(M) + sparse(bc6, length(M))[2]) ≈ L1*(BandedMatrix(bc6, length(M))[1]*vec(M) + BandedMatrix(bc6, length(M))[2])
end


text:
@with_kw struct Param

    #Temporal variables for simulation
    t_sim_init          :: Float64   = 0.0
    t_sim_final         :: Float64   = 0.0
    out_every_t         :: Float64   = 1.0
    deltat              :: Float64   = 0.05

    #Max possible sim
    max_runtime         ::String    = "70:00:00"

   #Spacial variables
   rmin            :: Real          =  1.0
   rmax            :: Float64       =  30.0
   rnodes          :: Integer       =  128

   thnodes         :: Integer       =  128

    #Initial Conditions
    A0              ::Float64       = 1.0
    σ               ::Float64       = 2.0
    r0              ::Float64       = 5.0
    ω               ::Float64       = 5.0
    m               ::Integer       = 0.0

    #Cylinder Variables
    ϵ               ::Real      = 0.1
    Ω               ::Real      = 0.3
    Z0              ::Real      = 0.001

    #Data variables
    folder          :: String    = "data/"
    fname           :: String    = "default_file_name"

end


text:
module WhiskerUtils
# Define types
include("whiskerstructs.jl")
include("conversions.jl")

using .WhiskerStructs
using .Conversions

export HydroWhisker, TopoWhisker, hydro2topo

end # module


text:
#== INCLUDE files ==#
include("aggregate.jl")
include("twoassets.jl")

include("solveHJB.jl")
include("solveKFE.jl")

#== Instance of prices ==#
pr = Prices(0.0, 0.04, 0.03, 0.09, 8.0)

#== Instance of TwoAssetProblem ==#
twoap = TwoAssetsProb2(pr; γ = 2.0, ρ = 0.06, ξ = 0.10, χ₀ = 0.08, χ₁= 3.0, τ = 0.0, T = 0.0)

#== Create and FiniteDifference Structure ==#
fde = FDExp(twoap; z = [.8, 1.3], λ = [-1/3 1/3; 1/3 -1/3],
fixedΔa = true, fixedΔb = false, bn = 100, bmin = -2.0, bmax = 40.0, an =70, amax = 70.0, invΔᴷ = 0.05)

#== Solve the HJB equation ==#
solve_hjb!(twoap, fde, fde.sol)

#== Solve the KF equation ==#
solve_fp!(fde, fde.sol, maxit = 4000, tolFK = 1e-9)

include("main_fig.jl");

#== Policies ==#
# fig_pol(fde)
#== Distributions ==#
fig_dens(fde)


text:
using Distributions
using LinearAlgebra
using NeuralArithmetic
using Flux
using Random

Random.seed!(0)

# define the learning task R²->R⁴
function task(x::Vector)
    a, b = x[1], x[2]
    [a+b, a*b, a/b, sqrt(b)]
end
task(x::Matrix) = mapslices(task, x, dims=1)

function generate_data(batchsize::Int)
    a = Uniform(-2,2)
    b = Uniform(0.01,2) # prevent zero/negative inputs to div/sqrt
    p = Product([a,b])
    x = Float32.(rand(p,batchsize))
    y = task(x)
    (x,y)
end

isize = 2     # input size
hsize = 6     # hidden size
osize = 4     # output size
β     = 1e-4  # L₁ regularization

batchsize = 100
nrsteps   = 50000
lr        = 1e-3

model = Chain(NPU(isize,hsize), NAU(hsize,osize))
ps    = Flux.params(model)
data  = [generate_data(batchsize) for _ in 1:nrsteps]
opt   = ADAM(lr)
mse(x,y) = Flux.mse(model(x),y)
loss(x,y) = mse(x,y) + β*norm(ps,1)

cb = Flux.throttle(()->(@info "training..." loss(data[1]...) mse(data[1]...)),1)
Flux.train!(loss, ps, data, opt, cb=cb)

@info "NPU" model[1].Re model[1].Im model[1].g
@info "NAU" model[2].W

using Plots
pyplot()
p1 = heatmap(model[1].Re, yflip=true, aspectratio=1,
             ylabel="NPU Re", clim=(-1,1), c=:bluesreds, colorbar=false)
p2 = heatmap(model[1].Im, yflip=true, aspectratio=1, title="Learned Solution",
             ylabel="NPU Im", clim=(-1,1), c=:bluesreds, colorbar=false)
p4 = heatmap(model[2].W, yflip=true, aspectratio=1,
             ylabel="NAU W", clim=(-1,1), c=:bluesreds)

Re = [1 -1; 1 1; 0 0.5; 1 0; 0 0; 0 1]
Im = zeros(6,2)
g  = ones(2)
W  = [0 0 0 1 0 1;
      0 1 0 0 0 0;
      1 0 0 0 0 0;
      0 0 1 0 0 0]
p5 = heatmap(Re, yflip=true, aspectratio=1,
             ylabel="NPU Re", clim=(-1,1), c=:bluesreds, colorbar=false)
p6 = heatmap(Im, yflip=true, aspectratio=1, title="Perfect Solution",
             ylabel="NPU Im", clim=(-1,1), c=:bluesreds, colorbar=false)
p8 = heatmap(W, yflip=true, aspectratio=1,
             ylabel="NAU W", clim=(-1,1), c=:bluesreds)
plt1 = plot(p1,p2,p4,p5,p6,p8,layout=(2,3), size=(800,370))


text:
using BSplines
struct QuadraticPP{T<:Number,U<:Real}
    a::Vector{T}
    b::Vector{T}
    c::Vector{T}
    x::Vector{U}
    #logx::Vector{T}
    QuadraticPP(T, U, n::Int) = new{T,U}(zeros(T, n), zeros(T, n), zeros(T, n - 1), zeros(U, n))
    QuadraticPP(a::Vector{T}, b::Vector{T}, c::Vector{T}, x::Vector{U}) where {T<:Real,U<:Real} = new{T,U}(a, b, c, x)
end

(p::QuadraticPP{T})(x::T) where {T} = evaluate(p, x)
Base.length(p::QuadraticPP) = length(pp.a)

function evaluate(self::QuadraticPP, z::T) where {T}
    if z <= self.x[1]
        return self.b[1] * (z - self.x[1]) + self.a[1]
    elseif z >= self.x[end]
        rightSlope = self.b[end]
        return rightSlope * (z - self.x[end]) + self.a[end]
    end
    i = searchsortedfirst(self.x, z)  # x[i-1]<z<=x[i]
    if z == self.x[i]
        return self.a[i]
    end
    if i > 1
        i -= 1
    end
    h = z - self.x[i]
    return self.a[i] + h * (self.b[i] + h * (self.c[i]))
end

function evaluateDerivative(self::QuadraticPP, z::T) where {T}
    if z <= self.x[1]
        return self.b[1]
    elseif z >= self.x[end]
        rightSlope = self.b[end]
        return rightSlope
    end
    i = searchsortedfirst(self.x, z)  # x[i-1]<z<=x[i]
    if z == self.x[i]
        return self.a[i]
    end
    if i > 1
        i -= 1
    end
    h = z - self.x[i]
    return self.b[i] + 2h * self.c[i]
end
function convert(::Type{QuadraticPP}, spl::BSplines.Spline)
    t = BSplines.knots(spl.basis)
    n = length(spl.basis) - 1
    T = typeof(spl.coeffs[1])
    a = zeros(T, n)
    b = zeros(T, n)
    c = zeros(T, n - 1)
    x = spl.basis.breakpoints
    α = spl.coeffs
    for i = 1:n
        a[i] = (t[i+2] - t[i+1]) / (t[i+3] - t[i+1]) * (α[i+1] - α[i]) + α[i]
        b[i] = 2 * (α[i+1] - α[i]) / (t[i+3] - t[i+1])
    end
    for i = 1:n-1
        c[i] = ((α[i+2] - α[i+1]) / (t[i+4] - t[i+2]) + (α[i] - α[i+1]) / (t[i+3] - t[i+1])) / (t[i+3] - t[i+2])
    end
    return QuadraticPP(a, b, c, x)
end

function convert(::Type{BSplines.Spline}, pp::QuadraticPP)
    basis = BSplines.BSplineBasis(3, pp.x)
    t = BSplines.knots(basis)
    n = length(pp.x)
    T = typeof(pp.x[1])
    α = zeros(T, n + 1)
    for i = 1:n
        α[i] = pp.a[i] - pp.b[i] / 2 * (t[i+2] - t[i+1])
    end
    α[n+1] = pp.b[n] / 2 * (t[n+3] - t[n+1]) + α[n]
    return BSplines.Spline(basis, α)
end


text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

using Plots;

"""
    demorevolve()

Demonstration of surface construction by revolving a profile curve.

# Examples:
```julia
julia> demorevolve()
```
"""
function demorevolve()
# Construct a test profile in the x-z plane
pnts = [3.0 5.5 5.5 1.5 1.5 4.0 4.5;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.5 1.5 4.5 3.0 7.5 6.0 8.5];
crv = nrbmak(pnts,vec([0 0 0 1/4 1/2 3/4 3/4 1 1 1]));

# rotate and vectrans by some arbitrary amounts.
xx = vecrotz(deg2rad(25)) * vecroty(deg2rad(15)) * vecrotx(deg2rad(20));
nrb = nrbtform(crv,vectrans([5.0;5.0])*xx);

# define axes of rotation
pnt = [5.0;5.0;0.0];
ver = xx * [0.0;0.0;1.0;1.0];
srf = nrbrevolve(nrb,ver[1:3],pnt);

# make and draw nurbs curve
p = nrbeval(srf,collect.([range(0.0,1.0,length=20),range(0.0,1.0,length=20)]));
Plots.pyplot();
plot(p[1,:,:], p[2,:,:], p[3,:,:],
     st = :surface,
     cbar = false,
     legend = false,
     camera = (-30, 30),
     title = "Construct of a 3D surface by revolution of a curve.",
     c=:copper);
# TODO: come back later as pyplot backend does not has light/shading yet
#       furthermore grid lines should be off and interpolated shading should be
#       on
end # demorevolve


text:
# test for the optimization of branch lengths
# with Case g
# Claudia January 2015

## include("../examples/case_g_example.jl");
## q1 = Quartet(1,["6","7","4","8"],[0.5,0.4,0.1]);
## q2 = Quartet(2,["6","7","10","8"],[0.5,0.4,0.1]);
## q3 = Quartet(3,["10","7","4","8"],[0.5,0.4,0.1]);
## q4 = Quartet(4,["6","10","4","8"],[0.5,0.4,0.1]);
## q5 = Quartet(5,["6","7","4","10"],[0.5,0.4,0.1]);

## d = DataCF([q1,q2,q3,q4,q5]);
## extractQuartet!(net,d)

## df = writeExpCF(d.quartet)
## writetable("CaseG_output.csv",df)

include("../src/types.jl")
include("../src/functions.jl")

df = readtable("CaseG_output.csv")
d2 = readTableCF(df)

tree = "((((6,4)1,(7)11#H1:::0.8)5,(11#H1:::0.2,8)),10);" # Case G different starting branch lengths
f = open("prueba_tree.txt","w")
write(f,tree)
close(f)
net = readTopologyUpdate("prueba_tree.txt");
printEdges(net)

net.ht
realht = [0.1,0.2,0.1,1.0]

@time optBL!(net,d2)
# with usual logPseudoLik:
#got 5.33118 at [0.247,0.18076,0.12026,0.2956] after 81 iterations (returned FTOL_REACHED)
#elapsed time: 8.654081812 seconds (8624832 bytes allocated)
#(5.331178696104555,[0.246999,0.180762,0.120265,0.295604])

# with new logPseudoLik:
#got 0.0 at [0.14088,0.19534,0.10493,0.59544] after 512 iterations (returned FTOL_REACHED)
#elapsed time: 6.552122524 seconds (95796384 bytes allocated, 0.56% gc time)
#(5.522931422394296e-8,[0.140883,0.195339,0.10493,0.595437])

# with 100*logPseudoLik and t in (0,10): takes longer, but finds it!
#got 0.0 at [0.10085,0.19991,0.1001,0.98557] after 1123 iterations (returned FTOL_REACHED)
#elapsed time: 53.20413795 seconds (107820188 bytes allocated, 0.13% gc time)
#(2.169777681982341e-9,[0.100853,0.199907,0.100098,0.985569])

# with new ftol, xtol values
#got 0.00073 at [0.2474,0.18189,0.11927,0.30014] after 48 iterations (returned XTOL_REACHED)
#elapsed time: 0.023274544 seconds (4256960 bytes allocated)
#

# with old ftol,xtol values
#got 0.0 at [0.13768,0.19568,0.10456,0.61459] after 561 iterations (returned FTOL_REACHED)
#elapsed time: 2.122653096 seconds (118537484 bytes allocated, 4.50% gc time)

@time optBL!(net,d2,false,1e-5,1e-6,1e-3,1e-4)
#got 0.0001 at [0.24708,0.18075,0.12033,0.29593] after 64 iterations (returned SUCCESS)
#elapsed time: 0.028643709 seconds (5502208 bytes allocated)

@allocated fmin,xmin=optBL!(net,d2)

# -------- different starting point ------

# (0.12,1.,1.,1.)
tree = "((((6,4)1,(7)11#H1:::0.88)5,(11#H1:::0.12,8)),10);" # Case G different starting branch lengths
f = open("prueba_tree.txt","w")
write(f,tree)
close(f)
net = readTopologyUpdate("prueba_tree.txt");
printEdges(net)

@time fmin,xmin=optBL!(net,d2)

# (0.1,1.,1.,1.)
tree = "((((6,4)1,(7)11#H1)5,(11#H1,8)),10);" # Case G different starting branch lengths
f = open("prueba_tree.txt","w")
write(f,tree)
close(f)
net = readTopologyUpdate("prueba_tree.txt");
printEdges(net)

@time optBL!(net,d2)
#got 0.01558 at [0.09049,0.20435,0.09535,0.91173] after 31 iterations (returned XTOL_REACHED)
#elapsed time: 0.014890177 seconds (2772976 bytes allocated)

# (0.05,1.,1.,1.)
tree = "((((6,4)1,(7)11#H1:::0.95)5,(11#H1:::0.05,8)),10);" # Case G different starting branch lengths
f = open("prueba_tree.txt","w")
write(f,tree)
close(f)
net = readTopologyUpdate("prueba_tree.txt");
printEdges(net)

@time fmin,xmin=optBL!(net,d2)

# (0.15,1.,1.,1.)
tree = "((((6,4)1,(7)11#H1:::0.85)5,(11#H1:::0.15,8)),10);" # Case G different starting branch lengths
f = open("prueba_tree.txt","w")
write(f,tree)
close(f)
net = readTopologyUpdate("prueba_tree.txt");
printEdges(net)

@time fmin,xmin=optBL!(net,d2)

# (0.35,1.,1.,1.)
tree = "((((6,4)1,(7)11#H1:::0.65)5,(11#H1:::0.35,8)),10);" # Case G different starting branch lengths
f = open("prueba_tree.txt","w")
write(f,tree)
close(f)
net = readTopologyUpdate("prueba_tree.txt");
printEdges(net)

@time fmin,xmin=optBL!(net,d2)

# (0.11,0.22,0.11,1.1)
tree = "((((6,4)1,(7)11#H1:::0.89)5,(11#H1:::0.11,8)),10);" # Case G different starting branch lengths
f = open("prueba_tree.txt","w")
write(f,tree)
close(f)
net = readTopologyUpdate("prueba_tree.txt");
printEdges(net)
printNodes(net)

@time fmin,xmin=optBL!(net,d2)

# ==================================================================================================

# test optBL with Case H
# Claudia January 2015

## include("../examples/case_h_example.jl");
## q1 = Quartet(1,["6","7","4","8"],[0.5,0.4,0.1]);
## q2 = Quartet(2,["6","7","10","8"],[0.5,0.4,0.1]);
## q3 = Quartet(3,["10","7","4","8"],[0.5,0.4,0.1]);
## q4 = Quartet(4,["6","10","4","8"],[0.5,0.4,0.1]);
## q5 = Quartet(5,["6","7","4","10"],[0.5,0.4,0.1]);

## d = DataCF([q1,q2,q3,q4,q5]);
## extractQuartet!(net,d)

## df = writeExpCF(d.quartet)
## writetable("CaseH_output.csv",df)

include("../src/types.jl")
include("../src/functions.jl")

df = readtable("CaseH_output.csv")
df = readtable("CaseH_output2.csv")
d2 = readTableCF(df)

# starting ht (gamma,t3,t5,t7)
ht = [0.2,1.,1.,1.]
ht = [0.05,0.0001,2.,1.] #strange case un debug18bad.txt

tree = string("((((6:0.1,4:1.5):",string(ht[2]),",#H1:::",string(ht[1]),"):",string(ht[3]),",7:0.2):",string(ht[4]),",(8)#H1:::",string(1-ht[1]),",10:0.1);") # Case H
f = open("prueba_tree.txt","w")
write(f,tree)
close(f)
net = readTopologyUpdate("prueba_tree.txt");
printEdges(net)

net.ht
realht = [0.1,0.1,1.,0.1]

@time optBL!(net,d2)
# with usual logPseudoLik:
#got 4.34098 at [0.1,0.10001,0.99999,0.1] after 151 iterations (returned FTOL_REACHED)
#elapsed time: 21.076511672 seconds (21355496 bytes allocated, 0.17% gc time)
#(4.340977085797292,[0.0999962,0.100005,0.999989,0.0999984])

# with new logPseudoLik:
#got 0.0 at [0.1,0.1,1.0,0.1] after 192 iterations (returned FTOL_REACHED)
#elapsed time: 0.973165584 seconds (16698980 bytes allocated, 3.44% gc time)
#(4.3376522534621724e-13,[0.1,0.0999979,1.0,0.100001])

# with 100*logPseudoLik and t in (0,10):
#got 0.0 at [0.1,0.1,1.0,0.1] after 208 iterations (returned FTOL_REACHED)
#elapsed time: 6.892366587 seconds (20646048 bytes allocated)
#(2.4010737722561867e-13,[0.0999999,0.0999992,1.0,0.1])

# with new ftol, xtol (start 0.2,1,1,1)
#got 0.04565 at [0.20363,0.64827,0.40823,0.11374] after 51 iterations (returned XTOL_REACHED)
#elapsed time: 0.057075544 seconds (4462816 bytes allocated, 53.14% gc time)

#got 0.04478 at [0.20178,0.64205,0.41684,0.11366] after 25 iterations (returned SUCCESS)
#elapsed time: 0.016354695 seconds (2226424 bytes allocated)

# with old ftol, xol
#got 0.0 at [0.1,0.1,1.0,0.1] after 176 iterations (returned FTOL_REACHED)
#elapsed time: 0.117537385 seconds (15143800 bytes allocated, 29.95% gc time)

@time optBL!(net,d2,false,1e-5,1e-6,1e-5,1e-6)
#got 1.0e-5 at [0.10043,0.1054,0.99408,0.10017] after 115 iterations (returned FTOL_REACHED)
#elapsed time: 0.055701792 seconds (10398360 bytes allocated)

@time optBL!(net,d2,false,1e-5,1e-6,1e-3,1e-4)
#got 1.0e-5 at [0.10054,0.10558,0.99394,0.10022] after 25 iterations (returned FTOL_REACHED)
#elapsed time: 0.013228257 seconds (2226504 bytes allocated)

# from debug18bad
tree = string("(4,6,(#H2:0.7392085405544356::0.046179825120885414,(7,(10,(8)#H2:0.0::0.9538201748791146):0.9803511144374873):2.212878358589699):0.00000038687);")
f = open("prueba_tree.txt","w")
write(f,tree)
close(f)
net = readTopologyUpdate("prueba_tree.txt");
printEdges(net)

net.ht
realht = [0.1,0.1,1.,0.1]

@time optBL!(net,d2)

# ==================================================================================================================================

# test optBL with Case J
# Claudia January 2015

## include("../examples/case_j_example.jl");
## q1 = Quartet(1,["6","7","4","8"],[0.5,0.4,0.1]);
## q2 = Quartet(2,["6","7","10","8"],[0.5,0.4,0.1]);
## q3 = Quartet(3,["10","7","4","8"],[0.5,0.4,0.1]);
## q4 = Quartet(4,["6","10","4","8"],[0.5,0.4,0.1]);
## q5 = Quartet(5,["6","7","4","10"],[0.5,0.4,0.1]);

## d = DataCF([q1,q2,q3,q4,q5]);
## extractQuartet!(net,d)

## df = writeExpCF(d.quartet)
## writetable("CaseJ_output.csv",df)

include("../src/types.jl")
include("../src/functions.jl")

df = readtable("CaseJ_output.csv")
d2 = readTableCF(df)

# starting ht (gamma,t3,t5,t7)
ht = [0.2,1.,1.,1.]

tree = string("((((6)#H1:::",string(1-ht[1]),",4:1.5):",string(ht[2]),",7:0.2):",string(ht[3]),",8:0.1,(#H1:::",string(ht[1]),",10));") # Case J
f = open("prueba_tree.txt","w")
write(f,tree)
close(f)
net = readTopologyUpdate("prueba_tree.txt");
printEdges(net)

net.ht
realht = [0.1,0.2,0.1,1.0]

@time optBL!(net,d2)
# with usual logPseudoLik:
#got 5.39161 at [0.09997,0.19999,0.1,1.00061] after 246 iterations (returned FTOL_REACHED)
#elapsed time: 59.297791699 seconds (71939204 bytes allocated, 0.11% gc time)
#(5.3916134028946034,[0.0999661,0.199991,0.0999991,1.00061])

# with new logPseudoLik:
#got 0.0 at [0.10002,0.20001,0.1,0.99958] after 260 iterations (returned FTOL_REACHED)
#elapsed time: 1.211645168 seconds (21966196 bytes allocated)
#(9.797186804450541e-12,[0.100022,0.200007,0.100001,0.999579])

# with 100*logPseudoLik and t in (0,10):
#got 0.0 at [0.1,0.2,0.1,1.00003] after 249 iterations (returned FTOL_REACHED)
#elapsed time: 6.656956697 seconds (23010204 bytes allocated, 0.56% gc time)
#(5.280560068867562e-12,[0.0999983,0.199999,0.0999999,1.00003])

# with ftol, xtol
#got 0.05209 at [0.26015,0.26161,0.10523,0.24363] after 42 iterations (returned XTOL_REACHED)
#elapsed time: 0.021402455 seconds (3634568 bytes allocated)

# with old ftol, xtol
#got 0.0 at [0.1,0.2,0.1,1.0] after 233 iterations (returned FTOL_REACHED)
#elapsed time: 0.112556852 seconds (19426488 bytes allocated)

@time optBL!(net,d2,false,1e-5,1e-6,1e-3,1e-4)
#got 0.00219 at [0.1303,0.20847,0.10099,0.6283] after 106 iterations (returned FTOL_REACHED)
#elapsed time: 0.060469424 seconds (8887448 bytes allocated)

@time optBL!(net,d2,false,1e-5,1e-6,1e-5,1e-6)
#got 0.00177 at [0.13073,0.20994,0.10121,0.64478] after 26 iterations (returned FTOL_REACHED)
#elapsed time: 0.012715507 seconds (2277912 bytes allocated)

# ==================================================================================================================================

# test optBL with Case F Bad Diamond I
# Claudia January 2015

## include("../examples/case_f_example.jl");
## parameters!(net)
## q1 = Quartet(1,["6","7","4","8"],[0.5,0.4,0.1]);
## q2 = Quartet(2,["6","7","10","8"],[0.5,0.4,0.1]);
## q3 = Quartet(3,["10","7","4","8"],[0.5,0.4,0.1]);
## q4 = Quartet(4,["6","10","4","8"],[0.5,0.4,0.1]);
## q5 = Quartet(5,["6","7","4","10"],[0.5,0.4,0.1]);

## d = DataCF([q1,q2,q3,q4,q5]);
## extractQuartet!(net,d)

## df = writeExpCF(d.quartet)
## writetable("CaseF_output.csv",df)

include("../src/types.jl")
include("../src/functions.jl")

df = readtable("CaseF_output.csv")
d2 = readTableCF(df)

# starting ht (gamma,t4,t5,t9)
ht = [0.1,1.,1.,1.]

tree = string("(((6:0.1,(4)11#H1:::",string(1-ht[1]),")1:",string(ht[3]),",(11#H1:::",string(ht[1]),",7))5:",string(ht[4])",8:0.1,10:0.1);") # Case F: bad diamond I
f = open("prueba_tree.txt","w")
write(f,tree)
close(f)
net = readTopologyUpdate("prueba_tree.txt");
printEdges(net)

net.ht
realht = [0.1,0.127,0.0285]

@time optBL!(net,d2)
# with usual logPseudoLik:
#got 5.40235 at [0.1,0.12689,0.02855] after 116 iterations (returned FTOL_REACHED)
#elapsed time: 15.648447251 seconds (17823396 bytes allocated, 0.18% gc time)
#(5.402353356033268,[0.1,0.126887,0.0285486])

# with new logPseudoLik:
#got 0.0 at [0.1,0.12689,0.02855] after 116 iterations (returned FTOL_REACHED)
#elapsed time: 0.494250509 seconds (9272420 bytes allocated)
#(3.6216219428084243e-12,[0.1,0.126887,0.0285488])

# with new logPseudoLik and new algorithm LN_COBYLA
#got 0.0 at [0.1,0.12689,0.02855] after 3277 iterations (returned FTOL_REACHED)
#elapsed time: 23.001580343 seconds (312597904 bytes allocated, 0.98% gc time)
#(1.920897287842076e-13,[0.0999998,0.126889,0.0285491])

# with 100*logPseudoLik and t in (0,10):
# got 0.0 at [0.1,0.12689,0.02855] after 3356 iterations (returned FTOL_REACHED)
# elapsed time: 90.405727385 seconds (281649264 bytes allocated, 0.24% gc time)
# (8.157692309723136e-12,[0.0999999,0.126889,0.028549])

# with ftol,xtol
#got 0.25041 at [0.06362,0.15339,0.04918] after 24 iterations (returned XTOL_REACHED)
#elapsed time: 0.045551191 seconds (2403576 bytes allocated)

#with old ftol, xtol
#got 0.0 at [0.1,0.12689,0.02855] after 2917 iterations (returned FTOL_REACHED)
#elapsed time: 1.303746554 seconds (229545768 bytes allocated, 11.98% gc time)

@time optBL!(net,d2,false,1e-5,1e-6,1e-3,1e-4)
#got 5.0e-5 at [0.09974,0.12735,0.02892] after 92 iterations (returned XTOL_REACHED)
#elapsed time: 0.079739994 seconds (7321368 bytes allocated, 45.28% gc time)

# ==================================================================================================================================

# test optBL with Case I Bad Diamond II
# Claudia January 2015

## include("../examples/case_i_example.jl");
## q1 = Quartet(1,["6","7","4","8"],[0.5,0.4,0.1]);
## q2 = Quartet(2,["6","7","10","8"],[0.5,0.4,0.1]);
## q3 = Quartet(3,["10","7","4","8"],[0.5,0.4,0.1]);
## q4 = Quartet(4,["6","10","4","8"],[0.5,0.4,0.1]);
## q5 = Quartet(5,["6","7","4","10"],[0.5,0.4,0.1]);

## d = DataCF([q1,q2,q3,q4,q5]);
## extractQuartet!(net,d)

## df = writeExpCF(d.quartet)
## writetable("CaseI_output.csv",df)

include("../src/types.jl")
include("../src/functions.jl")

df = readtable("CaseI_output.csv")
d2 = readTableCF(df)

# starting ht (gamma,t4,t6,t9,t10)
ht = [0.2,0.0,2.0,2.0,2.0]
ht = [0.2,0.0,0.5,0.5,0.5]

tree = string("((((8,10):",string(ht[2]),")#H1:::",string(1-ht[1]),",7):",string(ht[3]),",6,(4,#H1:",string(ht[4]),"::",string(ht[1]),"):",string(ht[5]),");") # Case I Bad diamond II
f = open("prueba_tree.txt","w")
write(f,tree)
close(f)
net = readTopologyUpdate("prueba_tree.txt");
printEdges(net)

net.ht
realht = [0.1,2.0,1.0,1.0,1.0]

@time optBL!(net,d2)
## with usual logPseudoLik:
# with ht = [0.2,0.0,0.5,0.5,0.5]
#got 3.69564 at [0.09753,1.97129,0.99532,0.46822,1.04214] after 267 iterations (returned FTOL_REACHED)
#elapsed time: 45.70869142 seconds (35023452 bytes allocated, 0.06% gc time)
#(3.6956398763931144,[0.0975269,1.97129,0.995318,0.468224,1.04214])

#ht: [0.2,1.0,2.0,2.0,2.0]
#got 3.69564 at [0.1034,2.03878,1.00653,3.58586,0.94929] after 176 iterations (returned FTOL_REACHED)
#elapsed time: 15.743171805 seconds (17639780 bytes allocated)
#(3.6956415275081724,[0.103403,2.03878,1.00653,3.58586,0.949288])

## with new logPseudoLik:
#ht: [0.2,1.0,2.0,2.0,2.0]
#got 0.0 at [0.10342,2.03906,1.00656,3.60666,0.94899] after 195 iterations (returned FTOL_REACHED)
#elapsed time: 1.11298624 seconds (17529096 bytes allocated, 2.98% gc time)
#(3.3635155002465492e-6,[0.103419,2.03906,1.00656,3.60666,0.948988])

# with ht = [0.2,0.0,0.5,0.5,0.5]
#got 0.0 at [0.09997,1.99969,0.99995,0.99284,1.00042] after 424 iterations (returned FTOL_REACHED)
#elapsed time: 1.733474741 seconds (36994324 bytes allocated, 1.90% gc time)
#(1.8811458391292683e-10,[0.0999737,1.99969,0.99995,0.992835,1.00042])

# with 100*logPseudoLik and t in (0,10):
# ht = [0.2,0.0,2.0,2.0,2.0]
#got 0.0 at [0.1,2.0,1.0,0.99994,1.0] after 753 iterations (returned FTOL_REACHED)
#elapsed time: 47.213517434 seconds (80265788 bytes allocated, 0.16% gc time)
#(1.2568189823957448e-12,[0.0999998,2.0,1.0,0.999941,1.0])

# with ftol,xtol
# ht = [0.2,0.0,2.0,2.0,2.0]
#got 0.00256 at [0.10406,2.0465,1.00637,3.33775,0.98475] after 58 iterations (returned XTOL_REACHED)
#elapsed time: 0.034984668 seconds (5205128 bytes allocated)

# with old ftol, xtol
#got 0.00034 at [0.1034,2.03863,1.00652,3.5823,0.94937] after 201 iterations (returned FTOL_REACHED)
#elapsed time: 0.094041162 seconds (17632400 bytes allocated)

@time optBL!(net,d2,false,1e-5,1e-6,1e-3,1e-4)
# ht = [0.2,0.0,2.0,2.0,2.0]
#got 0.00043 at [0.10287,2.02652,1.00596,5.93267,0.96293] after 51 iterations (returned SUCCESS)
#elapsed time: 0.027747237 seconds (4558192 bytes allocated)

@time optBL!(net,d2,false,1e-5,1e-6,1e-5,1e-6)
#got 0.00036 at [0.10311,2.03708,1.00602,3.43898,0.95781] after 35 iterations (returned FTOL_REACHED)
#elapsed time: 0.018274282 seconds (3167744 bytes allocated)

# -------------------5taxon tree------------------
## include("../tree_example_read.jl");
## printEdges(net)
## q1 = Quartet(1,["6","7","4","8"],[0.5,0.4,0.1]);
## q2 = Quartet(2,["6","7","10","8"],[0.5,0.4,0.1]);
## q3 = Quartet(3,["10","7","4","8"],[0.5,0.4,0.1]);
## q4 = Quartet(4,["6","10","4","8"],[0.5,0.4,0.1]);
## q5 = Quartet(5,["6","7","4","10"],[0.5,0.4,0.1]);

## d = DataCF([q1,q2,q3,q4,q5]);
## extractQuartet!(net,d)

## df = writeExpCF(d.quartet)
## writetable("Tree_output.csv",df)

include("../src/types.jl")
include("../src/functions.jl")

df = readtable("Tree_output.csv")
d2 = readTableCF(df)

# starting tree:
ht = [1.0,1.0]
tree = string("(((6:0.1,4:1.5)1:",string(ht[1]),",7:0.2)5:",string(ht[2]),",8:0.1,10:0.1);") # normal tree
f = open("prueba_tree.txt","w")
write(f,tree)
close(f)
net = readTopologyUpdate("prueba_tree.txt");

net.ht
realht = [0.2,0.1]

@time optBL!(net,d2)
#got 5.34957 at [0.2,0.1] after 28 iterations (returned FTOL_REACHED)
#elapsed time: 5.533522804 seconds (54162200 bytes allocated, 0.41% gc time)
#(5.349567420518451,[0.2,0.1])

#with ftol, xtol
#got 0.0 at [0.19999,0.09999] after 20 iterations (returned XTOL_REACHED)
#elapsed time: 0.007440742 seconds (1234840 bytes allocated)

#==============================================================================================
#================ Debugging optBL ============================================================
# compare expCF from wrong estimates with real estimates
# ===========================================================================================

# test optBL with Case I Bad Diamond II
# does not yield correct ht for one starting point

include("../src/types.jl")
include("../src/functions.jl")

df = readtable("CaseI_output.csv")
d2 = readTableCF(df)

# starting ht (gamma,t4,t6,t9,t10)
wronght = [0.1,1.0,1.0,3.6,1.0]
ht = wronght
realht = [0.1,2.0,1.0,1.0,1.0]

tree = string("((((8,10):",string(ht[2]),")#H1:::",string(1-ht[1]),",7):",string(ht[3]),",6,(4,#H1:",string(ht[4]),"::",string(ht[1]),"):",string(ht[5]),");") # Case I Bad diamond II
f = open("prueba_tree.txt","w")
write(f,tree)
close(f)
net = readTopologyUpdate("prueba_tree.txt");
printEdges(net)

net.ht

q1 = Quartet(1,["6","7","4","8"],[0.5,0.4,0.1]);
q2 = Quartet(2,["6","7","10","8"],[0.5,0.4,0.1]);
q3 = Quartet(3,["10","7","4","8"],[0.5,0.4,0.1]);
q4 = Quartet(4,["6","10","4","8"],[0.5,0.4,0.1]);
q5 = Quartet(5,["6","7","4","10"],[0.5,0.4,0.1]);

d = DataCF([q1,q2,q3,q4,q5]);
extractQuartet!(net,d)

wrongdf = writeExpCF(d.quartet)
writetable("CaseI_output_wrong.csv",wrongdf)

# test optBL with Case G
# does not yield correct ht for one starting point
include("../src/types.jl")
include("../src/functions.jl")

df = readtable("CaseG_output.csv")
d2 = readTableCF(df)

# starting ht (gamma,t3,t6,t9)
wronght = [0.14,0.2,0.1,0.6]
ht = wronght
realht = [0.1,0.2,0.1,1.0]

tree = string("((((6,4)1:",string(ht[2]),",(7)11#H1:::",string(1-ht[1]),")5:",string(ht[3]),",(11#H1:::",string(ht[1]),",8):",string(ht[4]),"),10);") # Case G different starting branch lengths
f = open("prueba_tree.txt","w")
write(f,tree)
close(f)
net = readTopologyUpdate("prueba_tree.txt");
printEdges(net

net.ht

q1 = Quartet(1,["6","7","4","8"],[0.5,0.4,0.1]);
q2 = Quartet(2,["6","7","10","8"],[0.5,0.4,0.1]);
q3 = Quartet(3,["10","7","4","8"],[0.5,0.4,0.1]);
q4 = Quartet(4,["6","10","4","8"],[0.5,0.4,0.1]);
q5 = Quartet(5,["6","7","4","10"],[0.5,0.4,0.1]);

d = DataCF([q1,q2,q3,q4,q5]);
extractQuartet!(net,d)

wrongdf = writeExpCF(d.quartet)
writetable("CaseG _output_wrong.csv",wrongdf)