text:
using KNITRO
using Test

const MPS_PROBLEM =
"""
NAME         lo1
OBJSENSE     MAX
ROWS
 N  obj
 E  c1
 G  c2
 L  c3
COLUMNS
    x1        obj       3
    x1        c1        3
    x1        c2        2
    x2        obj       1
    x2        c1        1
    x2        c2        1
    x2        c3        2
    x3        obj       5
    x3        c1        2
    x3        c2        3
    x4        obj       1
    x4        c2        1
    x4        c3        3
RHS
    rhs       c1        30
    rhs       c2        15
    rhs       c3        25
RANGES
BOUNDS
 UP bound     x2        10
ENDATA
"""

@testset "Instantiation Knitro C interface" begin
    # get KNITRO.KNITRO release version
    rel = KNITRO.get_release()
    @test isa(rel, String)

    @testset "Definition of model" begin
        m = KNITRO.Model()
        options = joinpath(dirname(@__FILE__), "..", "examples", "knitro.opt")
        KNITRO.KN_reset_params_to_defaults(m)

        KNITRO.KN_free(m)
        @test m.env.ptr_env == C_NULL
    end
end

# add generic callbacks for future tests

function evalAll(kc, cb, evalRequest, evalResult, userParams)
    x = evalRequest.x
    evalRequestCode = evalRequest.evalRequestCode

    if evalRequestCode == KNITRO.KN_RC_EVALFC
        # Evaluate nonlinear objective
        evalResult.obj[1] = x[1]^2 * x[3] + x[2]^3 * x[3]^2
    elseif evalRequestCode == KNITRO.KN_RC_EVALGA
        evalResult.objGrad[1] = 2 * x[1] * x[3]
        evalResult.objGrad[2] = 3 * x[2]^2 * x[3]^2
        evalResult.objGrad[3] = x[1]^2 + 2 * x[2]^3 * x[3]
    elseif evalRequestCode == KNITRO.KN_RC_EVALH
        evalResult.hess[1] = 2 * x[3]
        evalResult.hess[2] = 2 * x[1]
        evalResult.hess[3] = 6 * x[2] * x[3]^2
        evalResult.hess[4] = 6 * x[2]^2 * x[3]
        evalResult.hess[5] = 2 * x[2]^3
    elseif evalRequestCode == KNITRO.KN_RC_EVALHV
        vec = evalRequest.vec
        evalResult.hessVec[1] = (2 * x[3]) * vec[1] + (2 * x[1]) * vec[3]
        evalResult.hessVec[2] = (6 * x[2] * x[3]^2) * vec[2] + (6 * x[2]^2 * x[3]) * vec[3]
        evalResult.hessVec[3] = (2 * x[1]) * vec[1] + (6 * x[2]^2 * x[3]) * vec[2] + (2 * x[2]^3) * vec[3]

    elseif evalRequestCode == KNITRO.KN_RC_EVALH_NO_F
        evalResult.hess[1] = 0
        evalResult.hess[2] = 0
        evalResult.hess[3] = 0
        evalResult.hess[4] = 0
        evalResult.hess[5] = 0

    elseif evalRequestCode == KNITRO.KN_RC_EVALHV_NO_F
        vec = evalRequest.vec
        evalResult.hessVec[1] = 0
        evalResult.hessVec[2] = lambda_[4] * vec[3]
        evalResult.hessVec[3] = lambda_[4] * vec[2]

    else
        return KNITRO.KN_RC_CALLBACK_ERR
    end
    return 0
end

function callback(name)
    function callbackFn(kc, x, lambda_, userParams)
        println(name)
        return 0
    end
    return callbackFn
end

if KNITRO.KNITRO_VERSION >= v"12.0"
    @testset "Names getters" begin
        kc = KNITRO.KN_new()
        KNITRO.KN_add_vars(kc, 3)
        KNITRO.KN_add_cons(kc, 3)
        xnames = ["x1", "x2", "x3"]
        cnames = ["c1", "c2", "c3"]

        KNITRO.KN_set_var_names(kc, xnames)
        KNITRO.KN_set_con_names(kc, cnames)
        index = Cint(1)

        name = KNITRO.KN_get_var_names(kc, index)
        @test name == xnames[2]
        outnames = KNITRO.KN_get_var_names(kc, [index])
        @test outnames[1] == xnames[2]
        outnames = KNITRO.KN_get_var_names(kc)
        @test xnames == outnames

        name = KNITRO.KN_get_con_names(kc, index)
        @test name == cnames[2]
        outnames = KNITRO.KN_get_con_names(kc, [index])
        @test outnames[1] == cnames[2]
        outnames = KNITRO.KN_get_con_names(kc)
        @test cnames == outnames

        KNITRO.KN_free(kc)
    end
end
if KNITRO.KNITRO_VERSION >= v"12.1"
    @testset "MPS reader/writer" begin
        mps_name = joinpath(dirname(@__FILE__), "lp.mps")
        mps_name_out = joinpath(dirname(@__FILE__), "lp2.mps")
        open(mps_name, "w") do io
            write(io, MPS_PROBLEM)
        end
        kc = KNITRO.KN_new()
        KNITRO.KN_load_mps_file(kc, mps_name)
        KNITRO.KN_set_param(kc, "outlev", 0)
        KNITRO.KN_write_mps_file(kc, mps_name_out)
        status = KNITRO.KN_solve(kc)
        obj = KNITRO.get_objective(kc)
        KNITRO.KN_free(kc)
        @test status == 0
        @test obj ≈ 250.0 / 3.0

        # Resolve with dumped MPS file
        kc = KNITRO.KN_new()
        KNITRO.KN_load_mps_file(kc, mps_name_out)
        KNITRO.KN_set_param(kc, "outlev", 0)
        status = KNITRO.KN_solve(kc)
        obj = KNITRO.get_objective(kc)
        KNITRO.KN_free(kc)
        @test status == 0
        @test obj ≈ 250.0 / 3.0
    end
end

@testset "First problem" begin
    kc = KNITRO.KN_new()
    @test isa(kc, KNITRO.Model)
    # By default, kc does not have any callback
    @test !KNITRO.has_callbacks(kc)

    release = KNITRO.get_release()

    KNITRO.KN_reset_params_to_defaults(kc)

    options = joinpath(dirname(@__FILE__), "..", "examples", "test_knitro.opt")
    tuner1 = joinpath(dirname(@__FILE__), "..", "examples", "tuner-fixed.opt")
    tuner2 = joinpath(dirname(@__FILE__), "..", "examples", "tuner-explore.opt")
    KNITRO.KN_set_param(kc, "algorithm", 0)
    KNITRO.KN_set_param(kc, "cplexlibname", ".")
    KNITRO.KN_set_param(kc, "xtol", 1e-15)
    KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_ALG, KNITRO.KN_ALG_BAR_DIRECT)
    KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_CPLEXLIB, ".")
    KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_XTOL, 1e-15)

    @test KNITRO.KN_get_int_param(kc, "algorithm") == KNITRO.KN_ALG_BAR_DIRECT
    @test KNITRO.KN_get_double_param(kc, "xtol") == 1e-15
    @test KNITRO.KN_get_int_param(kc, KNITRO.KN_PARAM_ALG) == KNITRO.KN_ALG_BAR_DIRECT
    @test KNITRO.KN_get_double_param(kc, KNITRO.KN_PARAM_XTOL) == 1e-15
    @test KNITRO.KN_get_param_name(kc, KNITRO.KN_PARAM_XTOL)  == "xtol"

    @test KNITRO.KN_get_param_doc(kc, KNITRO.KN_PARAM_XTOL) == "# Step size tolerance used for terminating the optimization.\n"
    @test KNITRO.KN_get_param_type(kc, KNITRO.KN_PARAM_XTOL) == KNITRO.KN_PARAMTYPE_FLOAT
    @test KNITRO.KN_get_num_param_values(kc, KNITRO.KN_PARAM_XTOL) == 0

    @test KNITRO.KN_get_param_value_doc(kc, KNITRO.KN_PARAM_GRADOPT, 1) == "exact"

    @test KNITRO.KN_get_param_id(kc, "xtol") == KNITRO.KN_PARAM_XTOL

    # START: Some specific parameter settings
    KNITRO.KN_set_param(kc, "hessopt", 1)
    KNITRO.KN_set_param(kc, "presolve", 0)
    KNITRO.KN_set_param(kc, "outlev", 0)
    # END:   Some specific parameter settings

    # Perform a derivative check.
    KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_DERIVCHECK, KNITRO.KN_DERIVCHECK_ALL)

    function newpt_callback(kc, x, lambda_, user_data)
        a = KNITRO.KN_get_rel_feas_error(kc)
        b = KNITRO.KN_get_rel_opt_error(kc)
        return 0
    end

    # Define objective goal
    objGoal = KNITRO.KN_OBJGOAL_MAXIMIZE
    KNITRO.KN_set_obj_goal(kc, objGoal)

    # Add the variables and set their bounds.
    nV = 3
    KNITRO.KN_add_vars(kc, nV)
    KNITRO.KN_set_var_lobnds(kc, [0, 0.1, 0])
    KNITRO.KN_set_var_upbnds(kc, [0., 2, 2])
    # Define an initial point.
    KNITRO.KN_set_var_primal_init_values(kc, [1., 1, 1.5])
    KNITRO.KN_set_var_dual_init_values(kc,  [1., 1, 1, 1])

    # Add the constraints and set their bounds.
    nC = 1
    KNITRO.KN_add_cons(kc, nC)
    KNITRO.KN_set_con_lobnds(kc, [0.1])
    KNITRO.KN_set_con_upbnds(kc, [2 * 2 * 0.99])

    # Test getters.
    if KNITRO.KNITRO_VERSION >= v"12.0"
        xindex = Cint[0, 1, 2]
        @test KNITRO.KN_get_var_lobnds(kc, xindex) == [0, 0.1, 0]
        @test KNITRO.KN_get_var_upbnds(kc, xindex) == [0., 2, 2]

        cindex = Cint[0]
        @test KNITRO.KN_get_con_lobnds(kc, cindex) == [0.1]
        @test KNITRO.KN_get_con_upbnds(kc, cindex) == [2 * 2 * 0.99]
    end

    # Load quadratic structure x1*x2 for the constraint.
    KNITRO.KN_add_con_quadratic_struct(kc, 0, 1, 2, 1.0)

    # Define callback functions.
    cb = KNITRO.KN_add_objective_callback(kc, evalAll)
    @test KNITRO.has_callbacks(kc)
    KNITRO.KN_set_cb_grad(kc, cb, evalAll)
    KNITRO.KN_set_cb_hess(kc, cb, 5, evalAll,
                          hessIndexVars1=Int32[0, 0, 1, 1, 2],
                          hessIndexVars2=Int32[0, 2, 1, 2, 2])

    KNITRO.KN_set_newpt_callback(kc, newpt_callback)

    # Add complementarity constraints.
    KNITRO.KN_set_compcons(kc, [KNITRO.KN_CCTYPE_VARVAR], Int32[0], Int32[1])

    # Solve the problem.
    status = KNITRO.KN_solve(kc)
    @test status == 0

    # Restart using the previous solution.
    nStatus, objSol, x, lambda_ = KNITRO.KN_get_solution(kc)
    KNITRO.KN_set_var_primal_init_values(kc, x)
    KNITRO.KN_set_var_dual_init_values(kc, lambda_)
    status = KNITRO.KN_solve(kc)
    @test status == 0

    # Restart with new variable bounds
    KNITRO.KN_set_var_lobnds(kc, Float64[0., 0, 0])
    KNITRO.KN_set_var_upbnds(kc, Float64[2., 2, 2])
    status = KNITRO.KN_solve(kc)
    @test status == 0

    # Retrieve relevant solve information
    @test KNITRO.KN_get_number_FC_evals(kc) >= 1
    @test KNITRO.KN_get_number_GA_evals(kc) >= 1
    @test KNITRO.KN_get_number_H_evals(kc) >= 1
    @test KNITRO.KN_get_number_HV_evals(kc) == 0
    @test KNITRO.KN_get_number_iters(kc) >= 1
    @test KNITRO.KN_get_number_cg_iters(kc) >= 0
    @test KNITRO.KN_get_abs_feas_error(kc) < 1e-10
    @test KNITRO.KN_get_rel_feas_error(kc) < 1e-10
    @test KNITRO.KN_get_abs_opt_error(kc) < 1e-7
    @test KNITRO.KN_get_rel_opt_error(kc) < 1e-8
    @test KNITRO.KN_get_con_values(kc)[1] ≈ 3.96

    nStatus, objSol, x, lambda_ = KNITRO.KN_get_solution(kc)
    @test nStatus == 0
    @test x ≈ [0., 2.0, 1.98]

    @test objSol ≈ 31.363199 atol=1e-5

    # Test getters for primal and dual variables
    if KNITRO.KNITRO_VERSION >= v"12.0"
        xopt = KNITRO.KN_get_var_primal_values(kc, Cint[0, 1, 2])
        @test xopt == x
        rc = KNITRO.KN_get_var_dual_values(kc, Cint[0, 1, 2])
        @test rc == lambda_[2:4]
        dual = KNITRO.KN_get_con_dual_values(kc, Cint[0])
        @test dual == [lambda_[1]]
    end

    KNITRO.KN_free(kc)
end

@testset "Second problem test" begin
    kc = KNITRO.KN_new()

    function prettyPrinting(str, userParams)
        s = "KNITRO-Julia: " * str
        println(s)
        return length(s)
    end

    KNITRO.KN_set_puts_callback(kc, prettyPrinting)

    # START: Some specific parameter settings
    KNITRO.KN_set_param(kc, "outlev", 0)
    KNITRO.KN_set_param(kc, "presolve", 0)
    KNITRO.KN_set_param(kc, "ms_enable", 1)
    KNITRO.KN_set_param(kc, "ms_maxsolves", 5)
    KNITRO.KN_set_param(kc, "hessian_no_f", 1)
    KNITRO.KN_set_param(kc, "hessopt", KNITRO.KN_HESSOPT_PRODUCT)
    # END:   Some specific parameter settings

    # Define objective goal
    objGoal = KNITRO.KN_OBJGOAL_MAXIMIZE
    KNITRO.KN_set_obj_goal(kc, objGoal)

    # Add the variables and set their bounds.
    nV = 3
    KNITRO.KN_add_vars(kc, nV)
    KNITRO.KN_set_var_lobnds(kc, [0, 0.1, 0])
    KNITRO.KN_set_var_upbnds(kc, [0., 2, 2])

    # Define an initial point.
    KNITRO.KN_set_var_primal_init_values(kc, [1, 1, 1.5])
    KNITRO.KN_set_var_dual_init_values(kc, [1, 1, 1, 1.])

    # Add the constraints and set their lower bounds.
    nC = 1
    KNITRO.KN_add_cons(kc, nC)
    KNITRO.KN_set_con_lobnds(kc, [0.1])
    KNITRO.KN_set_con_upbnds(kc, [2 * 2 * 0.99])

    # Load quadratic structure x1*x2 for the constraint.
    KNITRO.KN_add_con_quadratic_struct(kc, 0, 1, 2, 1.0)

    # Define callback functions.
    cb = KNITRO.KN_add_objective_callback(kc, evalAll)

    KNITRO.KN_set_cb_grad(kc, cb, evalAll)
    KNITRO.KN_set_cb_hess(kc, cb, KNITRO.KN_DENSE_ROWMAJOR, evalAll)

    KNITRO.KN_set_ms_process_callback(kc, callback("ms callback"))

    function ms_initpt_callbackFn(kc, nSolveNumber, x, lambda_, userParams)
        x[:] = [1, 1, 1.1 + 0.1 * nSolveNumber]
        lambda_[:] = [1., 1, 1, 1]
        return 0
    end

    # Set multistart initial point callback
    KNITRO.KN_set_ms_initpt_callback(kc, ms_initpt_callbackFn)

    # Add complementarity constraints.
    KNITRO.KN_set_compcons(kc, [KNITRO.KN_CCTYPE_VARVAR], Int32[0], Int32[1])

    # Solve the problem.
    status = KNITRO.KN_solve(kc)

    @test status == 0
    nStatus, objSol, x, lambda_ = KNITRO.KN_get_solution(kc)
    @test nStatus == 0
    @test x ≈ [0., 2.0, 1.98] atol=1e-5

    @test objSol ≈ 31.363199 atol=1e-5

    KNITRO.KN_free(kc)
end

@testset "Third problem test" begin
    kc = KNITRO.KN_new()

    KNITRO.KN_set_param(kc, "outlev", 0)
    KNITRO.KN_set_param(kc, "presolve", KNITRO.KN_PRESOLVEDBG_NONE)

    # Define objective goal
    objGoal = KNITRO.KN_OBJGOAL_MAXIMIZE
    KNITRO.KN_set_obj_goal(kc, objGoal)

    # Add the variables and set their bounds.
    nV = 3
    KNITRO.KN_add_vars(kc, nV)
    KNITRO.KN_set_var_lobnds(kc, [0, 0.1, 0])
    KNITRO.KN_set_var_upbnds(kc, [0., 2, 2])

    # Define an initial point.
    KNITRO.KN_set_var_primal_init_values(kc, [1, 1, 1.5])
    KNITRO.KN_set_var_dual_init_values(kc, [1., 1, 1, 1])

    # Add the constraints and set their lower bounds.
    nC = 1
    KNITRO.KN_add_cons(kc, nC)
    KNITRO.KN_set_con_lobnds(kc, [0.1])
    KNITRO.KN_set_con_upbnds(kc, [2 * 2 * 0.99])

    # Load quadratic structure x1*x2 for the constraint.
    KNITRO.KN_add_con_quadratic_struct(kc, 0, 1, 2, 1.0)

    # Define callback functions.
    cb = KNITRO.KN_add_objective_callback(kc, evalAll)
    KNITRO.KN_set_cb_grad(kc, cb, evalAll)
    KNITRO.KN_set_cb_hess(kc, cb, KNITRO.KN_DENSE_ROWMAJOR, evalAll)

    KNITRO.KN_set_compcons(kc, [KNITRO.KN_CCTYPE_VARVAR], Int32[0], Int32[1])

    KNITRO.KN_set_var_honorbnds(kc,
                                [KNITRO.KN_HONORBNDS_ALWAYS,
                                KNITRO.KN_HONORBNDS_INITPT,
                                KNITRO.KN_HONORBNDS_NO])

    KNITRO.KN_set_var_scalings(kc, Int32[1,2,3], [1.,1,1])
    KNITRO.KN_set_con_scalings(kc, [0.5])
    KNITRO.KN_set_compcon_scalings(kc, [2.])
    KNITRO.KN_set_obj_scaling(kc, 10.)

    status = KNITRO.KN_solve(kc)
    @test status == 0

    # Retrieve derivatives values
    objGrad = KNITRO.KN_get_objgrad_values(kc)
    jac     = KNITRO.KN_get_jacobian_values(kc)
    hess    = KNITRO.KN_get_hessian_values(kc)

    nStatus, objSol, x, lambda_ = KNITRO.KN_get_solution(kc)
    @test nStatus == 0
    @test x ≈ [0., 2.0, 1.98]
    @test objSol ≈ 31.363199 atol=1e-5

    KNITRO.KN_free(kc)
end

@testset "Fourth problem test" begin
    kc = KNITRO.KN_new()

    # START: Some specific parameter settings
    KNITRO.KN_set_param(kc, "presolve", 0)
    KNITRO.KN_set_param(kc, "outlev", 0)
    KNITRO.KN_set_param(kc, "gradopt", 2)
    KNITRO.KN_set_param(kc, "hessopt", 2)
    # END:   Some specific parameter settings

    function evalF_evalGA(kc, cb, evalRequest, evalResult, userParams)
        x = evalRequest.x
        evalRequestCode = evalRequest.evalRequestCode

        if evalRequestCode == KNITRO.KN_RC_EVALFC
            # Evaluate nonlinear objective
            evalResult.obj[1] = x[1]^2 * x[3] + x[2]^3 * x[3]^2
        elseif evalRequestCode == KNITRO.KN_RC_EVALGA
            evalResult.objGrad[1] = 2 * x[1] * x[3]
            evalResult.objGrad[2] = 3 * x[2]^2 * x[3]^2
            evalResult.objGrad[3] = x[1]^2 + 2 * x[2]^3 * x[3]
        else
            return KNITRO.KN_RC_CALLBACK_ERR
        end
        return 0
    end

    # Define objective goal
    objGoal = KNITRO.KN_OBJGOAL_MAXIMIZE
    KNITRO.KN_set_obj_goal(kc, objGoal)

    # Add the variables and set their bounds.
    nV = 3
    KNITRO.KN_add_vars(kc, nV)
    KNITRO.KN_set_var_lobnds(kc, [0, 0.1, 0])
    KNITRO.KN_set_var_upbnds(kc, [0., 2, 2])
    KNITRO.KN_set_var_types(kc, [KNITRO.KN_VARTYPE_CONTINUOUS, KNITRO.KN_VARTYPE_INTEGER, KNITRO.KN_VARTYPE_INTEGER])

    # Define an initial point.
    KNITRO.KN_set_var_primal_init_values(kc, [1, 1, 1.5])
    KNITRO.KN_set_var_dual_init_values(kc, [1., 1, 1, 1])

    # Add the constraints and set their lower bounds.
    nC = 1
    KNITRO.KN_add_cons(kc, nC)
    KNITRO.KN_set_con_lobnds(kc, [0.1])
    KNITRO.KN_set_con_upbnds(kc, [2 * 2 * 0.99])

    # Load quadratic structure x1*x2 for the constraint.
    KNITRO.KN_add_con_quadratic_struct(kc, 0, 1, 2, 1.0)

    # Define callback functions.
    cb = KNITRO.KN_add_objective_callback(kc, evalF_evalGA)
    KNITRO.KN_set_cb_grad(kc, cb, evalF_evalGA)

    # Define complementarity constraints
    KNITRO.KN_set_compcons(kc, [KNITRO.KN_CCTYPE_VARVAR], Int32[0], Int32[1])

    # Set MIP parameters
    KNITRO.KN_set_mip_branching_priorities(kc, Int32[0, 1, 2])
    # not compatible with MPEC constraint as a variable cannot be involved in two different complementarity constraints.
    # KNITRO.KN_set_mip_intvar_strategies(kc, 2, KNITRO.KN_MIP_INTVAR_STRATEGY_MPEC)
    KNITRO.KN_set_mip_node_callback(kc, callback("mip_node"))

    # Set var, con and obj names
    KNITRO.KN_set_var_names(kc, ["myvar1", "myvar2", "myvar3"])
    KNITRO.KN_set_con_names(kc, ["mycon1"])
    KNITRO.KN_set_obj_name(kc, "myobj")

    # Set feasibility tolerances
    KNITRO.KN_set_var_feastols(kc, [0.1, 0.001, 0.1])
    KNITRO.KN_set_con_feastols(kc, [0.1])
    KNITRO.KN_set_compcon_feastols(kc, [0.1])

    # Set finite differences step size
    KNITRO.KN_set_cb_relstepsizes(kc, cb, [0.1, 0.001, 0.1])

    # Solve the problem.
    status = KNITRO.KN_solve(kc)
    # Test for return codes 0 for optimality, and KN_RC_MIP_EXH_FEAS for all nodes explored, assumed optimal
    @test status == 0 || status == KNITRO.KN_RC_MIP_EXH_FEAS

    @test KNITRO.KN_get_mip_number_nodes(kc) >= 1
    @test KNITRO.KN_get_mip_number_solves(kc) >= 1
    @test KNITRO.KN_get_mip_relaxation_bnd(kc) ≈ 31.0495680023
    @test KNITRO.KN_get_mip_lastnode_obj(kc) ≈ 32.0
    @test KNITRO.KN_get_con_values(kc)[1] == 4.0
    @test KNITRO.KN_get_mip_incumbent_obj(kc) ≈ 32.0
    @test KNITRO.KN_get_mip_incumbent_x(kc) == 0.0

    KNITRO.KN_free(kc)
end

@testset "Fifth problem test" begin
    # Test in this environment the setting of user params
    myParams = "stringUserParam"
    kc = KNITRO.KN_new()

    # START: Some specific parameter settings
    KNITRO.KN_set_param(kc, "outlev", 0)
    KNITRO.KN_set_param(kc, "gradopt", 1)
    # END:   Some specific parameter settings

    function evalR(kc, cb, evalRequest, evalResult, userParams)
        x = evalRequest.x
        # Each time we call this callback, the userParams should
        # be as specified.
        @test userParams == myParams
        evalResult.rsd[1] = x[1] * 1.309^x[2] - 2.138
        evalResult.rsd[2] = x[1] * 1.471^x[2] - 3.421
        evalResult.rsd[3] = x[1] * 1.49^x[2] - 3.597
        evalResult.rsd[4] = x[1] * 1.565^x[2] - 4.34
        evalResult.rsd[5] = x[1] * 1.611^x[2] - 4.882
        evalResult.rsd[6] = x[1] * 1.68^x[2] - 5.66
        return 0
    end

    function evalJ(kc, cb, evalRequest, evalResult, userParams)
        x = evalRequest.x
        # Each time we call this callback, the userParams should
        # be as specified.
        @test userParams == myParams
        evalResult.rsdJac[1] = 1.309^x[2]
        evalResult.rsdJac[2] = x[1] * log(1.309) * 1.309^x[2]
        evalResult.rsdJac[3] = 1.471^x[2]
        evalResult.rsdJac[4] = x[1] * log(1.471) * 1.471^x[2]
        evalResult.rsdJac[5] = 1.49^x[2]
        evalResult.rsdJac[6] = x[1] * log(1.49) * 1.49^x[2]
        evalResult.rsdJac[7] = 1.565^x[2]
        evalResult.rsdJac[8] = x[1] * log(1.565) * 1.565^x[2]
        evalResult.rsdJac[9] = 1.611^x[2]
        evalResult.rsdJac[10] = x[1] * log(1.611) * 1.611^x[2]
        evalResult.rsdJac[11] = 1.68^x[2]
        evalResult.rsdJac[12] = x[1] * log(1.68) * 1.68^x[2]
        return 0
    end

    # Add the variables and set their bounds.
    nV = 2
    KNITRO.KN_add_vars(kc, nV)
    KNITRO.KN_set_var_lobnds(kc, [ -1.0, -1.0 ])
    KNITRO.KN_set_var_upbnds(kc, [ 1.0, 1.0 ])
    KNITRO.KN_set_var_primal_init_values(kc, [ 1.0, 5.0 ])

    # Add the residuals
    KNITRO.KN_add_rsds(kc, 6)

    # Define callbacks
    cb = KNITRO.KN_add_lsq_eval_callback(kc, evalR)
    nnzJ = 12
    KNITRO.KN_set_cb_rsd_jac(kc, cb, nnzJ, evalJ,
                            jacIndexRsds=Int32[ 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5 ],
                            jacIndexVars=Int32[ 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1 ])
    KNITRO.KN_set_cb_user_params(kc, cb, myParams)

    # Solve the problem.
    status = KNITRO.KN_solve(kc)
    @test status == 0

    jac = KNITRO.KN_get_rsd_jacobian_values(kc)

    nStatus, objSol, x, lambda_ = KNITRO.KN_get_solution(kc)
    @test nStatus == 0

    @test objSol ≈ 21.5848 atol=1e-3
    @test x ≈ [1., 1.] atol=1e-5

    KNITRO.KN_free(kc)
end

@testset "User callback test (issue #110)" begin
    kc = KNITRO.KN_new()

    KNITRO.KN_set_param(kc, "outlev", 0)
    # Define objective goal
    objGoal = KNITRO.KN_OBJGOAL_MAXIMIZE
    KNITRO.KN_set_obj_goal(kc, objGoal)

    # Add the variables and set their bounds.
    nV = 3
    KNITRO.KN_add_vars(kc, nV)
    KNITRO.KN_set_var_lobnds(kc, [0, 0.1, 0])
    KNITRO.KN_set_var_upbnds(kc, [0., 2, 2])

    # Define an initial point.
    KNITRO.KN_set_var_primal_init_values(kc, [1, 1, 1.5])
    KNITRO.KN_set_var_dual_init_values(kc, [1., 1, 1, 1])

    # Add the constraints and set their lower bounds.
    nC = 1
    KNITRO.KN_add_cons(kc, nC)
    KNITRO.KN_set_con_lobnds(kc, [0.1])
    KNITRO.KN_set_con_upbnds(kc, [2 * 2 * 0.99])

    # Load quadratic structure x1*x2 for the constraint.
    KNITRO.KN_add_con_quadratic_struct(kc, 0, 1, 2, 1.0)

    # Define callback functions.
    cb = KNITRO.KN_add_eval_callback_all(kc, evalAll)
    KNITRO.KN_set_cb_grad(kc, cb, evalAll)
    KNITRO.KN_set_cb_hess(kc, cb, KNITRO.KN_DENSE_ROWMAJOR, evalAll)

    KNITRO.KN_set_compcons(kc, [KNITRO.KN_CCTYPE_VARVAR], Int32[0], Int32[1])

    function newpt_callback(kc, x, lambda_, user_data)
        if KNITRO.KN_get_number_iters(kc) > 1
            return KNITRO.KN_RC_USER_TERMINATION
        end
        return 0
    end

    KNITRO.KN_set_newpt_callback(kc, newpt_callback)

    status = KNITRO.KN_solve(kc)
    @test status == KNITRO.KN_RC_USER_TERMINATION

    nStatus, objSol, x, lambda_ = KNITRO.KN_get_solution(kc)
    @test nStatus == KNITRO.KN_RC_USER_TERMINATION
    @test KNITRO.KN_get_number_iters(kc) == 2

    KNITRO.KN_free(kc)
end

@testset "Handling exception in callbacks" begin
    function eval_kn(kc, cb, evalRequest, evalResult, userParams)
        # Generate exception in callback
        throw(LoadError)
        return 0
    end

    kc = KNITRO.KN_new()
    KNITRO.KN_set_param(kc, "outlev", 0)
    KNITRO.KN_add_vars(kc, 1)
    KNITRO.KN_set_var_primal_init_values(kc, [0.0])
    cb = KNITRO.KN_add_objective_callback(kc, eval_kn)
    nstatus = KNITRO.KN_solve(kc)
    @test nstatus == KNITRO.KN_RC_CALLBACK_ERR
    KNITRO.KN_free(kc)
end

@testset "Knitro evaluation exception" begin
    function eval_kn(kc, cb, evalRequest, evalResult, userParams)
        x = evalRequest.x
        evalResult.obj[1] = sqrt(x[1])
        return 0
    end

    kc = KNITRO.KN_new()
    KNITRO.KN_set_param(kc, "outlev", 0)
    KNITRO.KN_add_vars(kc, 1)
    # Start from a non-evaluable point
    KNITRO.KN_set_var_primal_init_values(kc, [-1.0])
    cb = KNITRO.KN_add_objective_callback(kc, eval_kn)
    nstatus = KNITRO.KN_solve(kc)
    nStatus, objSol, x, lambda_ = KNITRO.KN_get_solution(kc)
    @test x ≈ [0.] atol=1e-5
    KNITRO.KN_free(kc)
end

@testset "Knitro violation information" begin
    if KNITRO.KNITRO_VERSION < v"12.4"
        return 0
    end
    #*------------------------------------------------------------------*
    #*     FUNCTION callbackEvalFC                                      *
    #*------------------------------------------------------------------*
    # The signature of this function matches KNITRO.KN_eval_callback in knitro.h.
    # Only "obj" and "c" are set in the KNITRO.KN_eval_result structure.
    function callbackEvalFC(kc, cb, evalRequest, evalResult, userParams)
        x = evalRequest.x

        # Evaluate nonlinear term in objective
        evalResult.obj[1] = x[1] * x[2] * x[3] * x[4]

        # Evaluate nonlinear terms in constraints
        evalResult.c[1] = x[1] * x[1] * x[1]
        evalResult.c[2] = x[1] * x[1] * x[4]

        return 0
    end

    #*------------------------------------------------------------------*
    #*     main                                                         *
    #*------------------------------------------------------------------*

    # Create a new Knitro solver instance.
    kc = KNITRO.KN_new()
    KNITRO.KN_set_param(kc, "outlev", 0)

    # Initialize Knitro with the problem definition.

    # Add the variables and specify initial values for them.
    # Note: any unset lower bounds are assumed to be
    # unbounded below and any unset upper bounds are
    # assumed to be unbounded above.
    xIndices = KNITRO.KN_add_vars(kc, 4)
    for x in xIndices
        KNITRO.KN_set_var_primal_init_values(kc, x, 0.8)
    end

    # Add the constraints and set the rhs and coefficients
    KNITRO.KN_add_cons(kc, 3)
    KNITRO.KN_set_con_eqbnds(kc, [1.0, 0.0, 0.0])

    # Coefficients for 2 linear terms
    lconIndexCons = Int32[1, 2]
    lconIndexVars = Int32[2, 1]
    lconCoefs = [-1.0, -1.0]
    KNITRO.KN_add_con_linear_struct(kc, lconIndexCons, lconIndexVars, lconCoefs)

    # Coefficients for 2 quadratic terms

    # 1st term:  x1^2 term in c0
    # 2nd term:  x3^2 term in c2
    qconIndexCons = Int32[0, 2]
    qconIndexVars1 = Int32[1, 3]
    qconIndexVars2 = Int32[1, 3]
    qconCoefs = [1.0, 1.0]

    KNITRO.KN_add_con_quadratic_struct(kc, qconIndexCons, qconIndexVars1, qconIndexVars2, qconCoefs)

    # Add callback to evaluate nonlinear(non-quadratic) terms in the model:
    #    x0*x1*x2*x3  in the objective
    #    x0^3         in first constraint c0
    #    x0^2*x3      in second constraint c1
    cb = KNITRO.KN_add_eval_callback(kc, true, Int32[0, 1], callbackEvalFC)

    # Set minimize or maximize(if not set, assumed minimize)
    KNITRO.KN_set_obj_goal(kc, KNITRO.KN_OBJGOAL_MAXIMIZE)

    # Solve the problem.
    #
    # Return status codes are defined in "knitro.h" and described
    # in the Knitro manual.
    nStatus = KNITRO.KN_solve(kc)

    println("Knitro converged with final status = ", nStatus)

    # An example of obtaining solution information.
    nStatus, objSol, x, lambda_ =  KNITRO.KN_get_solution(kc)
    println("  optimal objective value  = ", objSol)
    println("  optimal primal values x  = ", x)
    println("  feasibility violation    = ", KNITRO.KN_get_abs_feas_error(kc))
    println("  KKT optimality violation = ", KNITRO.KN_get_abs_opt_error(kc))

    varbndInfeas, varintInfeas, varviols = KNITRO.KN_get_var_viols(kc, Cint[0,1,2,3])
    println("Variables bound violations = ", varbndInfeas)
    println("Variables integrality violations = ", varintInfeas)
    println("Variables violation values = ", varviols)

    coninfeas, conviols = KNITRO.KN_get_con_viols(kc, Cint[0,1,2])
    println("Constraints bound violations = ", coninfeas)
    println("Constraints violation values = ", conviols)

    println(KNITRO.KN_get_presolve_error(kc))
    @testset "Example HS40 nlp1noderivs" begin
        @test varbndInfeas == [0,0,0,0]
        @test varintInfeas == [0,0,0,0]
        @test varviols ≈ [0.,0.,0.,0.] atol=1e-6
        @test coninfeas == [0,0,0]
        @test conviols ≈ [0.,0.,0.] atol=1e-6
        @test KNITRO.KN_get_abs_feas_error(kc) == max(conviols...)
    end

    # Delete the Knitro solver instance.
    KNITRO.KN_free(kc)
end

@testset "Knitro structural manipulation" begin
    if KNITRO.KNITRO_VERSION < v"12.4"
        return 0
    end
    #*------------------------------------------------------------------*
    #*     FUNCTION callbackEvalFC                                      *
    #*------------------------------------------------------------------*
    # The signature of this function matches KNITRO.KN_eval_callback in knitro.h.
    # Only "obj" and "c" are set in the KNITRO.KN_eval_result structure.
    function callbackEvalFC(kc, cb, evalRequest, evalResult, userParams)
        x = evalRequest.x

        # Evaluate nonlinear term in objective
        evalResult.obj[1] = x[1] * x[2] * x[3] * x[4]

        # Evaluate nonlinear terms in constraints
        evalResult.c[1] = x[1] * x[1] * x[1]
        evalResult.c[2] = x[1] * x[1] * x[4]

        return 0
    end

    #*------------------------------------------------------------------*
    #*     FUNCTION callbackEvalGA                                      *
    #*------------------------------------------------------------------*
    # The signature of this function matches KNITRO.KN_eval_callback in knitro.h.
    # Only "objGrad" and "jac" are set in the KNITRO.KN_eval_result structure.
    function callbackEvalGA(kc, cb, evalRequest, evalResult, userParams)
        x = evalRequest.x

        # Evaluate nonlinear term in objective gradient
        evalResult.objGrad[1] = x[2] * x[3] * x[4]
        evalResult.objGrad[2] = x[1] * x[3] * x[4]
        evalResult.objGrad[3] = x[1] * x[2] * x[4]
        evalResult.objGrad[4] = x[1] * x[2] * x[3]

        # Evaluate nonlinear terms in constraint gradients(Jacobian)
        evalResult.jac[1] = 3.0 * x[1] * x[1] # derivative of x0^3 term  wrt x0
        evalResult.jac[2] = 2.0 * x[1] * x[4] # derivative of x0^2 * x3 term  wrt x0
        evalResult.jac[3] = x[1] * x[1]       # derivative of x0^2 * x3 terms wrt x3

        return 0
    end

    #*------------------------------------------------------------------*
    #*     FUNCTION callbackEvalH                                       *
    #*------------------------------------------------------------------*
    # The signature of this function matches KNITRO.KN_eval_callback in knitro.h.
    # Only "hess" or "hessVec" are set in the KNITRO.KN_eval_result structure.
    function callbackEvalH(kc, cb, evalRequest, evalResult, userParams)
        x = evalRequest.x
        lambda_ = evalRequest.lambda
        # Scale objective component of hessian by sigma
        sigma = evalRequest.sigma

        # Evaluate nonlinear term in the Hessian of the Lagrangian.
        # Note: If sigma=0, some computations can be avoided.
        if sigma > 0.0 # Evaluate the full Hessian of the Lagrangian
            evalResult.hess[1] = lambda_[1] * 6.0 * x[1] + lambda_[2] * 2.0 * x[4]
            evalResult.hess[2] = sigma * x[3] * x[4]
            evalResult.hess[3] = sigma * x[2] * x[4]
            evalResult.hess[4] = sigma * x[2] * x[3] + lambda_[2] * 2.0 * x[1]
            evalResult.hess[5] = sigma * x[1] * x[4]
            evalResult.hess[6] = sigma * x[1] * x[3]
            evalResult.hess[7] = sigma * x[1] * x[2]
        else # sigma=0, do not include objective component
            evalResult.hess[1] = lambda_[1] * 6.0 * x[1] + lambda_[2] * 2.0 * x[4]
            evalResult.hess[2] = 0.0
            evalResult.hess[3] = 0.0
            evalResult.hess[4] = lambda_[2] * 2.0 * x[1]
            evalResult.hess[5] = 0.0
            evalResult.hess[6] = 0.0
            evalResult.hess[7] = 0.0
        end

        return 0
    end

    #*------------------------------------------------------------------*
    #*     main                                                         *
    #*------------------------------------------------------------------*

    # Create a new Knitro solver instance.
    kc = KNITRO.KN_new()
    KNITRO.KN_set_param(kc, "outlev", 6)

    # Initialize Knitro with the problem definition.

    # Add the variables and specify initial values for them.
    # Note: any unset lower bounds are assumed to be
    # unbounded below and any unset upper bounds are
    # assumed to be unbounded above.
    xIndices = KNITRO.KN_add_vars(kc, 4)
    for x in xIndices
        KNITRO.KN_set_var_primal_init_values(kc, x, 0.8)
    end

    # Add the constraints and set the rhs and coefficients
    KNITRO.KN_add_cons(kc, 3)
    KNITRO.KN_set_con_eqbnds(kc, [1.0, 0.0, 0.0])

    # Coefficients for 2 linear terms
    lconIndexCons = Int32[1, 2]
    lconIndexVars = Int32[2, 1]
    lconCoefs = [-1.0, -1.0]
    KNITRO.KN_add_con_linear_struct(kc, lconIndexCons, lconIndexVars, lconCoefs)

    # Coefficients for 2 quadratic terms

    # 1st term:  x1^2 term in c0
    # 2nd term:  x3^2 term in c2
    qconIndexCons = Int32[0, 2]
    qconIndexVars1 = Int32[1, 3]
    qconIndexVars2 = Int32[1, 3]
    qconCoefs = [1.0, 1.0]

    KNITRO.KN_add_con_quadratic_struct(kc, qconIndexCons, qconIndexVars1, qconIndexVars2, qconCoefs)

    # Add callback to evaluate nonlinear(non-quadratic) terms in the model:
    #    x0*x1*x2*x3  in the objective
    #    x0^3         in first constraint c0
    #    x0^2*x3      in second constraint c1
    cb = KNITRO.KN_add_eval_callback(kc, true, Int32[0, 1], callbackEvalFC)

    # Set obj. gradient and nonlinear jac provided through callbacks.
    # Mark objective gradient as dense, and provide non-zero sparsity
    # structure for constraint Jacobian terms.
    cbjacIndexCons = Int32[0, 1, 1]
    cbjacIndexVars = Int32[0, 0, 3]
    KNITRO.KN_set_cb_grad(kc, cb, callbackEvalGA, jacIndexCons=cbjacIndexCons, jacIndexVars=cbjacIndexVars)

    # Set nonlinear Hessian provided through callbacks. Since the
    # Hessian is symmetric, only the upper triangle is provided.
    # The upper triangular Hessian for nonlinear callback structure is:
    #    # lambda0*6*x0 + lambda1*2*x3     x2*x3    x1*x3    x1*x2 + lambda1*2*x0
    #              0                    0      x0*x3         x0*x2
    #                                            0           x0*x1
    #                                                         0
    #(7 nonzero elements)
    cbhessIndexVars1 = Int32[0, 0, 0, 0, 1, 1, 2]
    cbhessIndexVars2 = Int32[0, 1, 2, 3, 2, 3, 3]
    KNITRO.KN_set_cb_hess(kc, cb, length(cbhessIndexVars1), callbackEvalH,
                        hessIndexVars1=cbhessIndexVars1,  hessIndexVars2=cbhessIndexVars2)

    # Set minimize or maximize(if not set, assumed minimize)
    KNITRO.KN_set_obj_goal(kc, KNITRO.KN_OBJGOAL_MAXIMIZE)

    # Solve the initial problem.
    #
    # Return status codes are defined in "knitro.h" and described
    # in the Knitro manual.
    nStatus = KNITRO.KN_solve(kc)

    # An example of obtaining solution information.
    nStatus_origin, objSol_origin, x_origin, lambda_origin =  KNITRO.KN_get_solution(kc)

    # =============== MODIFY PROBLEM AND RE-SOLVE ===========
    # Add 0.5x3 linear term to c2
    KNITRO.KN_add_con_linear_struct(kc, 2, 3, 0.5)
    # Change -x2 to 5x2 in c1
    KNITRO.KN_chg_con_linear_term(kc, 1, 2, 5.0)
    # Now add a new linear constraint x1 + 2x2 + x3 <= 2.5 (c3) and re-solve
    c3 = KNITRO.KN_add_con(kc)
    KNITRO.KN_set_con_upbnd(kc, c3, 2.5)
    KNITRO.KN_add_con_linear_struct(kc, c3, Int32[1,2,3], [1.,2.,1.])

    # Add a constant to the objective
    KNITRO.KN_add_obj_constant(kc, 100.)

    # Tell Knitro to try a "warm-start" since it is starting from the solution
    # of the previous solve, which may be a good initial point for the solution
    # of the slightly modified problem.
    KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_STRAT_WARM_START, KNITRO.KN_STRAT_WARM_START_YES)

    nStatus = KNITRO.KN_solve(kc)
    nStatus, objSol, x, lambda_ =  KNITRO.KN_get_solution(kc)

    # =============== MODIFY PROBLEM BACK TO THE ORIGINAL ONE AND RE-SOLVE AGAIN ===========
    # Remove 0.5x3 from c2
    KNITRO.KN_del_con_linear_term(kc, 2, 3)
    # Change 5x2 back to -x2 in c1
    KNITRO.KN_chg_con_linear_term(kc, 1, 2, -1.0)
    # Remove new constraint c3
    KNITRO.KN_del_con_linear_struct(kc, c3, Int32[1,2,3])
    # Remove the constant in the objective
    KNITRO.KN_del_obj_constant(kc)

    nStatus = KNITRO.KN_solve(kc)
    nStatus, objSol, x, lambda_ =  KNITRO.KN_get_solution(kc)

    # Test that the solution is the same than the original problem
    @testset begin
        @test nStatus == nStatus_origin
        @test objSol ≈ objSol_origin atol=1e-6
        @test x ≈ x_origin atol=1e-6
        @test lambda_[1:length(lambda_origin)] ≈ lambda_origin atol=1e-1
    end

    # Delete the Knitro solver instance.
    KNITRO.KN_free(kc)
end


text:
using Distributed
@everywhere using Mamba

## Data
dogs = Dict{Symbol, Any}(
  :Y =>
    [0 0 1 0 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
     0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1
     0 0 0 0 0 1 1 0 1 1 0 0 1 1 0 1 0 1 1 1 1 1 1 1 1
     0 1 1 0 0 1 1 1 1 0 1 0 1 0 1 1 1 1 1 1 1 1 1 1 1
     0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
     0 0 0 0 0 0 1 1 1 1 0 0 1 0 1 1 1 1 1 1 1 1 1 1 1
     0 0 0 0 0 1 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1
     0 0 0 0 0 0 0 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1
     0 0 0 0 0 1 0 1 0 1 1 0 1 0 0 0 1 1 1 1 1 0 1 1 0
     0 0 0 0 1 0 0 1 1 0 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1
     0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
     0 0 0 0 0 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1
     0 0 0 1 1 0 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
     0 0 0 0 1 0 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
     0 0 0 1 0 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
     0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
     0 1 0 1 0 0 0 1 0 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1
     0 0 0 0 1 0 1 0 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1
     0 1 0 0 0 0 1 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
     0 0 0 0 1 1 0 1 0 1 1 0 1 0 1 1 1 1 1 1 1 1 1 1 1
     0 0 0 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
     0 0 1 0 1 0 1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1
     0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
     0 0 0 0 0 0 0 0 1 1 1 0 1 0 0 0 1 1 0 1 1 1 1 1 1
     0 0 0 0 0 0 1 1 0 1 1 1 0 1 0 1 1 1 1 1 1 1 1 1 1
     0 0 1 0 1 1 1 0 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1
     0 0 0 0 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
     0 0 0 1 0 1 0 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1
     0 0 0 0 1 1 0 0 1 1 1 0 1 0 1 0 1 0 1 1 1 1 1 1 1
     0 0 0 0 1 1 1 1 1 1 0 1 0 1 1 1 1 1 1 1 1 1 1 1 1]
)
dogs[:Dogs] = size(dogs[:Y], 1)
dogs[:Trials] = size(dogs[:Y], 2)

dogs[:xa] = mapslices(cumsum, dogs[:Y], dims=2)
dogs[:xs] = mapslices(x -> collect(1:25) .- x, dogs[:xa], dims=2)
dogs[:y] = 1 .- dogs[:Y][:, 2:25]

## Model Specification
model = Model(

  y = Stochastic(2,
    (Dogs, Trials, alpha, xa, beta, xs) ->
      UnivariateDistribution[(
        p = exp(alpha * xa[i, j] + beta * xs[i, j]);
        Bernoulli(p)) for i in 1:Dogs, j in 1:Trials-1
      ],
    false
  ),

  alpha = Stochastic(
    () -> Truncated(Flat(), -Inf, -1e-5)
  ),

  A = Logical(
    alpha -> exp(alpha)
  ),

  beta = Stochastic(
    () -> Truncated(Flat(), -Inf, -1e-5)
  ),

  B = Logical(
    beta -> exp(beta)
  )

)

## Initial Values
inits = [
  Dict(:y => dogs[:y], :alpha => -1, :beta => -1),
  Dict(:y => dogs[:y], :alpha => -2, :beta => -2)
]

## Sampling Scheme
scheme = [Slice([:alpha, :beta], 1.0)]
setsamplers!(model, scheme)

## MCMC Simulations
sim = mcmc(model, dogs, inits, 10000, burnin=2500, thin=2, chains=2)
describe(sim)

ext: jl


text:
function psoptimize(f, Opts::PSOOptions)

    # Initialize Parameters
    Time0 = time()
    Iterations = 0

    StallTime0 = Time0
    StallIterations = 0
    FStall = Inf

    # Initialize Swarm with Uniform Distribution
    swarm = initswarmuniform(f, Opts)

    # Compute Minimum Neighborhood Size
    MinNeighborhoodSize = max(2, floor(Opts.SwarmSize * Opts.MinNeighborsFraction))

    # Begin loop
    ExitFlag = 0
    while ExitFlag == 0

        # Update Iteration Counter
        Iterations += 1

        # Update Swarm Velocities
        updatevelocities!(swarm)

        # Update Positions
        step!(swarm)

        # Enforce Bounds
        if !(all(isinf.(Opts.LowerBounds)) && all(isinf.(Opts.UpperBounds)))
            enforcebnds!(swarm, Opts)
        end

        # Evaluate Objective Function
        feval!(f, swarm, Opts)

        # Update Global Best
        flag = setgbest!(swarm)

        # Update Stall Counter and Neighborhood
        if flag
            swarm.c = max(0, swarm.c - 1)
            swarm.N = MinNeighborhoodSize
        else
            swarm.c += 1
            swarm.N = min(swarm.N + MinNeighborhoodSize, Opts.SwarmSize - 1)
        end

        # Update Inertia
        if swarm.c < 2
            swarm.W *= 2.0
        elseif swarm.c > 5
            swarm.W /= 2.0
        end

        # Ensure Inertia is in InertiaRange
        if swarm.W > Opts.InertiaRange[2]
            swarm.W = Opts.InertiaRange[2]
        elseif swarm.W < Opts.InertiaRange[1]
            swarm.W = Opts.InertiaRange[1]
        end

        # Track Stalling
        if FStall - swarm.b > Opts.FunctionTolerance
            FStall = swarm.b
            StallIterations = 0
            StallTime0 = time()
        else
            StallIterations += 1
        end

        # Stopping Criteria
        if StallIterations >= Opts.MaxStallIterations
            ExitFlag = 1
        elseif Iterations >= Opts.MaxIterations
            ExitFlag = 2
        elseif swarm.b <= Opts.ObjectiveLimit
            ExitFlag = 3
        elseif time() - StallTime0 >= Opts.MaxStallTime
            ExitFlag = 4
        elseif time() - Time0 >= Opts.MaxTime
            ExitFlag = 5
        end

        # Output Status
        if Opts.Display == "iter" && Iterations % Opts.DisplayInterval == 0
            printstatus(swarm, time() - Time0, Iterations, StallIterations)
        end
    end

    if Opts.HybridFcn
        if Opts.Display == "iter"
            println("Beginning optimization with hybrid function...")
        end

        res = optimize(f, MVector{Opts.NDims}(swarm.d), Opts.HybridOptimizer, Opts.HybridOptimizerOpts; autodiff = Opts.HybridAutoDiff)
        if Optim.minimum(res) < swarm.b
            fbest = Optim.minimum(res)
            xbest = SVector{Opts.NDims}(Optim.minimizer(res))
        else
            fbest = swarm.b
            xbest = swarm.d
        end

        return PSOSolution(swarm, fbest, xbest, Iterations, time() - Time0, ExitFlag, res, swarm.b, swarm.d)

    else
        return PSOSolution(swarm, swarm.b, swarm.d, Iterations, time() - Time0, ExitFlag)
    end
end

max_stars_repo_name: GrantHecht/JuSwarm

max_forks_repo_name: GrantHecht/JuSwarm


text:
using ITensors, Test

@testset "not" begin
  i = Index(2, "i")
  j = Index(2, "j")
  k = Index(2, "k")

  A = randomITensor(i, j, k')

  Ap = prime(A, not("j"))

  @test hassameinds(Ap, (i', j, k''))

  Ap = prime(A; tags=!ts"j")

  @test hassameinds(Ap, (i', j, k''))

  At = addtags(A, "x", not("k"))

  @test hassameinds(At, (addtags(i, "x"), addtags(j, "x"), k'))

  Ap2 = prime(A, 2, not(i))

  @test hassameinds(Ap2, (i, j'', k'''))

  Ap2 = prime(A, 2; inds=!i)

  @test hassameinds(Ap2, (i, j'', k'''))

  Ap3 = prime(A, 3, not(i, k'))

  @test hassameinds(Ap3, (i, j''', k'))

  Ap3 = prime(A, 3, !(i, k'))

  @test hassameinds(Ap3, (i, j''', k'))

  At2 = settags(A, "y", not(IndexSet(j, k')))

  @test hassameinds(At2, (settags(i, "y"), j, k'))

  At2 = settags(A, "y"; inds=!IndexSet(j, k'))

  @test hassameinds(At2, (settags(i, "y"), j, k'))

  B = filterinds(A; plev=not(0))

  @test hassameinds(B, (k',))

  @test_throws MethodError !"x"

  @test_throws MethodError !1
end

nothing


text:
struct LinearAnalysis
    options::Dict{Symbol, Any}
end

LinearAnalysis(; kwargs...) = LinearAnalysis(Dict{Symbol, Any}(kwargs))

add_intercept_column(x::AbstractVector{T}) where {T} = [ones(T, length(x)) x]

# TODO: add multidimensional version
function (l::LinearAnalysis)(le::Entry)
    return splitapply(le) do entry
        x, y = entry.positional
        weights = get(entry.named, :weights, similar(x, 0))
        npoints = get(l.options, :npoints, 200)
        interval = get(l.options, :interval, length(weights) > 0 ? nothing : :confidence)
        dropcollinear = false
        lin_model = try
            GLM.lm(add_intercept_column(x), collect(y); wts=weights, dropcollinear)
        catch e
            @warn "linear fit not possible"
            nothing
        end
        isnothing(lin_model) && return Entry[]
        x̂ = range(extrema(x)..., length=npoints)
        pred = GLM.predict(lin_model, add_intercept_column(x̂); interval)
        if !isnothing(interval)
            ŷ, lower, upper = map(vec, pred) # GLM prediction returns matrices
            default_plottype = LinesFill
            positional, named = (x̂, ŷ), (; lower, upper)
        else
            ŷ = vec(pred) # GLM prediction returns matrix
            default_plottype = Lines
            positional, named = (x̂, ŷ), (;)
        end
        plottype = AbstractPlotting.plottype(entry.plottype, default_plottype)
        return Entry(entry; plottype, positional, named)
    end
end

"""
    linear(; interval)

Compute a linear fit of `y ~ 1 + x`. An optional named mapping `weights` determines the weights.
Use `interval` to specify what type of interval the shaded band should represent.
Valid values of interval are `:confidence` delimiting the uncertainty of the predicted
relationship, and `:prediction` delimiting estimated bounds for new data points.
"""
linear(; kwargs...) = Layer((LinearAnalysis(; kwargs...),))


text:
using LinearAlgebra
using SpecialFunctions
using OffsetArrays

const Γ = gamma

# Coefficients in the Jacobi polynomial recurrence relations.
a(α, β, n) = 2/(2n+α+β) * √(n * (n+α+β) * (n+α) * (n+β) / (2n+α+β-1) / (2n+α+β+1))
b(α, β, n) = -(α^2 - β^2) / (2n+α+β) / (2n+α+β+2)

" Evaluates the Jacobi polynomial Pₙᵅᵝ(x). "
function jacobi(x, α, β, n::Int)
    Pᵅᵝ = n <= 1 ? OffsetArray(zeros(2), 0:1) : OffsetArray(zeros(n+1), 0:n)
    Pᵅᵝ[0] = √(2.0^-(α+β+1) * Γ(α+β+2) / Γ(α+1) / Γ(β+1))
    Pᵅᵝ[1] = Pᵅᵝ[0]/2 * √((α+β+3) / (α+1) / (β+1)) * ((α+β+2)*x + α - β)
    for n′ in 1:n-1
        Pᵅᵝ[n′+1] = ((x - b(α,β,n′)) * Pᵅᵝ[n′] - a(α,β,n′) * Pᵅᵝ[n′-1]) / a(α, β, n′+1)
    end
    return Pᵅᵝ[n]
end

" Evaluates the derivative of the Jacobi polynomial Pₙᵅᵝ(x). "
∂jacobi(x, α, β, n::Int) = n == 0 ? 0 : √(n * (n+α+β+1)) * jacobi(x, α+1, β+1, n-1)

" Guassian quadrature points and weights for the Jacobi polynomial Pₙᵅᵝ. "
function jacobi_gauss_quadrature(α, β, N)
    N == 0 && return [(α-β) / (α+β+2)], [2]

    # Form symmetric matrix from recurrence.
    dv = OffsetArray(zeros(N+1), 0:N)  # diagonal vector
    ev = OffsetArray(zeros(N+1), 0:N)  # sub/super-diagonal vector

    for n in 0:N
        dv[n] = b(α, β, n)
        ev[n] = a(α, β, n)
    end

    # Create full matrix combining the two.
    # Need to pass arrays that are not offset.
    J = SymTridiagonal(dv[0:N], ev[1:N])
    (α + β) ≈ 0 && (J[1, 1] = 0)

    # Compute quadrature points and weights by eigenvalue solve.
    x, V = eigen(J)
    w = @. V[1, :]^2 * 2^(α+β+1) / (α+β+1)
    @. w *= factorial(α) * factorial(β) / factorial(α+β)

    return x, w
end

" Guass-Labatto quadrature points for the Jacobi Polynomial Pₙᵅᵝ. "
function jacobi_gauss_lobatto(α, β, N)
    N == 0 && error("What are you doing? Gauss-Lobatto points only make sense if N >= 1.")
    N == 1 && return [-1, 1]

    x = zeros(N+1)
    x[1], x[N+1] = -1, 1

    x_GQ, _ = jacobi_gauss_quadrature(α+1, β+1, N-2)
    x[2:N] .= x_GQ

    return x
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Rotations from an Earth-Centered Inertial (ECI) reference frame to another
#   ECI reference frame.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export r_eci_to_eci

"""
    r_ecef_to_eci([T,] ECIo, ECIf, JD_UTC::Number [, eop_data])
    r_ecef_to_eci([T,] ECIo, JD_UTCo::Number, ECIf, JD_UTCf::Number [, eop_data])

Compute the rotation from an Earth-Centered Inertial (`ECI`) reference frame to
another ECI reference frame. If the origin and destination frame contain only
one *of date* frame, then the first signature is used and `JD_UTC` is the epoch
of this frame. On the other hand, if the origin and destination frame contain
two *of date* frame`¹`, e.g. TOD => MOD, then the second signature must be used
in which `JD_UTCo` is the epoch of the origin frame and `JD_UTCf` is the epoch
of the destination frame.

The rotation description that will be used is given by `T`, which can be `DCM`
or `Quaternion`. The origin ECI frame is selected by the input `ECIo` and the
destination ECI frame is selected by the input `ECIf`. The model used to compute
the rotation is specified by the selection of the origin and destination frames.
Currently, there are two models supported: IAU-76/FK5 and IAU-2006 with 2010
conventions (CIO and equinox approaches).

`¹`: TEME is an *of date* frame.

# Rotation description

The rotations that aligns the origin ECI frame with the destination ECI frame
can be described by Direction Cosine Matrices or Quaternions. This is selected
by the parameter `T`.

The possible values are:

- `DCM`: The rotation will be described by a Direction Cosine Matrix.
- `Quaternion`: The rotation will be described by a Quaternion.

If no value is specified, then it falls back to `DCM`.

# Conversion model

The model that will be used to compute the rotation is automatically inferred
given the selection of the origin and destination frames. **Notice that mixing
IAU-76/FK5 and IAU-2006/2010 frames is not supported.**

# ECI Frame

The supported ECI frames for both origin `ECIo` and destination `ECIf` are:

- `TEME()`: ECI will be selected as the True Equator Mean Equinox (TEME)
    reference frame.
- `TOD()`: ECI will be selected as the True of Date (TOD).
- `MOD()`: ECI will be selected as the Mean of Date (MOD).
- `J2000()`: ECI will be selected as the J2000 reference frame.
- `GCRF()`: ECI will be selected as the Geocentric Celestial Reference Frame
    (GCRF).
- `CIRS()`: ECEF will be selected as the Celestial Intermediate Reference System
    (CIRS).
- `ERS()`: ECI will be selected as the Earth Reference System (ERS).
- `MOD06()`: ECI will be selected as the Mean of Date (MOD) according to the
    definition in IAU-2006/2010 theory.
- `MJ2000()`: ECI will be selected as the J2000 mean equatorial frame (MJ2000).

!!! note
    The frames `MOD()` and `MOD06()` are virtually the same. However, we
    selected different names to make clear which theory are being used since
    mixing transformation between frames from IAU-76/FK5 and IAU-2006/2010 must
    be performed with caution.

# EOP Data

The conversion between the frames depends on EOP Data (see
[`get_iers_eop`](@ref) and [`read_iers_eop`](@ref)). If IAU-76/FK5 model is
used, then the type of `eop_data` must be [`EOPData_IAU1980`](@ref). Otherwise,
if IAU-2006/2010 model is used, then the type of `eop_data` must be
[`EOPData_IAU2000A`](@ref). The following table shows the requirements for EOP
data given the selected frames.

|   Model                     |   ECIo   |   ECIf   |    EOP Data   | Function Signature |
|:----------------------------|:---------|:---------|:--------------|:-------------------|
| IAU-76/FK5                  | `GCRF`   | `J2000`  | EOP IAU1980   | First              |
| IAU-76/FK5                  | `GCRF`   | `MOD`    | EOP IAU1980   | First              |
| IAU-76/FK5                  | `GCRF`   | `TOD`    | EOP IAU1980   | First              |
| IAU-76/FK5                  | `GCRF`   | `TEME`   | EOP IAU1980   | First              |
| IAU-76/FK5                  | `J2000`  | `GCRF`   | EOP IAU1980   | First              |
| IAU-76/FK5                  | `J2000`  | `MOD`    | Not required  | First              |
| IAU-76/FK5                  | `J2000`  | `TOD`    | Not required  | First              |
| IAU-76/FK5                  | `J2000`  | `TEME`   | Not required  | First              |
| IAU-76/FK5                  | `MOD`    | `GCRF`   | EOP IAU1980   | First              |
| IAU-76/FK5                  | `MOD`    | `J2000`  | Not required  | First              |
| IAU-76/FK5                  | `MOD`    | `TOD`    | Not required  | Second             |
| IAU-76/FK5                  | `MOD`    | `TEME`   | Not required  | Second             |
| IAU-76/FK5                  | `TOD`    | `GCRF`   | EOP IAU1980   | First              |
| IAU-76/FK5                  | `TOD`    | `J2000`  | Not required  | First              |
| IAU-76/FK5                  | `TOD`    | `MOD`    | Not required  | Second             |
| IAU-76/FK5                  | `TOD`    | `TEME`   | Not required  | Second             |
| IAU-76/FK5                  | `TEME`   | `GCRF`   | EOP IAU1980   | First              |
| IAU-76/FK5                  | `TEME`   | `J2000`  | Not required  | First              |
| IAU-76/FK5                  | `TEME`   | `MOD`    | Not required  | Second             |
| IAU-76/FK5                  | `TEME`   | `TOD`    | Not required  | Second             |
| IAU-2006/2010 CIO-based     | `GCRF`   | `CIRS`   | Not required¹ | First              |
| IAU-2006/2010 CIO-based     | `CIRS`   | `CIRS`   | Not required¹ | Second             |
| IAU-2006/2010 Equinox-based | `GCRF`   | `MJ2000` | Not required  | First²             |
| IAU-2006/2010 Equinox-based | `GCRF`   | `MOD06`  | Not required  | First              |
| IAU-2006/2010 Equinox-based | `GCRF`   | `ERS`    | Not required³ | First              |
| IAU-2006/2010 Equinox-based | `MJ2000` | `GCRF`   | Not required  | First²             |
| IAU-2006/2010 Equinox-based | `MJ2000` | `MOD06`  | Not required  | First              |
| IAU-2006/2010 Equinox-based | `MJ2000` | `ERS`    | Not required³ | First              |
| IAU-2006/2010 Equinox-based | `MOD06`  | `GCRF`   | Not required  | First              |
| IAU-2006/2010 Equinox-based | `MOD06`  | `MJ2000` | Not required  | First              |
| IAU-2006/2010 Equinox-based | `MOD06`  | `ERS`    | Not required³ | First              |
| IAU-2006/2010 Equinox-based | `ERS`    | `GCRF`   | Not required³ | First              |
| IAU-2006/2010 Equinox-based | `ERS`    | `MJ2000` | Not required³ | First              |
| IAU-2006/2010 Equinox-based | `ERS`    | `MOD06`  | Not required³ | First              |

`¹`: In this case, the terms that account for the free-core nutation and time
dependent effects of the Celestial Intermediate Pole (CIP) position with respect
to the GCRF will not be available, reducing the precision.

`²`: The transformation between GCRF and MJ2000 is a constant rotation matrix
called bias. Hence, the date does not modify it. However, this signature was
kept to avoid complications in the API.

`³`: In this case, the terms that corrects the nutation in obliquity and in
longitude due to the free core nutation will not be available, reducing the
precision.

## MOD and TOD

In this function, if EOP corrections are not provided, then MOD and TOD frames
will be computed considering the original IAU-76/FK5 theory. Otherwise, the
corrected frame will be used.

# Returns

The rotation description represented by `T` that rotates the origin ECI
reference frame into alignment with the destination ECI reference frame.

# Examples

```julia-repl
julia> eop_IAU1980 = get_iers_eop(Val(:IAU1980));

julia> r_eci_to_eci(DCM, GCRF(), J2000(), date_to_jd(1986, 6, 19, 21, 35, 0), eop_IAU1980)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
  1.0          -4.71326e-12   1.53474e-9
  4.71332e-12   1.0          -3.53979e-9
 -1.53474e-9    3.53979e-9    1.0

julia> r_eci_to_eci(Quaternion, TEME(), GCRF(), date_to_jd(1986, 6, 19, 21, 35, 0), eop_IAU1980)
Quaternion{Float64}:
  + 0.999999 + 1.83013e-5⋅i + 0.000665304⋅j + 0.000665304⋅k

julia> r_eci_to_eci(TOD(), date_to_jd(1986,6,19,21,35,0), TOD(), date_to_jd(1987,5,19,3,0,0), eop_IAU1980)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 1.0          -0.000224088  -9.73787e-5
 0.000224087   1.0          -5.80065e-6
 9.738e-5      5.77883e-6    1.0

julia> r_eci_to_eci(Quaternion, TOD(), JD_J2000, MOD(), JD_J2000, eop_IAU1980)
Quaternion{Float64}:
  + 1.0 - 1.40025e-5⋅i + 1.34736e-5⋅j + 1.34736e-5⋅k

julia> r_eci_to_eci(J2000(), TEME(), date_to_jd(1986,6,19,21,35,0))
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
  0.999995    0.0030265    0.00133055
 -0.00302645  0.999995    -3.86125e-5
 -0.00133066  3.45854e-5   0.999999

julia> eop_IAU2000A = get_iers_eop(Val(:IAU2000A));

julia> r_eci_to_eci(CIRS(), GCRF(), date_to_jd(1986,6,19,21,35,0), eop_IAU2000A)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 0.999999     3.88389e-8  -0.00133066
 7.18837e-9   1.0          3.45897e-5
 0.00133066  -3.45897e-5   0.999999

julia> r_eci_to_eci(Quaternion, CIRS(), GCRF(), date_to_jd(1986,6,19,21,35,0), eop_IAU2000A)
Quaternion{Float64}:
  + 1.0 + 1.72949e-5⋅i + 0.000665332⋅j + 0.000665332⋅k
```
"""
@inline function r_eci_to_eci(
    T_ECIo::T_ECIs,
    T_ECIf::T_ECIs,
    JD_UTC::Number,
    eop_data::EOPData_IAU1980
)
    return r_eci_to_eci(DCM, T_ECIo, T_ECIf, JD_UTC, eop_data)
end

@inline function r_eci_to_eci(
    T_ECIo::T_ECIs_IAU_2006,
    T_ECIf::T_ECIs_IAU_2006,
    JD_UTC::Number,
    eop_data::EOPData_IAU2000A
)
    return r_eci_to_eci(DCM, T_ECIo, T_ECIf, JD_UTC, eop_data)
end

# Specializations for those cases in which we have two *of dates* frames.
@inline function r_eci_to_eci(
    T_ECIo::T_ECIs_of_date,
    JD_UTCo::Number,
    T_ECIf::T_ECIs_of_date,
    JD_UTCf::Number,
    eop_data::EOPData_IAU1980
)
    return r_eci_to_eci(DCM, T_ECIo, JD_UTCo, T_ECIf, JD_UTCf, eop_data)
end

@inline function r_eci_to_eci(
    T_ECIo::Val{:CIRS},
    JD_UTCo::Number,
    T_ECIf::Val{:CIRS},
    JD_UTCf::Number,
    eop_data::EOPData_IAU2000A
)
    return r_eci_to_eci(DCM, T_ECIo, JD_UTCo, T_ECIf, JD_UTCf, eop_data)
end

@inline function r_eci_to_eci(
    T_ECIo::T_ECIs_IAU_2006_Equinox_of_date,
    JD_UTCo::Number,
    T_ECIf::T_ECIs_IAU_2006_Equinox_of_date,
    JD_UTCf::Number,
    eop_data::EOPData_IAU2000A
)
    return r_eci_to_eci(DCM, T_ECIo, JD_UTCo, T_ECIf, JD_UTCf, eop_data)
end

# Specializations for those cases that EOP Data is not needed.
@inline function r_eci_to_eci(
    T_ECIo::Val{:J2000},
    T_ECIf::Union{Val{:MOD}, Val{:TOD}, Val{:TEME}},
    JD_UTC::Number
)
    return r_eci_to_eci(DCM, T_ECIo, T_ECIf, JD_UTC)
end

@inline function r_eci_to_eci(
    T_ECIo::Union{Val{:MOD}, Val{:TOD}, Val{:TEME}},
    T_ECIf::Val{:J2000},
    JD_UTC::Number
)
    return r_eci_to_eci(DCM, T_ECIo, T_ECIf, JD_UTC)
end

@inline function r_eci_to_eci(
    T_ECIo::T_ECIs_of_date,
    JD_UTCo::Number,
    T_ECIf::T_ECIs_of_date,
    JD_UTCf::Number
)
    return r_eci_to_eci(DCM, T_ECIo, JD_UTCo, T_ECIf, JD_UTCf)
end

@inline function r_eci_to_eci(
    T_ECIo::T_ECIs_IAU_2006,
    T_ECIf::T_ECIs_IAU_2006,
    JD_UTC::Number
)
    return r_eci_to_eci(DCM, T_ECIo, T_ECIf, JD_UTC)
end

 @inline function r_eci_to_eci(
    T_ECIo::Val{:CIRS},
    JD_UTCo::Number,
    T_ECIf::Val{:CIRS},
    JD_UTCf::Number
 )
    return r_eci_to_eci(DCM, T_ECIo, JD_UTCo, T_ECIf, JD_UTCf)
end

@inline function r_eci_to_eci(
    T_ECIo::T_ECIs_IAU_2006_Equinox_of_date,
    T_ECIf::T_ECIs_IAU_2006_Equinox_of_date,
    JD_UTC::Number
)
    return r_eci_to_eci(DCM, T_ECIo, T_ECIf, JD_UTC)
end

################################################################################
#                                  IAU-76/FK5
################################################################################

#                                GCRF <=> J2000
# ==============================================================================

function r_eci_to_eci(
    T::T_ROT,
    ::Val{:GCRF},
    ::Val{:J2000},
    JD_UTC::Number,
    eop_data::EOPData_IAU1980
)
    milliarcsec_to_rad = π / 648000000

    # Get the time in UT1 and TT.
    JD_UT1 = jd_utc_to_ut1(JD_UTC, eop_data)
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Get the EOP data related to the desired epoch.
    δΔϵ_1980 = eop_data.dEps(JD_UTC) * milliarcsec_to_rad
    δΔψ_1980 = eop_data.dPsi(JD_UTC) * milliarcsec_to_rad

    # In this case, we need to convert GCRF back to PEF and then convert to
    # J2000, which is the same conversion from PEF to GCRF **without** the EOP
    # data.
    #
    # TODO: Can I simplify the rotation from TOD with corrections to TOD without
    # corrections?
    r_MOD_GCRF  = r_gcrf_to_mod_fk5(T, JD_TT)
    r_PEF_MOD   = r_mod_to_pef_fk5(T, JD_UT1, JD_TT, δΔϵ_1980, δΔψ_1980)
    r_MOD_PEF   = r_pef_to_mod_fk5(T, JD_UT1, JD_TT, 0, 0)
    r_J2000_MOD = inv_rotation(r_MOD_GCRF)

    # Compose the full rotation.
    return compose_rotation(r_MOD_GCRF, r_PEF_MOD, r_MOD_PEF, r_J2000_MOD)
end

@inline function r_eci_to_eci(
    T::T_ROT,
    T_ECIo::Val{:J2000},
    T_ECIf::Val{:GCRF},
    JD_UTC::Number,
    eop_data::EOPData_IAU1980
)
    return inv_rotation(r_eci_to_eci(T, T_ECIf, T_ECIo, JD_UTC, eop_data))
end

#                                 GCRF <=> MOD
# ==============================================================================

function r_eci_to_eci(
    T::T_ROT,
    ::Val{:GCRF},
    ::Val{:MOD},
    JD_UTC::Number,
    eop_data::EOPData_IAU1980
)
    # Get the time in TT.
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Return the rotation.
    return r_gcrf_to_mod_fk5(T, JD_TT)
end

@inline function r_eci_to_eci(
    T::T_ROT,
    T_ECIo::Val{:MOD},
    T_ECIf::Val{:GCRF},
    JD_UTC::Number,
    eop_data::EOPData_IAU1980
)
    return inv_rotation(r_eci_to_eci(T, T_ECIf, T_ECIo, JD_UTC, eop_data))
end

#                                 GCRF <=> TOD
# ==============================================================================

function r_eci_to_eci(
    T::T_ROT,
    ::Val{:GCRF},
    ::Val{:TOD},
    JD_UTC::Number,
    eop_data::EOPData_IAU1980
)
    milliarcsec_to_rad = π / 648000000

    # Get the time in TT.
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Get the EOP data related to the desired epoch.
    δΔϵ_1980 = eop_data.dEps(JD_UTC) * milliarcsec_to_rad
    δΔψ_1980 = eop_data.dPsi(JD_UTC) * milliarcsec_to_rad

    # Return the rotation.
    r_MOD_GCRF = r_gcrf_to_mod_fk5(T, JD_TT)
    r_TOD_MOD  = r_mod_to_tod_fk5(T, JD_TT, δΔϵ_1980, δΔψ_1980)

    # Compose the full rotation.
    return compose_rotation(r_MOD_GCRF, r_TOD_MOD)
end

@inline function r_eci_to_eci(
    T::T_ROT,
    T_ECIo::Val{:TOD},
    T_ECIf::Val{:GCRF},
    JD_UTC::Number,
    eop_data::EOPData_IAU1980
)
    return inv_rotation(r_eci_to_eci(T, T_ECIf, T_ECIo, JD_UTC, eop_data))
end

#                                GCRF <=> TEME
# ==============================================================================

function r_eci_to_eci(
    T::T_ROT,
    ::Val{:GCRF},
    ::Val{:TEME},
    JD_UTC::Number,
    eop_data::EOPData_IAU1980
)
    milliarcsec_to_rad = π / 648000000

    # Get the time in TT.
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Get the EOP data related to the desired epoch.
    δΔϵ_1980 = eop_data.dEps(JD_UTC) * milliarcsec_to_rad
    δΔψ_1980 = eop_data.dPsi(JD_UTC) * milliarcsec_to_rad

    # Return the rotation.
    r_MOD_GCRF = r_gcrf_to_mod_fk5(T, JD_TT)
    r_TEME_MOD = r_mod_to_teme(T, JD_TT, δΔϵ_1980, δΔψ_1980)

    # Compose the full rotation.
    return compose_rotation(r_MOD_GCRF, r_TEME_MOD)
end

@inline function r_eci_to_eci(
    T::T_ROT,
    T_ECIo::Val{:TEME},
    T_ECIf::Val{:GCRF},
    JD_UTC::Number,
    eop_data::EOPData_IAU1980
)
    return inv_rotation(r_eci_to_eci(T, T_ECIf, T_ECIo, JD_UTC, eop_data))
end

#                                J2000 <=> MOD
# ==============================================================================

function r_eci_to_eci(
    T::T_ROT,
    ::Val{:J2000},
    ::Val{:MOD},
    JD_UTC::Number,
    eop_data::EOPData_IAU1980
)
    milliarcsec_to_rad = π / 648000000

    # Get the time in UT1 and TT.
    JD_UT1 = jd_utc_to_ut1(JD_UTC, eop_data)
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Get the EOP data related to the desired epoch.
    δΔϵ_1980 = eop_data.dEps(JD_UTC) * milliarcsec_to_rad
    δΔψ_1980 = eop_data.dPsi(JD_UTC) * milliarcsec_to_rad

    # In this case, we need to convert J2000 back to PEF and then convert to
    # MOD. This is necessary because we need to apply EOP corrections to convert
    # to MOD and just a `r_gcrf_to_mod_fk5` would yield a frame according to the
    # original IAU-76/FK5 theory.
    #
    # TODO: Can I simplify this rotation?
    r_MOD_J2000 = r_gcrf_to_mod_fk5(T, JD_TT)
    r_PEF_MOD   = r_mod_to_pef_fk5(T, JD_UT1, JD_TT, 0, 0)
    r_MOD_PEF   = r_pef_to_mod_fk5(T, JD_UT1, JD_TT, δΔϵ_1980, δΔψ_1980)

    # Compose the full rotation.
    return compose_rotation(r_MOD_J2000, r_PEF_MOD, r_MOD_PEF)
end

function r_eci_to_eci(T::T_ROT, ::Val{:J2000}, ::Val{:MOD}, JD_UTC::Number)
    # Get the time in TT.
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Compute and return the rotation.
    r_MOD_J2000 = r_gcrf_to_mod_fk5(T, JD_TT)

    return r_MOD_J2000
end

@inline function r_eci_to_eci(
    T::T_ROT,
    T_ECIo::Val{:MOD},
    T_ECIf::Val{:J2000},
    JD_UTC::Number,
    eop_data::EOPData_IAU1980
)
    return inv_rotation(r_eci_to_eci(T, T_ECIf, T_ECIo, JD_UTC, eop_data))
end

@inline function r_eci_to_eci(
    T::T_ROT,
    T_ECIo::Val{:MOD},
    T_ECIf::Val{:J2000},
    JD_UTC::Number
)
    return inv_rotation(r_eci_to_eci(T, T_ECIf, T_ECIo, JD_UTC))
end

#                                J2000 <=> TOD
# ==============================================================================

function r_eci_to_eci(
    T::T_ROT,
    ::Val{:J2000},
    ::Val{:TOD},
    JD_UTC::Number,
    eop_data::EOPData_IAU1980
)
    milliarcsec_to_rad = π / 648000000

    # Get the time in UT1 and TT.
    JD_UT1 = jd_utc_to_ut1(JD_UTC, eop_data)
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Get the EOP data related to the desired epoch.
    δΔϵ_1980 = eop_data.dEps(JD_UTC) * milliarcsec_to_rad
    δΔψ_1980 = eop_data.dPsi(JD_UTC) * milliarcsec_to_rad

    # In this case, we need to convert J2000 back to PEF and then convert to
    # TOD. This is necessary because we need to apply EOP corrections to convert
    # to TOD.
    #
    # TODO: Can I simplify this rotation?
    r_MOD_J2000 = r_gcrf_to_mod_fk5(T, JD_TT)
    r_PEF_MOD   = r_mod_to_pef_fk5(T, JD_UT1, JD_TT, 0, 0)
    r_MOD_PEF   = r_pef_to_mod_fk5(T, JD_UT1, JD_TT, δΔϵ_1980, δΔψ_1980)
    r_TOD_MOD   = r_mod_to_tod_fk5(T, JD_TT, δΔϵ_1980, δΔψ_1980)

    # Compose the full rotation.
    return compose_rotation(r_MOD_J2000, r_PEF_MOD, r_MOD_PEF, r_TOD_MOD)
end

function r_eci_to_eci(T::T_ROT, ::Val{:J2000}, ::Val{:TOD}, JD_UTC::Number)
    # Get the time in TT.
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Compute and return the composed rotation.
    r_MOD_J2000 = r_gcrf_to_mod_fk5(T, JD_TT)
    r_TOD_MOD   = r_mod_to_tod_fk5(T, JD_TT, 0, 0)

    return compose_rotation(r_MOD_J2000, r_TOD_MOD)
end

@inline function r_eci_to_eci(
    T::T_ROT,
    T_ECIo::Val{:TOD},
    T_ECIf::Val{:J2000},
    JD_UTC::Number,
    eop_data::EOPData_IAU1980
)
    return inv_rotation(r_eci_to_eci(T, T_ECIf, T_ECIo, JD_UTC, eop_data))
end

@inline function r_eci_to_eci(
    T::T_ROT,
    T_ECIo::Val{:TOD},
    T_ECIf::Val{:J2000},
    JD_UTC::Number
)
    return inv_rotation(r_eci_to_eci(T, T_ECIf, T_ECIo, JD_UTC))
end

#                                J2000 <=> TEME
# ==============================================================================

@inline function r_eci_to_eci(
    T::T_ROT,
    T_ECIo::Val{:J2000},
    T_ECEFf::Val{:TEME},
    JD_UTC::Number,
    eop_data::EOPData_IAU1980
)
    return r_eci_to_eci(T, T_ECIo, T_ECEFf, JD_UTC)
end

@inline function r_eci_to_eci(
    T::T_ROT,
    T_ECIo::Val{:TEME},
    T_ECEFf::Val{:J2000},
    JD_UTC::Number,
    eop_data::EOPData_IAU1980
)
    return r_eci_to_eci(T, T_ECIo, T_ECEFf, JD_UTC)
end

function r_eci_to_eci(T::T_ROT, ::Val{:J2000}, ::Val{:TEME}, JD_UTC::Number)
    # Get the time in TT.
    JD_TT  = jd_utc_to_tt(JD_UTC)

    # Return the rotation.
    return r_gcrf_to_teme(T, JD_TT, 0, 0)
end

@inline function r_eci_to_eci(
    T::T_ROT,
    T_ECIo::Val{:TEME},
    T_ECIf::Val{:J2000},
    JD_UTC::Number
)
    return inv_rotation(r_eci_to_eci(T, T_ECIf, T_ECIo, JD_UTC))
end

#                          Between MOD, TOD, and TEME
# ==============================================================================

function r_eci_to_eci(
    T::T_ROT,
    T_ECIo::T_ECIs_of_date,
    JD_UTCo::Number,
    T_ECIf::T_ECIs_of_date,
    JD_UTCf::Number,
    eop_data::EOPData_IAU1980
)
    # In this case, we convert origin to GCRF and then convert back to the
    # destination. This is necessary because the user may want to change the
    # epoch.
    r_GCRF_ECIo = r_eci_to_eci(T, T_ECIo,     Val(:GCRF), JD_UTCo, eop_data)
    r_ECIf_GCRF = r_eci_to_eci(T, Val(:GCRF), T_ECIf,     JD_UTCf, eop_data)

    # Return the full rotation.
    return compose_rotation(r_GCRF_ECIo, r_ECIf_GCRF)
end

function r_eci_to_eci(
    T::T_ROT,
    T_ECIo::T_ECIs_of_date,
    JD_UTCo::Number,
    T_ECIf::T_ECIs_of_date,
    JD_UTCf::Number
)
    # In this case, in which we do not have EOP data, we convert origin to J2000
    # and then convert back to the destination. This is necessary because the
    # user may want to change the epoch.
    r_GCRF_ECIo = r_eci_to_eci(T, T_ECIo,      Val(:J2000), JD_UTCo)
    r_ECIf_GCRF = r_eci_to_eci(T, Val(:J2000), T_ECIf,      JD_UTCf)

    # Return the full rotation.
    return compose_rotation(r_GCRF_ECIo, r_ECIf_GCRF)
end

################################################################################
#                           IAU-2006/2010 CIO-based
################################################################################

#                                GCRF <=> CIRS
# ==============================================================================

function r_eci_to_eci(
    T::T_ROT,
    ::Val{:GCRF},
    ::Val{:CIRS},
    JD_UTC::Number,
    eop_data::EOPData_IAU2000A
)
    milliarcsec_to_rad = π / 648000000

    # Get the time in TT.
    JD_TT = jd_utc_to_tt(JD_UTC)

    # Get the EOP data related to the desired epoch.
    dX = eop_data.dX(JD_UTC) * milliarcsec_to_rad
    dY = eop_data.dY(JD_UTC) * milliarcsec_to_rad

    # Compute and return the rotation.
    return r_gcrf_to_cirs_iau2006(T, JD_TT, dX, dY)
end

@inline function r_eci_to_eci(
    T::T_ROT,
    T_ECIo::Val{:CIRS},
    T_ECIf::Val{:GCRF},
    JD_UTC::Number,
    eop_data::EOPData_IAU2000A
)
    return inv_rotation(r_eci_to_eci(T, T_ECIf, T_ECIo, JD_UTC, eop_data))
end

function r_eci_to_eci(T::T_ROT, ::Val{:GCRF}, ::Val{:CIRS}, JD_UTC::Number)
    # Get the time in TT.
    JD_TT = jd_utc_to_tt(JD_UTC)

    # Compute and return the rotation.
    return r_gcrf_to_cirs_iau2006(T, JD_TT)
end

@inline function r_eci_to_eci(
    T::T_ROT,
    T_ECIo::Val{:CIRS},
    T_ECIf::Val{:GCRF},
    JD_UTC::Number
)
    return inv_rotation(r_eci_to_eci(T, T_ECIf, T_ECIo, JD_UTC))
end

#                                 Between CIRS
# ==============================================================================

function r_eci_to_eci(
    T::T_ROT,
    T_ECIo::Val{:CIRS},
    JD_UTCo::Number,
    T_ECIf::Val{:CIRS},
    JD_UTCf::Number,
    eop_data::EOPData_IAU2000A
)
    # In this case, we convert origin to GCRF and then convert back to the
    # destination. This is necessary because the user may want to change the
    # epoch.
    r_GCRF_ECIo = r_eci_to_eci(T, T_ECIo,     Val(:GCRF), JD_UTCo, eop_data)
    r_ECIf_GCRF = r_eci_to_eci(T, Val(:GCRF), T_ECIf,     JD_UTCf, eop_data)

    # Return the full rotation.
    return compose_rotation(r_GCRF_ECIo, r_ECIf_GCRF)
end

function r_eci_to_eci(
    T::T_ROT,
    T_ECIo::Val{:CIRS},
    JD_UTCo::Number,
    T_ECIf::Val{:CIRS},
    JD_UTCf::Number
)
    # In this case, we convert origin to GCRF and then convert back to the
    # destination. This is necessary because the user may want to change the
    # epoch.
    r_GCRF_ECIo = r_eci_to_eci(T, T_ECIo,     Val(:GCRF), JD_UTCo)
    r_ECIf_GCRF = r_eci_to_eci(T, Val(:GCRF), T_ECIf,     JD_UTCf)

    # Return the full rotation.
    return compose_rotation(r_GCRF_ECIo, r_ECIf_GCRF)
end

################################################################################
#                         IAU-2006/2010 equinox-based
################################################################################

#                               GCRF <=> MJ2000
# ==============================================================================

function r_eci_to_eci(
    T::T_ROT,
    ::Val{:MJ2000},
    ::Val{:GCRF},
    JD_UTC::Number,
    eop_data::EOPData_IAU2000A
)
    return r_mj2000_to_gcrf_iau2006(T)
end

function r_eci_to_eci(
    T::T_ROT,
    ::Val{:MJ2000},
    ::Val{:GCRF},
    JD_UTC::Number
)
    return r_mj2000_to_gcrf_iau2006(T)
end

function r_eci_to_eci(
    T::T_ROT,
    ::Val{:GCRF},
    ::Val{:MJ2000},
    JD_UTC::Number,
    eop_data::EOPData_IAU2000A
)
    return r_gcrf_to_mj2000_iau2006(T)
end

function r_eci_to_eci(T::T_ROT, ::Val{:GCRF}, ::Val{:MJ2000}, JD_UTC::Number)
    return r_gcrf_to_mj2000_iau2006(T)
end

#                                 GCRF <=> MOD
# ==============================================================================

function r_eci_to_eci(
    T::T_ROT,
    ::Val{:MOD06},
    ::Val{:GCRF},
    JD_UTC::Number,
    eop_data::EOPData_IAU2000A
)
    return r_eci_to_eci(T, Val(:MOD06), Val(:GCRF), JD_UTC)
end

function r_eci_to_eci(T::T_ROT, ::Val{:MOD06}, ::Val{:GCRF}, JD_UTC::Number)
    # Get the time in TT.
    JD_TT = jd_utc_to_tt(JD_UTC)

    # Compute and return the composed rotation.
    r_MJ2000_MOD  = r_mod_to_mj2000_iau2006(T, JD_TT)
    r_GCRF_MJ2000 = r_mj2000_to_gcrf_iau2006(T)

    return compose_rotation(r_MJ2000_MOD, r_GCRF_MJ2000)
end

function r_eci_to_eci(
    T::T_ROT,
    ::Val{:GCRF},
    ::Val{:MOD06},
    JD_UTC::Number,
    eop_data::EOPData_IAU2000A
)
    return r_eci_to_eci(T, Val(:GCRF), Val(:MOD06), JD_UTC)
end

function r_eci_to_eci(
    T::T_ROT,
    ::Val{:GCRF},
    ::Val{:MOD06},
    JD_UTC::Number
)
    return inv_rotation(r_eci_to_eci(T, Val(:MOD06), Val(:GCRF), JD_UTC))
end

#                                 GCRF <=> ERS
# ==============================================================================

function r_eci_to_eci(
    T::T_ROT,
    ::Val{:ERS},
    ::Val{:GCRF},
    JD_UTC::Number,
    eop_data::EOPData_IAU2000A
)
    milliarcsec_to_rad = π / 648000000

    # Get the time in TT.
    JD_TT = jd_utc_to_tt(JD_UTC)

    # Obtain the correction of the nutation in obliquity and longitude.
    δΔϵ_2000, δΔΨ_2000 = deps_dpsi(eop_data, JD_UTC)
    δΔϵ_2000 *= milliarcsec_to_rad
    δΔΨ_2000 *= milliarcsec_to_rad

    # Compute and return the composed rotation.
    r_MOD_ERS     = r_ers_to_mod_iau2006(T, JD_TT, δΔϵ_2000, δΔΨ_2000)
    r_MJ2000_MOD  = r_mod_to_mj2000_iau2006(T, JD_TT)
    r_GCRF_MJ2000 = r_mj2000_to_gcrf_iau2006(T)

    return compose_rotation(r_MOD_ERS, r_MJ2000_MOD, r_GCRF_MJ2000)
end

function r_eci_to_eci(T::T_ROT, ::Val{:ERS}, ::Val{:GCRF}, JD_UTC::Number)
    # Get the time in TT.
    JD_TT = jd_utc_to_tt(JD_UTC)

    # Compute and return the composed rotation.
    r_MOD_ERS     = r_ers_to_mod_iau2006(T, JD_TT, 0, 0)
    r_MJ2000_MOD  = r_mod_to_mj2000_iau2006(T, JD_TT)
    r_GCRF_MJ2000 = r_mj2000_to_gcrf_iau2006(T)

    return compose_rotation(r_MOD_ERS, r_MJ2000_MOD, r_GCRF_MJ2000)
end

function r_eci_to_eci(
    T::T_ROT,
    ::Val{:GCRF},
    ::Val{:ERS},
    JD_UTC::Number,
    eop_data::EOPData_IAU2000A
)
    return inv_rotation(r_eci_to_eci(T, Val(:ERS), Val(:GCRF), JD_UTC, eop_data))
end

function r_eci_to_eci(
    T::T_ROT,
    ::Val{:GCRF},
    ::Val{:ERS},
    JD_UTC::Number
)
    return inv_rotation(r_eci_to_eci(T, Val(:ERS), Val(:GCRF), JD_UTC))
end

#                                MJ2000 <=> MOD
# ==============================================================================

function r_eci_to_eci(
    T::T_ROT,
    ::Val{:MOD06},
    ::Val{:MJ2000},
    JD_UTC::Number,
    eop_data::EOPData_IAU2000A
)
    return r_eci_to_eci(T, Val(:MOD06), Val(:MJ2000), JD_UTC)
end

function r_eci_to_eci(T::T_ROT, ::Val{:MOD06}, ::Val{:MJ2000}, JD_UTC::Number)
    # Get the time in TT.
    JD_TT = jd_utc_to_tt(JD_UTC)

    # Compute the rotation.
    return r_mod_to_mj2000_iau2006(T, JD_TT)
end

function r_eci_to_eci(
    T::T_ROT,
    ::Val{:MJ2000},
    ::Val{:MOD06},
    JD_UTC::Number,
    eop_data::EOPData_IAU2000A
)
    return r_eci_to_eci(T, Val(:MJ2000), Val(:MOD06), JD_UTC)
end

function r_eci_to_eci(
    T::T_ROT,
    ::Val{:MJ2000},
    ::Val{:MOD06},
    JD_UTC::Number
)
    return inv_rotation(r_eci_to_eci(T, Val(:MOD06), Val(:MJ2000), JD_UTC))
end

#                                 MJ2000 <=> ERS
# ==============================================================================

function r_eci_to_eci(
    T::T_ROT,
    ::Val{:ERS},
    ::Val{:MJ2000},
    JD_UTC::Number,
    eop_data::EOPData_IAU2000A
)
    milliarcsec_to_rad = π / 648000000

    # Get the time in TT.
    JD_TT = jd_utc_to_tt(JD_UTC)

    # Obtain the correction of the nutation in obliquity and longitude.
    δΔϵ_2000, δΔΨ_2000 = deps_dpsi(eop_data, JD_UTC)
    δΔϵ_2000 *= milliarcsec_to_rad
    δΔΨ_2000 *= milliarcsec_to_rad

    # Compute and return the composed rotation.
    r_MOD_ERS    = r_ers_to_mod_iau2006(T, JD_TT, δΔϵ_2000, δΔΨ_2000)
    r_MJ2000_MOD = r_mod_to_mj2000_iau2006(T, JD_TT)

    return compose_rotation(r_MOD_ERS, r_MJ2000_MOD)
end

function r_eci_to_eci(T::T_ROT, ::Val{:ERS}, ::Val{:MJ2000}, JD_UTC::Number)
    # Get the time in TT.
    JD_TT = jd_utc_to_tt(JD_UTC)

    # Compute and return the composed rotation.
    r_MOD_ERS    = r_ers_to_mod_iau2006(T, JD_TT, 0, 0)
    r_MJ2000_MOD = r_mod_to_mj2000_iau2006(T, JD_TT)

    return compose_rotation(r_MOD_ERS, r_MJ2000_MOD)
end

function r_eci_to_eci(
    T::T_ROT,
    ::Val{:MJ2000},
    ::Val{:ERS},
    JD_UTC::Number,
    eop_data::EOPData_IAU2000A
)
    return inv_rotation(r_eci_to_eci(T, Val(:ERS), Val(:MJ2000), JD_UTC, eop_data))
end

function r_eci_to_eci(
    T::T_ROT,
    ::Val{:MJ2000},
    ::Val{:ERS},
    JD_UTC::Number
)
    return inv_rotation(r_eci_to_eci(T, Val(:ERS), Val(:MJ2000), JD_UTC))
end

#                             Between ERS and MOD
# ==============================================================================

function r_eci_to_eci(
    T::T_ROT,
    T_ECIo::T_ECIs_IAU_2006_Equinox_of_date,
    JD_UTCo::Number,
    T_ECIf::T_ECIs_IAU_2006_Equinox_of_date,
    JD_UTCf::Number,
    eop_data::EOPData_IAU2000A
)
    # In this case, we convert origin to GCRF and then convert back to the
    # destination. This is necessary because the user may want to change the
    # epoch.
    r_GCRF_ECIo = r_eci_to_eci(T, T_ECIo,     Val(:GCRF), JD_UTCo, eop_data)
    r_ECIf_GCRF = r_eci_to_eci(T, Val(:GCRF), T_ECIf,     JD_UTCf, eop_data)

    # Return the full rotation.
    return compose_rotation(r_GCRF_ECIo, r_ECIf_GCRF)
end

function r_eci_to_eci(
    T::T_ROT,
    T_ECIo::T_ECIs_IAU_2006_Equinox_of_date,
    JD_UTCo::Number,
    T_ECIf::T_ECIs_IAU_2006_Equinox_of_date,
    JD_UTCf::Number
)
    # In this case, we convert origin to GCRF and then convert back to the
    # destination. This is necessary because the user may want to change the
    # epoch. Notice that, differently from IAU-76/FK5, we can convert to GCRF
    # without using EOP data with minor degradation in precsion.
    r_GCRF_ECIo = r_eci_to_eci(T, T_ECIo,      Val(:GCRF), JD_UTCo)
    r_ECIf_GCRF = r_eci_to_eci(T, Val(:GCRF), T_ECIf,      JD_UTCf)

    # Return the full rotation.
    return compose_rotation(r_GCRF_ECIo, r_ECIf_GCRF)
end


text:
##
cd(@__DIR__())
using Underscores
using Base.Iterators: product
plen = 25
codes = @_ map(parse(Int, _), eachline("input.txt"));

##
i, invalid = Iterators.filter(enumerate(codes[plen+1:end])) do (i, c)
    foldl(product(codes[i:i+plen-1], codes[i:i+plen-1]); init=true) do p, (x, y)
        x >= y && return p
        p && x+y!=c
    end
end |> collect |> first

## PARTE 2
i = findfirst(==(invalid), codes);

##
filter(!isnothing, map(2:i-1) do l
    j = findfirst(1:i-l) do j
        sum(codes[j:j+l-1]) == invalid
    end
    isnothing(j) && return nothing
    extrema(codes[j:j+l-1]) |> sum
end)


text:
"""
Abstract type for dimension metadata wrappers.
"""
abstract type Metadata{K,V} <: AbstractDict{K,V} end

(::Type{T})() where T <: Metadata = T(Dict())

val(metadata::Metadata) = metadata.val

Base.get(metadata::Metadata, args...) = get(val(metadata), args...)
Base.getindex(metadata::Metadata, args...) = getindex(val(metadata), args...)
Base.setindex!(metadata::Metadata, args...) = setindex!(val(metadata), args...)
Base.keys(metadata::Metadata) = keys(val(metadata))
Base.iterate(metadata::Metadata, args...) = iterate(val(metadata), args...)
Base.IteratorSize(::Metadata{M})	where M = IteratorSize(M)
Base.IteratorEltype(::Metadata{M}) where M = IteratorEltype(M)
Base.eltype(::Metadata{M}) where M = eltype(M)
Base.length(metadata::Metadata) = length(val(metadata))

abstract type AbstractDimMetadata{K,V} <: Metadata{K,V} end

abstract type AbstractArrayMetadata{K,V} <: Metadata{K,V} end

struct ArrayMetadata{K,V} <: AbstractArrayMetadata{K,V}
    val::Dict{K,V}
end

struct DimMetadata{K,V} <: AbstractDimMetadata{K,V}
    val::Dict{K,V}
end


text:
using Embeddings
countries = ["Afghanistan", "Algeria", "Angola", "Arabia", "Argentina", "Australia", "Bangladesh", "Brazil", "Britain", "Canada", "China", "Colombia", "Congo", "Egypt", "England", "Ethiopia", "France", "Germany", "Ghana", "India", "Indonesia", "Iran", "Iraq", "Ireland", "Italy", "Japan", "Kenya", "Korea", "Madagascar", "Malaysia", "Mexico", "Morocco", "Mozambique", "Myanmar", "Nepal", "Nigeria", "Pakistan", "Peru", "Philippines", "Poland", "Russia", "South", "Spain", "Sudan", "Tanzania", "Thailand", "Uganda", "Ukraine", "Usa", "Uzbekistan", "Venezuela", "Vietnam", "Wales", "Yemen"]
usa_cities = ["Albuquerque", "Atlanta", "Austin", "Baltimore", "Boston", "Charlotte", "Chicago", "Columbus", "Dallas", "Denver", "Detroit", "Francisco", "Fresno", "Houston", "Indianapolis", "Jacksonville", "Las", "Louisville", "Memphis", "Mesa", "Milwaukee", "Nashville", "Omaha", "Philadelphia", "Phoenix", "Portland", "Raleigh", "Sacramento", "San", "Seattle", "Tucson", "Vegas", "Washington"]
world_capitals = ["Accra", "Algiers", "Amman", "Ankara", "Antananarivo", "Athens", "Baghdad", "Baku", "Bangkok", "Beijing", "Beirut", "Berlin", "Bogotá", "Brasília", "Bucharest", "Budapest", "Cairo", "Caracas", "Damascus", "Dhaka", "Hanoi", "Havana", "Jakarta", "Kabul", "Kampala", "Khartoum", "Kinshasa", "Kyiv", "Lima", "London", "Luanda", "Madrid", "Manila", "Minsk", "Moscow", "Nairobi", "Paris", "Pretoria", "Pyongyang", "Quito", "Rabat", "Riyadh", "Rome", "Santiago", "Seoul", "Singapore", "Stockholm", "Taipei", "Tashkent", "Tehran", "Tokyo", "Vienna", "Warsaw", "Yaoundé"]
animals = ["alpaca","camel","cattle","dog","dove","duck","ferret","goldfish","goose","rat","llama","mouse","pigeon","yak"]
china_cities = ["Shanghai", "Tianjin", "Shenzhen", "Chengdu", "Changsha", "Nanjing", "Chongqing", "Hefei", "Qinghai", "Jilin", "Heilongjiang"]
sports = ["archery","badminton","basketball","boxing","cycling","diving","equestrian","fencing","field","football","golf","gymnastics","handball","hockey","judo","kayak","pentathlon","polo","rowing","rugby","sailing","shooting","soccer","swimming","taekwondo","tennis","triathlon","volleyball","weightlifting","wrestling"]

words_by_class = [countries, usa_cities, world_capitals, china_cities, animals, sports]
all_words = reduce(vcat, words_by_class)
embedding_table = load_embeddings(Word2Vec; keep_words = all_words)
@assert Set(all_words) == Set(embedding_table.vocab)

embeddings = embedding_table.embeddings
all_words = embedding_table.vocab
classes = map(all_words) do word
    findfirst(col -> word ∈ col, [countries, usa_cities, world_capitals, china_cities, animals, sports])
end;

##
using MultivariateStats
using Plots
plotly()

#Direct projection -- no DR -- just throw away the information in the other axies
xs=embeddings
scatter(xs[1,:], xs[2,:], xs[3,:]; hover=all_words, zcolor=classes, color = :viridis, legend = false)

##
M = fit(PCA, embeddings; maxoutdim=3)
xs = transform(M, embeddings)
scatter(xs[1,:], xs[2,:], xs[3,:]; hover=all_words, zcolor=classes, color = :viridis, legend = false)

M = fit(ICA, embeddings, 3; maxiter = 1000, tol = 1)
xs = transform(M, embeddings)
scatter(xs[1,:], xs[2,:], xs[3,:]; hover=all_words, zcolor=classes, color = :viridis, legend = false)

using TSne
xs = tsne(embeddings', 3, 3, 1000, 20.0)'
scatter(xs[1,:], xs[2,:], xs[3,:]; hover=all_words, zcolor=classes, color = :viridis, legend = false)


text:
@testset "MPS based search finds the correct low energy spectrum" begin
    instance = "$(@__DIR__)/instances/pathological/cross_3_4_dd.txt"

    ig = ising_graph(instance)

    expected_energies = [-23.301855, -23.221513, -23.002799, -22.922457]

    max_states = 200
    to_show = length(expected_energies)

    β = 2.
    dβ = β/8.0

    Dcut = 32
    var_ϵ = 1E-8
    max_sweeps = 4

    @testset "without any purification" begin
        igp = prune(ig)
        schedule = fill(dβ, Int(ceil(β/dβ)))
        ψ = MPS(igp, Dcut, var_ϵ, max_sweeps, schedule)
        states, lprob, _ = solve(ψ, max_states)
        @test sort(energy.(states[1:to_show], Ref(igp))) ≈ expected_energies
    end
end