text:
import Base: convert

export InputConfig, InputTol, SolverParams, PreallocatedData

# problem: min 1/2 x'Qx + c'x + c0     s.t.  Ax = b,  lvar ≤ x ≤ uvar
abstract type Abstract_QM_FloatData{T<:Real} end

mutable struct QM_FloatData{T<:Real} <: Abstract_QM_FloatData{T}
    Q     :: SparseMatrixCSC{T,Int} # size nvar * nvar
    AT    :: SparseMatrixCSC{T,Int} # size ncon * nvar, using Aᵀ is easier to form systems
    b     :: Vector{T} # size ncon
    c     :: Vector{T} # size nvar
    c0    :: T
    lvar  :: Vector{T} # size nvar
    uvar  :: Vector{T} # size nvar
end

mutable struct QM_IntData
    ilow   :: Vector{Int} # indices of finite elements in lvar
    iupp   :: Vector{Int} # indices of finite elements in uvar
    irng   :: Vector{Int} # indices of finite elements in both lvar and uvar
    ifree  :: Vector{Int} # indices of infinite elements in both lvar and uvar
    ncon   :: Int # number of equality constraints after SlackModel! (= size of b)
    nvar   :: Int # number of variables
    nlow   :: Int # length(ilow)
    nupp   :: Int # length(iupp)
end

"""
Abstract type for tuning the parameters of the different solvers.
Each solver has its own `SolverParams` type.

The `SolverParams` currently implemented within RipQP are:

- [`RipQP.K2LDLParams`](@ref)
- [`RipQP.K2_5LDLParams`](@ref)

"""
abstract type SolverParams end

"""
Type to specify the configuration used by RipQP.

- `mode :: Symbol`: should be `:mono` to use the mono-precision mode, or `:multi` to use
    the multi-precision mode (start in single precision and gradually transitions
    to `T0`)
- `scaling :: Bool`: activate/deactivate scaling of A and Q in `QM0`
- `normalize_rtol :: Bool = true` : if `true`, the primal and dual tolerance for the stopping criteria
    are normalized by the initial primal and dual residuals
- `kc :: Int`: number of centrality corrections (set to `-1` for automatic computation)
- `refinement :: Symbol` : should be `:zoom` to use the zoom procedure, `:multizoom` to use the zoom procedure
    with multi-precision (then `mode` should be `:multi`), `ref` to use the QP refinement procedure, `multiref`
    to use the QP refinement procedure with multi_precision (then `mode` should be `:multi`), or `none` to avoid
    refinements
- `sp :: SolverParams` : choose a solver to solve linear systems that occurs at each iteration and during the
    initialization, see [`RipQP.SolverParams`](@ref)
- `solve_method :: Symbol` : method used to solve the system at each iteration, use `solve_method = :PC` to
    use the Predictor-Corrector algorithm (default), and use `solve_method = :IPF` to use the Infeasible Path
    Following algorithm

The constructor

    iconf = InputConfig(; mode :: Symbol = :mono, scaling :: Bool = true,
                        normalize_rtol :: Bool = true, kc :: I = 0,
                        refinement :: Symbol = :none, max_ref :: I = 0,
                        sp :: SolverParams = K2LDLParams(),
                        solve_method :: Symbol = :PC) where {I<:Integer}

returns a `InputConfig` struct that shall be used to solve the input `QuadraticModel` with RipQP.
"""
struct InputConfig{I<:Integer}
    mode                :: Symbol
    scaling             :: Bool
    normalize_rtol      :: Bool # normalize the primal and dual tolerance to the initial starting primal and dual residuals
    kc                  :: I # multiple centrality corrections, -1 = automatic computation

    # QP refinement
    refinement          :: Symbol
    max_ref             :: I # maximum number of refinements

    # Functions to choose formulations
    sp                  :: SolverParams
    solve_method        :: Symbol
end

function InputConfig(; mode :: Symbol = :mono, scaling :: Bool = true, normalize_rtol :: Bool = true,
                      kc :: I = 0, refinement :: Symbol = :none, max_ref :: I = 0, sp :: SolverParams = K2LDLParams(),
                      solve_method :: Symbol = :PC) where {I<:Integer}

    mode == :mono || mode == :multi || error("mode should be :mono or :multi")
    refinement == :zoom || refinement == :multizoom || refinement == :ref || refinement == :multiref ||
        refinement == :none || error("not a valid refinement parameter")
    solve_method == :IPF && kc != 0 && error("IPF method should not be used with centrality corrections")

    return InputConfig{I}(mode, scaling, normalize_rtol, kc, refinement, max_ref, sp, solve_method)
end

"""
Type to specify the tolerances used by RipQP.

- `max_iter :: Int`: maximum number of iterations
- `ϵ_pdd`: relative primal-dual difference tolerance
- `ϵ_rb`: primal tolerance
- `ϵ_rc`: dual tolerance
- `max_iter32`, `ϵ_pdd32`, `ϵ_rb32`, `ϵ_rc32`: same as `max_iter`, `ϵ_pdd`, `ϵ_rb` and
    `ϵ_rc`, but used for switching from single precision to double precision. They are
    only usefull when `mode=:multi`
- `max_iter64`, `ϵ_pdd64`, `ϵ_rb64`, `ϵ_rc64`: same as `max_iter`, `ϵ_pdd`, `ϵ_rb` and
    `ϵ_rc`, but used for switching from double precision to quadruple precision. They
    are only usefull when `mode=:multi` and `T0=Float128`
- `ϵ_rbz` : primal transition tolerance for the zoom procedure, (used only if `refinement=:zoom`)
- `ϵ_Δx`: step tolerance for the current point estimate (note: this criterion
    is currently disabled)
- `ϵ_μ`: duality measure tolerance (note: this criterion is currently disabled)
- `max_time`: maximum time to solve the QP

The constructor

    itol = InputTol(;max_iter :: I = 200, max_iter32 :: I = 40, max_iter64 :: I = 180,
                     ϵ_pdd :: T = 1e-8, ϵ_pdd32 :: T = 1e-2, ϵ_pdd64 :: T = 1e-4,
                     ϵ_rb :: T = 1e-6, ϵ_rb32 :: T = 1e-4, ϵ_rb64 :: T = 1e-5, ϵ_rbz :: T = 1e-3,
                     ϵ_rc :: T = 1e-6, ϵ_rc32 :: T = 1e-4, ϵ_rc64 :: T = 1e-5,
                     ϵ_Δx :: T = 1e-16, ϵ_μ :: T = 1e-9) where {T<:Real, I<:Integer}

returns a `InputTol` struct that initializes the stopping criteria for RipQP.
The 32 and 64 characters refer to the stopping criteria in `:multi` mode for the transitions from `Float32` to `Float64`
and `Float64` to `Float128` (if the input `QuadraticModel` is in `Float128`) respectively.
"""
struct InputTol{T<:Real, I<:Integer}
    # maximum number of iterations
    max_iter        :: I
    max_iter32      :: I # only in multi mode
    max_iter64      :: I # only in multi mode with T0 = Float128

    # relative primal-dual gap tolerance
    ϵ_pdd           :: T
    ϵ_pdd32         :: T # only in multi mode
    ϵ_pdd64         :: T # only in multi mode with T0 = Float128

    # primal residual
    ϵ_rb            :: T
    ϵ_rb32          :: T # only in multi mode
    ϵ_rb64          :: T # only in multi mode with T0 = Float128
    ϵ_rbz           :: T # only when using zoom refinement

    # dual residual
    ϵ_rc            :: T
    ϵ_rc32          :: T # only in multi mode
    ϵ_rc64          :: T # only in multi mode with T0 = Float128

    # unused residuals (for now)
    ϵ_μ             :: T
    ϵ_Δx            :: T

    # maximum time for resolution
    max_time        :: T
end

function InputTol(;max_iter :: I = 200, max_iter32 :: I = 40, max_iter64 :: I = 180,
                   ϵ_pdd :: T = 1e-8, ϵ_pdd32 :: T = 1e-2, ϵ_pdd64 :: T = 1e-4,
                   ϵ_rb :: T = 1e-6, ϵ_rb32 :: T = 1e-4, ϵ_rb64 :: T = 1e-5, ϵ_rbz :: T = 1e-5,
                   ϵ_rc :: T = 1e-6, ϵ_rc32 :: T = 1e-4, ϵ_rc64 :: T = 1e-5,
                   ϵ_Δx :: T = 1e-16, ϵ_μ :: T = 1e-9, max_time :: T = 1200.) where {T<:Real, I<:Integer}

    return InputTol{T, I}(max_iter, max_iter32, max_iter64, ϵ_pdd, ϵ_pdd32, ϵ_pdd64, ϵ_rb, ϵ_rb32, ϵ_rb64, ϵ_rbz,
                           ϵ_rc, ϵ_rc32, ϵ_rc64, ϵ_μ, ϵ_Δx, max_time)
end

mutable struct Tolerances{T<:Real}
    pdd              :: T  # primal-dual difference (relative)
    rb               :: T  # primal residuals tolerance
    rc               :: T  # dual residuals tolerance
    tol_rb           :: T  # ϵ_rb * (1 + ||r_b0||)
    tol_rc           :: T  # ϵ_rc * (1 + ||r_c0||)
    μ                :: T  # duality measure
    Δx               :: T
    normalize_rtol   :: Bool # true if normalize_rtol=true, then tol_rb, tol_rc = ϵ_rb, ϵ_rc
end

mutable struct Point{T<:Real}
    x    :: Vector{T} # size nvar
    y    :: Vector{T} # size ncon
    s_l  :: Vector{T} # size nlow (useless zeros corresponding to infinite lower bounds are not stored)
    s_u  :: Vector{T} # size nupp (useless zeros corresponding to infinite upper bounds are not stored)
end

convert(::Type{Point{T}}, pt) where {T<:Real} = Point(convert(Array{T}, pt.x), convert(Array{T}, pt.y),
                                                      convert(Array{T}, pt.s_l), convert(Array{T}, pt.s_u))

mutable struct Residuals{T<:Real}
    rb      :: Vector{T} # primal residuals Ax - b
    rc      :: Vector{T} # dual residuals -Qx + Aᵀy + s_l - s_u
    rbNorm  :: T # ||rb||
    rcNorm  :: T # ||rc||
    n_Δx    :: T # ||Δx||
end

convert(::Type{Residuals{T}}, res) where {T<:Real} = Residuals(convert(Array{T}, res.rb), convert(Array{T}, res.rc),
                                                               convert(T, res.rbNorm), convert(T, res.rcNorm),
                                                               convert(T, res.n_Δx))

# LDLFactorization conversion function
convertldl(T :: DataType, K_fact) = LDLFactorizations.LDLFactorization(K_fact.__analyzed, K_fact.__factorized, K_fact.__upper,
                                                                       K_fact.n, K_fact.parent, K_fact.Lnz, K_fact.flag,
                                                                       K_fact.P, K_fact.pinv, K_fact.Lp, K_fact.Cp,
                                                                       K_fact.Ci, K_fact.Li, convert(Array{T}, K_fact.Lx),
                                                                       convert(Array{T}, K_fact.d), convert(Array{T}, K_fact.Y),
                                                                       K_fact.pattern, T(K_fact.r1), T(K_fact.r2),
                                                                       T(K_fact.tol), K_fact.n_d)

mutable struct Regularization{T<:Real}
    ρ        :: T       # curent top-left regularization parameter
    δ        :: T       # cureent bottom-right regularization parameter
    ρ_min    :: T       # ρ minimum value
    δ_min    :: T       # δ minimum value
    regul    :: Symbol  # regularization mode (:classic, :dynamic, or :none)
end

convert(::Type{Regularization{T}}, regu::Regularization{T0}) where {T<:Real, T0<:Real} =
    Regularization(T(regu.ρ), T(regu.δ), T(regu.ρ_min), T(regu.δ_min), regu.regul)

mutable struct IterData{T<:Real}
    Δxy         :: Vector{T} # Newton step [Δx; Δy]
    Δs_l        :: Vector{T}
    Δs_u        :: Vector{T}
    x_m_lvar    :: Vector{T} # x - lvar
    uvar_m_x    :: Vector{T} # uvar - x
    Qx          :: Vector{T}
    ATy         :: Vector{T} # Aᵀy
    Ax          :: Vector{T}
    xTQx_2      :: T # xᵀQx
    cTx         :: T # cᵀx
    pri_obj     :: T # 1/2 xᵀQx + cᵀx + c0
    dual_obj    :: T # -1/2 xᵀQx + yᵀb + s_lᵀlvar - s_uᵀuvar + c0
    μ           :: T # duality measure (s_lᵀ(x-lvar) + s_uᵀ(uvar-x)) / (nlow+nupp)
    pdd         :: T # primal dual difference (relative) pri_obj - dual_obj / pri_obj
    l_pdd       :: Vector{T} # list of the 5 last pdd
    mean_pdd    :: T # mean of the 5 last pdd
    qp          :: Bool # true if qp false if lp
end

convert(::Type{IterData{T}}, itd :: IterData{T0}) where {T<:Real, T0<:Real} =
    IterData(convert(Array{T}, itd.Δxy),
             convert(Array{T}, itd.Δs_l),
             convert(Array{T}, itd.Δs_u),
             convert(Array{T}, itd.x_m_lvar),
             convert(Array{T}, itd.uvar_m_x),
             convert(Array{T}, itd.Qx),
             convert(Array{T}, itd.ATy),
             convert(Array{T}, itd.Ax),
             convert(T, itd.xTQx_2),
             convert(T, itd.cTx),
             convert(T, itd.pri_obj),
             convert(T, itd.dual_obj),
             convert(T, itd.μ),
             convert(T, itd.pdd),
             convert(Array{T}, itd.l_pdd),
             convert(T, itd.mean_pdd),
             itd.qp
             )

abstract type PreallocatedData{T<:Real} end

mutable struct StopCrit{T}
    optimal     :: Bool
    small_Δx    :: Bool
    small_μ     :: Bool
    tired       :: Bool
    max_iter    :: Int
    max_time    :: T
    start_time  :: T
    Δt          :: T
end

mutable struct Counters
    c_catch  :: Int # safety try:cath
    c_pdd    :: Int # maximum number of δ_min reductions when pdd does not change
    k        :: Int # iter count
    km       :: Int # iter relative to precision: if k+=1 and T==Float128, km +=16  (km+=4 if T==Float64 and km+=1 if T==Float32)
    kc       :: Int # maximum corrector steps
    max_ref  :: Int # maximum number of refinements
    c_ref    :: Int # current number of refinements
end


text:
# :LN_BOBYQA
@testset "f_3p_1im restricted :LN_BOBYQA" begin
    res0 = [get_right_endpoint(
        [3., 8., 2.1],
        i,
        (x::Vector{Float64}) -> f_3p_1im(x) - 9.,
        Val(:CICO_ONE_PASS);
        scan_bound = 4.,
        local_alg = :LN_BOBYQA
    ) for i in 1:3]
    @test res0[1][3] == :SCAN_BOUND_REACHED
    @test res0[2][3] == :SCAN_BOUND_REACHED
    @test res0[3][3] == :SCAN_BOUND_REACHED
end

@testset "f_3p_1im_dep :LN_BOBYQA" begin
    res0 = [get_right_endpoint(
        [3., 2., 2.1],
        i,
        (x::Vector{Float64}) -> f_3p_1im_dep(x) - 9.,
        Val(:CICO_ONE_PASS);
        local_alg = :LN_BOBYQA,
        loss_tol = 1e-3
    ) for i in 1:3]
    @test isapprox(res0[1][1], 5., atol=1e-2) # XXX: bad tolerance
    @test res0[1][3] == :BORDER_FOUND_BY_SCAN_TOL
    @test isapprox(res0[2][1], 2.0+2.0*sqrt(2.), atol=1e-2) # XXX: bad tolerance
    @test res0[2][3] == :BORDER_FOUND_BY_SCAN_TOL
    @test res0[3][3] == :SCAN_BOUND_REACHED
end

@testset "f_3p_1im_dep scan_bound :LN_BOBYQA" begin
    res0 = [get_right_endpoint(
        [3., 2., 2.1],
        i,
        (x::Vector{Float64}) -> f_3p_1im_dep(x) - 9.,
        Val(:CICO_ONE_PASS);
        scan_bound = [4.,10.,10.][i],
        local_alg = :LN_BOBYQA
    ) for i in 1:3]
    @test res0[1][3] == :SCAN_BOUND_REACHED
    @test isapprox(res0[2][1], 2.0+2.0*sqrt(2.), atol=1e-2) # XXX: bad tolerance
    @test res0[2][3] == :BORDER_FOUND_BY_SCAN_TOL
    @test res0[3][3] == :SCAN_BOUND_REACHED
end

# :LN_NEWUOA
@testset "f_3p_1im restricted :LN_NEWUOA" begin
    res0 = [get_right_endpoint(
        [3., 8., 2.1],
        i,
        (x::Vector{Float64}) -> f_3p_1im(x) - 9.,
        Val(:CICO_ONE_PASS);
        scan_bound = 4.,
        local_alg = :LN_NEWUOA
    ) for i in 1:3]
    @test res0[1][3] == :SCAN_BOUND_REACHED # XXX: fail
    @test res0[2][3] == :SCAN_BOUND_REACHED
    @test res0[3][3] == :SCAN_BOUND_REACHED # XXX: fail
end

@testset "f_3p_1im_dep :LN_NEWUOA" begin
    res0 = [get_right_endpoint(
        [3., 2., 2.1],
        i,
        (x::Vector{Float64}) -> f_3p_1im_dep(x) - 9.,
        Val(:CICO_ONE_PASS);
        local_alg = :LN_NEWUOA
    ) for i in 1:3]
    @test isapprox(res0[1][1], 5., atol=1e-2)
    @test res0[1][3] == :BORDER_FOUND_BY_SCAN_TOL
    @test isapprox(res0[2][1], 2.0+2.0*sqrt(2.), atol=1e-2)
    @test res0[2][3] == :BORDER_FOUND_BY_SCAN_TOL
    @test res0[3][3] == :SCAN_BOUND_REACHED
end

@testset "f_3p_1im_dep scan_bound :LN_NEWUOA" begin
    res0 = [get_right_endpoint(
        [3., 2., 2.1],
        i,
        (x::Vector{Float64}) -> f_3p_1im_dep(x) - 9.,
        Val(:CICO_ONE_PASS);
        scan_bound = [4.,10.,10.][i],
        local_alg = :LN_NEWUOA
    ) for i in 1:3]
    @test res0[1][3] == :SCAN_BOUND_REACHED
    @test isapprox(res0[2][1], 2.0+2.0*sqrt(2.), atol=1e-2)
    @test res0[2][3] == :BORDER_FOUND_BY_SCAN_TOL
    @test res0[3][3] == :SCAN_BOUND_REACHED
end

# :LN_SBPLX
@testset "f_3p_1im restricted :LN_SBPLX" begin
    res0 = [get_right_endpoint(
        [3., 8., 2.1],
        i,
        (x::Vector{Float64}) -> f_3p_1im(x) - 9.,
        Val(:CICO_ONE_PASS);
        scan_bound = 4.,
        local_alg = :LN_SBPLX
    ) for i in 1:3]
    @test res0[1][3] == :SCAN_BOUND_REACHED
    @test res0[2][3] == :SCAN_BOUND_REACHED
    @test res0[3][3] == :SCAN_BOUND_REACHED
end

@testset "f_3p_1im_dep :LN_SBPLX" begin
    res0 = [get_right_endpoint(
        [3., 2., 2.1],
        i,
        (x::Vector{Float64}) -> f_3p_1im_dep(x) - 9.,
        Val(:CICO_ONE_PASS);
        local_alg = :LN_SBPLX
    ) for i in 1:3]
    @test isapprox(res0[1][1], 5., atol=1e-2)
    @test res0[1][3] == :BORDER_FOUND_BY_SCAN_TOL
    @test isapprox(res0[2][1], 2.0+2.0*sqrt(2.), atol=1e-2) # XXX: bad tolerance
    @test res0[2][3] == :BORDER_FOUND_BY_SCAN_TOL
    @test res0[3][3] == :SCAN_BOUND_REACHED
end

@testset "f_3p_1im_dep scan_bound :LN_SBPLX" begin
    res0 = [get_right_endpoint(
        [3., 2., 2.1],
        i,
        (x::Vector{Float64}) -> f_3p_1im_dep(x) - 9.,
        Val(:CICO_ONE_PASS);
        scan_bound = [4.,10.,10.][i],
        local_alg = :LN_SBPLX
    ) for i in 1:3]
    @test res0[1][3] == :SCAN_BOUND_REACHED
    @test isapprox(res0[2][1], 2.0+2.0*sqrt(2.), atol=1e-2) # XXX: bad tolerance
    @test res0[2][3] == :BORDER_FOUND_BY_SCAN_TOL
    @test res0[3][3] == :SCAN_BOUND_REACHED
end

ext: jl


text:
"""
dyn_frwd takes care of "forward" construction of double-encoding sequences during the HMM-parameterized step of CAMEOS.
	Orchestrates the construction of the tensor containing appropriate scores for encodings across positions in gene x and y.
	Calculates the costs of substitutions and transitions between states and is the place where recursion back to the
	start of the sequence (through code in dyn_back) is initiated. "dyn" from "Dyn"amic programming.
"""

include("math.jl")
include("dyn_back.jl")

FORCE_MET = true #whether we need all proteins to start with a M. Generally true, unless dealing with fragments.
trinucs = ["AAA", "AAC", "AAG", "AAT", "ACA", "ACC", "ACG", "ACT", "AGA", "AGC", "AGG", "AGT", "ATA", "ATC", "ATG", "ATT", "CAA", "CAC", "CAG",
				 "CAT", "CCA", "CCC", "CCG", "CCT", "CGA", "CGC", "CGG", "CGT", "CTA", "CTC", "CTG", "CTT", "GAA", "GAC", "GAG", "GAT", "GCA", "GCC",
				 "GCG", "GCT", "GGA", "GGC", "GGG", "GGT", "GTA", "GTC", "GTG", "GTT", "TAA", "TAC", "TAG", "TAT", "TCA", "TCC", "TCG", "TCT", "TGA",
				 "TGC", "TGG", "TGT", "TTA", "TTC", "TTG", "TTT"]

function recurse_back(X_hmm, Y_hmm, short_node, far_node, cur_state, cur_score, last_letter, full_length, cumulative_score, path, bbt, thresh=1.0, temperature=1.0)

	short_mm::Float64, short_md::Float64, short_im::Float64,
	short_dm::Float64, short_dd::Float64 = state_probs(X_hmm, short_node - 1)[[1, 3, 4, 6, 7]]

	far_mm::Float64, far_md::Float64, far_im::Float64,
	far_dm::Float64, far_dd::Float64 = state_probs(Y_hmm, far_node - 1)[[1, 3, 4, 6, 7]]

	short_mi::Float64, short_ii::Float64 = state_probs(X_hmm, short_node)[[2, 5]]
	far_mi::Float64, far_ii::Float64 = state_probs(Y_hmm, far_node)[[2, 5]]

	all_vals::Array{Float64, 1} = [short_mm, short_md, short_im, short_dm, short_dd, short_mi, short_ii, far_mm, far_md, far_im, far_dm, far_dd, far_mi, far_ii]

	next_state = 0
	next_codon = ""
	next_stop = 0

	# 1) MxMy, 2) IxMy, 3) DxMy, 4) MxIy, 5) DxIy, 6) MxDy, 7) IxDy
	if cur_state == 1 #Meaning we're in MxMy
		next_state, next_letter, next_codon, cur_score, cumulative_score = dyn_back.mxmy(trinucs, last_letter, short_node, far_node, cumulative_score, cur_score, all_vals, bbt, thresh, temperature)
	elseif cur_state == 2
		next_state, next_letter, next_codon, cur_score, cumulative_score = dyn_back.ixmy(trinucs, last_letter, short_node, far_node, cumulative_score, cur_score, all_vals, bbt, thresh, temperature)
	elseif cur_state == 3
		next_state, next_letter, next_codon, cur_score, cumulative_score = dyn_back.dxmy(trinucs, last_letter, short_node, far_node, cumulative_score, cur_score, all_vals, bbt, thresh, temperature)
	elseif cur_state == 4
		next_state, next_letter, next_codon, cur_score, cumulative_score = dyn_back.mxiy(trinucs, last_letter, short_node, far_node, cumulative_score, cur_score, all_vals, bbt, thresh, temperature)
	elseif cur_state == 5
		next_state, next_letter, next_codon, cur_score, cumulative_score = dyn_back.dxiy(trinucs, last_letter, short_node, far_node, cumulative_score, cur_score, all_vals, bbt, thresh, temperature)
	elseif cur_state == 6
		next_state, next_letter, next_codon, cur_score, cumulative_score = dyn_back.mxdy(trinucs, last_letter, short_node, far_node, cumulative_score, cur_score, all_vals, bbt, thresh, temperature)
	elseif cur_state == 7
		next_state, next_letter, next_codon, cur_score, cumulative_score = dyn_back.ixdy(trinucs, last_letter, short_node, far_node, cumulative_score, cur_score, all_vals, bbt, thresh, temperature)
	end

	#Depending on the state, we move zero or one aa closer to start in both frames.
	if next_state == 1 || next_state == 3 || next_state == 4 || next_state == 5 || next_state == 6
		short_node -= 1
	end
	if next_state == 1 || next_state == 2 || next_state == 3 || next_state == 6 || next_state == 7
		far_node -= 1
	end

	push!(path, (short_node, far_node, next_state))
	if short_node == 1 || far_node == 1 #then we're at a base case of recursion being done.
		return full_length, path, cumulative_score, cur_score
	else
		new_full_length = next_codon * full_length
		return recurse_back(X_hmm, Y_hmm, short_node, far_node, next_state, cur_score, string("ACGT"[next_letter]), new_full_length, cumulative_score, path, bbt, thresh, temperature)
	end
end

function state_transition_calcs(tensor::Array{Float32, 4}, short_node::Int, far_node::Int, X_hmm::types.wholeHMM, Y_hmm::types.wholeHMM, the_keys::Array{Int, 1}, the_vals::Array{Float64, 1},
                      mxmy_calcs::Array{Float64, 1}, mxiy_calcs::Array{Float64, 1}, mxdy_calcs::Array{Float64, 1}, ixmy_calcs::Array{Float64, 1},
                      dxmy_calcs::Array{Float64, 1}, dxiy_calcs::Array{Float64, 1}, ixdy_calcs::Array{Float64, 1},
                      mxmy_keys::Array{Int, 1}, mxiy_keys::Array{Int, 1}, mxdy_keys::Array{Int, 1}, ixmy_keys::Array{Int, 1},
                      dxmy_keys::Array{Int, 1}, ixdy_keys::Array{Int, 1}, dxiy_keys::Array{Int, 1})
	mxmy_back = tensor[short_node - 1, far_node - 1, 1:7, 1:64]
	mxiy_back = tensor[short_node - 1, far_node, 1:5, 1:64]
	mxdy_back = tensor[short_node - 1, far_node - 1, [1, 2, 3, 6, 7], 1:64]
	ixmy_back = tensor[short_node, far_node - 1, [1, 2, 4, 6, 7], 1:64]
	dxmy_back = tensor[short_node - 1, far_node - 1, [1, 3, 4, 5, 6], 1:64]
	ixdy_back = tensor[short_node, far_node - 1, [1, 2, 6, 7], 1:64]
	dxiy_back = tensor[short_node - 1, far_node, [1, 3, 4, 5], 1:64]

	short_mm::Float64, short_md::Float64, short_im::Float64,
	short_dm::Float64, short_dd::Float64 = state_probs(X_hmm, short_node - 1)[[1, 3, 4, 6, 7]]

	far_mm::Float64, far_md::Float64, far_im::Float64,
	far_dm::Float64, far_dd::Float64 = state_probs(Y_hmm, far_node - 1)[[1, 3, 4, 6, 7]]

	short_mi::Float64, short_ii::Float64 = state_probs(X_hmm, short_node)[[2, 5]]
	far_mi::Float64, far_ii::Float64 = state_probs(Y_hmm, far_node)[[2, 5]]

	#In all these cases we are finding the best option for the next state, and adding the cost of the state transition as well.
	mxmy_calcs[1] = maximum(mxmy_back[1, 1:end] .+ (short_mm + far_mm))
	mxmy_calcs[2] = maximum(mxmy_back[2, 1:end] .+ (short_im + far_mm))
	mxmy_calcs[3] = maximum(mxmy_back[3, 1:end] .+ (short_dm + far_mm))
	mxmy_calcs[4] = maximum(mxmy_back[4, 1:end] .+ (short_mm + far_im))
	mxmy_calcs[5] = maximum(mxmy_back[5, 1:end] .+ (short_dm + far_im))
	mxmy_calcs[6] = maximum(mxmy_back[6, 1:end] .+ (short_mm + far_dm))
	mxmy_calcs[7] = maximum(mxmy_back[7, 1:end] .+ (short_im + far_dm))
	the_keys[1], the_vals[1] = top_n(mxmy_calcs, mxmy_keys) #these are all correct.

	ixmy_calcs[1] = maximum(ixmy_back[1, 1:end] .+ (short_mi + far_mm))
	ixmy_calcs[2] = maximum(ixmy_back[2, 1:end] .+ (short_ii + far_mm))
	ixmy_calcs[3] = maximum(ixmy_back[3, 1:end] .+ (short_mi + far_im))
	ixmy_calcs[4] = maximum(ixmy_back[4, 1:end] .+ (short_mi + far_dm))
	ixmy_calcs[5] = maximum(ixmy_back[5, 1:end] .+ (short_ii + far_dm))
	the_keys[2], the_vals[2] = top_n(ixmy_calcs, ixmy_keys)

	dxmy_calcs[1]= maximum(dxmy_back[1, 1:end] .+ (short_md + far_mm))
	dxmy_calcs[2] = maximum(dxmy_back[2, 1:end] .+ (short_dd + far_mm))
	dxmy_calcs[3] = maximum(dxmy_back[3, 1:end] .+ (short_md + far_im))
	dxmy_calcs[4] = maximum(dxmy_back[4, 1:end] .+ (short_dd + far_im))
	dxmy_calcs[5] = maximum(dxmy_back[5, 1:end] .+ (short_md + far_dm))
	the_keys[3], the_vals[3] = top_n(dxmy_calcs, dxmy_keys)

	mxiy_calcs[1] = maximum(mxiy_back[1, 1:end] .+ (short_mm + far_mi))
	mxiy_calcs[2] = maximum(mxiy_back[2, 1:end] .+ (short_im + far_mi))
	mxiy_calcs[3] = maximum(mxiy_back[3, 1:end] .+ (short_dm + far_mi))
	mxiy_calcs[4] = maximum(mxiy_back[4, 1:end] .+ (short_mm + far_ii))
	mxiy_calcs[5] = maximum(mxiy_back[5, 1:end] .+ (short_dm + far_ii))
	the_keys[4], the_vals[4] = top_n(mxiy_calcs, mxiy_keys)

	dxiy_calcs[1] = maximum(dxiy_back[1, 1:end] .+ (short_md + far_mi))
	dxiy_calcs[2] = maximum(dxiy_back[2, 1:end] .+ (short_dd + far_mi))
	dxiy_calcs[3] = maximum(dxiy_back[3, 1:end] .+ (short_md + far_ii))
	dxiy_calcs[4] = maximum(dxiy_back[4, 1:end] .+ (short_dd + far_ii))
	the_keys[5], the_vals[5] = top_n(dxiy_calcs, dxiy_keys)

	mxdy_calcs[1] = maximum(mxdy_back[1, 1:end] .+ (short_mm + far_md))
	mxdy_calcs[2] = maximum(mxdy_back[2, 1:end] .+ (short_im + far_md))
	mxdy_calcs[3] = maximum(mxdy_back[3, 1:end] .+ (short_dm + far_md))
	mxdy_calcs[4] = maximum(mxdy_back[4, 1:end] .+ (short_mm + far_dd))
	mxdy_calcs[5] = maximum(mxdy_back[5, 1:end] .+ (short_im + far_dd))
	the_keys[6], the_vals[6] = top_n(mxdy_calcs, mxdy_keys)

	ixdy_calcs[1] = maximum(ixdy_back[1, 1:end] .+ (short_mi + far_md))
	ixdy_calcs[2] = maximum(ixdy_back[2, 1:end] .+ (short_ii + far_md))
	ixdy_calcs[3] = maximum(ixdy_back[3, 1:end] .+ (short_mi + far_dd))
	ixdy_calcs[4] = maximum(ixdy_back[4, 1:end] .+ (short_ii + far_dd))
	the_keys[7], the_vals[7] = top_n(ixdy_calcs, ixdy_keys)

end

function opt_tensor_core(short_hmm::types.wholeHMM, short_node::Int, far_hmm::types.wholeHMM, far_node::Int,
							assumed_letter::Int, next_letter::Int,
							short_aa_match_probs::Dict{Char, Float64}, short_aa_insert_probs::Dict{Char, Float64},
						 	far_aa_match_probs::Dict{Char, Float64}, far_aa_insert_probs::Dict{Char, Float64},
							aa7_score::Array{Float64, 2},
							codon_array::Array{Char, 1}, codon_array_rev::Array{Char, 1},
							short_aa::Char, far_aa::Char, short_aa_match::Float64,
							far_aa_match::Float64, short_aa_insert::Float64, far_aa_insert::Float64, frame::AbstractString, max_i::Int, max_j::Int)

	#for dinuc in dinucs
	flag = false
	for dinuc in 1:16
		short_aa = codon_array[assumed_letter + dinuc]
		far_aa = codon_array[(dinuc - 1)*4 + next_letter]

		if (short_node < max_i + 0)
			if short_aa == '*'
				continue
			end
			short_aa_match = short_aa_match_probs[short_aa]
			short_aa_insert = short_aa_insert_probs[short_aa]
		else
			if short_aa != '*'
				continue
			end
			short_aa_match = -0.01
			short_aa_insert = -1.5
		end

		if (far_node < max_j + 0)
			if far_aa == '*'
				continue
			end
			far_aa_match = far_aa_match_probs[far_aa]
			far_aa_insert = far_aa_insert_probs[far_aa]
		else
			if far_aa != '*'
				continue
			end
			far_aa_match = -0.01
			far_aa_insert = -1.5
		end

		if FORCE_MET && short_node == 2 && short_aa != 'M'
			continue
		end

		if FORCE_MET && far_node == 2 && far_aa != 'M'
			continue
		end #just don't consider these options... should only invoke assumed_letter is 'A' otherwise no option will work here.

		flag = true

		aa7_score[1, dinuc] = short_aa_match + far_aa_match
		aa7_score[2, dinuc] = short_aa_insert + far_aa_match
		aa7_score[3, dinuc] = short_aa_insert + far_aa_match
		aa7_score[4, dinuc] = short_aa_match + far_aa_insert
		aa7_score[5, dinuc] = short_aa_insert + far_aa_insert
		aa7_score[6, dinuc] = short_aa_match + far_aa_insert
		aa7_score[7, dinuc] = short_aa_insert + far_aa_insert
	end
	return flag
end

function initialize_4d(X_len::Int, Y_len::Int, x_init::Array{Float32, 2}, y_init::Array{Float32, 2}, min_over=50)
	tensor::Array{Float32, 4} = zeros(Float32, X_len, Y_len, 7, 64)
	tensor[1, (Y_len - min_over + 1):end, 1, 1:end] = fill(-1e12, 1, min_over, 64)
	tensor[(X_len - min_over + 1):end, 1, 1, 1:end] = fill(-1e12, min_over, 1, 64)
	if min_over != Y_len
		tensor[1, 1:(Y_len - min_over), 1, 1:end] = x_init
	end
	if min_over != X_len
		tensor[1:(X_len - min_over), 1, 1, 1:end] = y_init
	end
	tensor[1, 1:end, 2:7, 1:end] = fill(-1e12, 1, Y_len, 6, 64)
	tensor[1:end, 1, 2:7, 1:end] = fill(-1e12, X_len, 1, 6, 64)
	return tensor
end

function fill_in_tensor(X_hmm::types.wholeHMM, Y_hmm::types.wholeHMM, X_len::Int, Y_len::Int, x_init::Array{Float32, 2}, y_init::Array{Float32, 2}, X_prot::AbstractString, Y_prot::AbstractString, min_len::Int, frame::AbstractString)
	letter_tensors = Dict{AbstractString, Array{Float32, 4}}("A" => initialize_4d(X_len, Y_len, x_init, y_init, min_len),
										 "C" => initialize_4d(X_len, Y_len, x_init, y_init, min_len),
										 "G" => initialize_4d(X_len, Y_len, x_init, y_init, min_len),
										 "T" => initialize_4d(X_len, Y_len, x_init, y_init, min_len))

	maximal_scores = fill(-1.0e12, 7, 64, 4)

	codon_table = Dict{AbstractString, Char}("TGT"=>'C',"GAC"=>'D',"TTC"=>'F',"TAG"=>'*',"GTG"=>'V',"CCT"=>'P', "GCT"=>'A',"GGC"=>'G',"CGG"=>'R',
										 "ATT"=>'I',"GAT"=>'D',"CAG"=>'Q',"ATG"=>'M',"CTC"=>'L',"TCT"=>'S',"CGT"=>'R',"ACG"=>'T',"AGA"=>'R',
										 "TGG"=>'W',"TCG"=>'S',"TTA"=>'L',"AGT"=>'S',"CGA"=>'R',"TGC"=>'C',"CAA"=>'Q',"TTG"=>'L',"AAT"=>'N',
										 "AAC"=>'N',"TAA"=>'*',"TAC"=>'Y',"CCA"=>'P',"ACT"=>'T',"TAT"=>'Y',"CAC"=>'H',"CCG"=>'P',"GCA"=>'A',
										 "GAA"=>'E',"ACA"=>'T',"GTA"=>'V',"AGG"=>'R',"AGC"=>'S',"TCA"=>'S',"GGT"=>'G',"GCC"=>'A',"TGA"=>'*',
										 "GTC"=>'V',"TTT"=>'F',"CAT"=>'H',"AAG"=>'K',"AAA"=>'K',"CTT"=>'L',"CTA"=>'L',"ATC"=>'I',"ACC"=>'T',
										 "TCC"=>'S',"CCC"=>'P',"GAG"=>'E',"GCG"=>'A',"CTG"=>'L',"CGC"=>'R',"GTT"=>'V',"GGA"=>'G',"ATA"=>'I',
										 "GGG"=>'G')

	codon_array::Array{Char, 1} =	['K', 'N', 'K', 'N', 'T', 'T', 'T', 'T', 'R', 'S', 'R', 'S', 'I', 'I', 'M', 'I', 'Q', 'H', 'Q', 'H', 'P',
										 'P', 'P', 'P', 'R', 'R', 'R', 'R', 'L', 'L', 'L', 'L', 'E', 'D', 'E', 'D', 'A', 'A', 'A', 'A', 'G', 'G',
										 'G', 'G', 'V', 'V', 'V', 'V', '*', 'Y', '*', 'Y', 'S', 'S', 'S', 'S', '*', 'C', 'W', 'C', 'L', 'F', 'L', 'F']

	codon_array_rev::Array{Char, 1} = ['F', 'V', 'L', 'I', 'C', 'G', 'R', 'S', 'S', 'A', 'P', 'T', 'Y', 'D', 'H', 'N', 'L', 'V', 'L', 'M', 'W', 'G',
											'R', 'R', 'S', 'A', 'P', 'T', '*', 'E', 'Q', 'K', 'F', 'V', 'L', 'I', 'C', 'G', 'R', 'S', 'S', 'A', 'P', 'T',
											'Y', 'D', 'H', 'N', 'L', 'V', 'L', 'I', '*', 'G', 'R', 'R', 'S', 'A', 'P', 'T', '*', 'E', 'Q', 'K']

	dinucs = ["AA", "AC", "AG", "AT", "CA", "CC", "CG", "CT", "GA", "GC", "GG", "GT", "TA", "TC", "TG", "TT"]

	mxmy_keys::Array{Int, 1} = [111, 211, 311, 411, 511, 611, 711]
	mxiy_keys::Array{Int, 1} = [112, 212, 312, 412, 512]
	mxdy_keys::Array{Int, 1} = [113, 213, 313, 613, 713]
	ixmy_keys::Array{Int, 1} = [121, 221, 421, 621, 721]
	dxmy_keys::Array{Int, 1} = [131, 331, 431, 531, 631]
	ixdy_keys::Array{Int, 1} = [123, 223, 623, 723]
	dxiy_keys::Array{Int, 1} = [132, 332, 432, 532]

	aas = ['A', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'K', 'L', 'M', 'N', 'P', 'Q', 'R', 'S',
				 'T', 'V', 'W', 'Y']

	short_aa::Char = 'A'
	far_aa::Char = 'A'
	short_aa_match::Float64 = 0.0
	far_aa_match::Float64 = 0.0
	short_aa_insert::Float64 = 0.0
	far_aa_insert::Float64 = 0.0

	base_counter::Int = 1
	set_up_counter::Int = 0

	mxmy_calcs::Array{Float64, 1} = zeros(7)
	mxiy_calcs::Array{Float64, 1} = zeros(5)
	mxdy_calcs::Array{Float64, 1} = zeros(5)
	ixmy_calcs::Array{Float64, 1} = zeros(5)
	dxmy_calcs::Array{Float64, 1} = zeros(5)
	dxiy_calcs::Array{Float64, 1} = zeros(4)
	ixdy_calcs::Array{Float64, 1} = zeros(4)

	the_keys::Array{Int, 1} = zeros(Int, 7)
	the_vals::Array{Float64, 1} = zeros(Float64, 7)

	aa7_scores::Array{Float64, 2} = fill(-1000.0, 7, 16)

	min_i::Int = 2
	max_i::Int = X_len
	min_j::Int = 2
	max_j::Int = Y_len

	for i in min_i:max_i
		short_node = i

		if short_node < X_len
			short_aa_match_probs::Dict{Char, Float64} = emit_match_probs(X_hmm, short_node - 1, aas)
			short_aa_insert_probs::Dict{Char, Float64} = emit_insert_probs(X_hmm, short_node - 1, aas)
		else
			short_aa_match_probs = Dict{Char, Float64}(aa => -1e10 for aa in aas)
			short_aa_insert_probs = Dict{Char, Float64}(aa => -1e10 for aa in aas)
		end

		for j in min_j:max_j
			did_set = false
			far_node = j

			if far_node < Y_len
				far_aa_match_probs::Dict{Char, Float64} = emit_match_probs(Y_hmm, far_node - 1, aas)
				far_aa_insert_probs::Dict{Char, Float64} = emit_insert_probs(Y_hmm, far_node - 1, aas)
			else
				far_aa_match_probs = Dict{Char, Float64}(aa => -1e10 for aa in aas)
				far_aa_insert_probs = Dict{Char, Float64}(aa => -1e10 for aa in aas)
			end

			maximal_scores = fill(-1.0e12, 7, 64, 4)

			set_up_counter = 0

			for set_up in bases
				#This basically fills up the values of the_keys and the_vals with top scores for transitions to every other state.
				state_transition_calcs(letter_tensors[set_up], i, j, X_hmm, Y_hmm, the_keys, the_vals,
										 mxmy_calcs, mxiy_calcs, mxdy_calcs, ixmy_calcs, dxmy_calcs, dxiy_calcs, ixdy_calcs,
										 mxmy_keys, mxiy_keys, mxdy_keys, ixmy_keys, dxmy_keys, ixdy_keys, dxiy_keys)

				base_counter = 0
				for next_letter in 1:4
					base_counter += 1

					aa7_scores = fill(-100000.0, 7, 16)
					#Scores each of the states.
					fxn_call = opt_tensor_core(X_hmm, short_node, Y_hmm, far_node, 16 * set_up_counter, base_counter,
												short_aa_match_probs, short_aa_insert_probs, far_aa_match_probs, far_aa_insert_probs,
												aa7_scores, codon_array, codon_array_rev, short_aa, far_aa,
												short_aa_match, far_aa_match, short_aa_insert, far_aa_insert, frame, max_i, max_j)

					total_score = broadcast(+, the_vals, aa7_scores) #So this is a 7x16 matrix.

					maximal_scores[1:end, (1:16) .+ (next_letter-1)*16, set_up_counter + 1] = total_score
				end
				set_up_counter += 1
			end
			next_letter = 0
			for mat_letter in bases
				next_letter += 1
				for mat_num in 1:7
					options = vec(maximal_scores[mat_num, (1:16) .+ (next_letter-1)*16, 1:4])
					letter_tensors[mat_letter][i, j, mat_num, 1:end] = options[1:64]
				end
			end
		end
	end
	return letter_tensors
end

function sample_start_point_on_range(X_hmm, X_prot, X_len, Y_hmm, Y_prot, Y_len, num_samples, X_range, Y_range, thresh=1.0, soft_temp=1.0, first_softmax_temp=0.01, gen_bbt=false)
	#Go along all endpoints in all letters and probabilistically choose a start point based on scores.
	#NOW WITH X_RANGE and Y_RANGE TO INDICATE ACCEPTABLE PARTS OF MATRIX FROM WHICH TO SAMPLE.
	vecs_to_sample_x = Dict{AbstractString, Array{Float64, 1}}() #is the score a Float64?
	vecs_to_sample_y = Dict{AbstractString, Array{Float64, 1}}()
	letter_legal_x = Dict{AbstractString, Array{Float64, 1}}()
	letter_legal_y = Dict{AbstractString, Array{Float64, 1}}()
	best_letter_count = 0

	tail_offset = 0.0

	term_x = false
	term_y = false
	if length(X_range) == 1 && length(Y_range) > 1
		term_x = true #this means we want to go from last letter of x backward.
	elseif length(Y_range) == 1 && length(X_range) > 1
		term_y = true
	else
		return false #this will raise an error.
	end

	if !(gen_bbt == false)
		bbt = gen_bbt
	end

	MIN_STOP_OFFSET = 2

	L_values = Dict{Any, Any}()

	for letter_iterator in ["A", "C", "G", "T"]
		top_L = bbt[letter_iterator][1:end, 1:end, 1:end, 1:end]
		best_letter_count += 1
		for integer in 1:64
			top_L[1:X_len, Y_len:Y_len, 1:7, integer:integer] += reshape(repeat(tail_prot_give_log_2d(X_hmm, X_prot, X_len), 1, 7), X_len, 1, 7)
			top_L[X_len:X_len, 1:Y_len, 1:7, integer:integer] += reshape(repeat(tail_prot_give_log_2d(Y_hmm, Y_prot, Y_len), 1, 7), 1, Y_len, 7)
		end
		L_values[letter_iterator] = top_L[1:end, 1:end, 1:end, 1:end]

		final_x_points = top_L[X_len, Y_range, 1:7, 1:64]
		final_y_points = top_L[X_range, Y_len, 1:7, 1:64]

		vec_to_sample_x = vec(final_x_points)
		vec_to_sample_y = vec(final_y_points)
		vecs_to_sample_x[letter_iterator] = vec_to_sample_x
		vecs_to_sample_y[letter_iterator] = vec_to_sample_y
	end

	ADJUSTED_Y_LEN = length(Y_range)
	ADJUSTED_X_LEN = length(X_range)

	all_scores = Float64[]
	for letter_iterator in ["A", "C", "G", "T"]
		for i in vecs_to_sample_x[letter_iterator]
			if term_x
				push!(all_scores, i)
			end
		end
		for j in vecs_to_sample_y[letter_iterator]
			if term_y
				push!(all_scores, j)
			end
		end
	end

	first_steps = math.softmax(all_scores, num_samples, first_softmax_temp)
	if term_x
		sample_x_len = length(vecs_to_sample_x["A"])
		sample_y_len = 0
	elseif term_y
		sample_x_len = 0
		sample_y_len = length(vecs_to_sample_y["A"])
	end
	per_letter = sample_x_len + sample_y_len

	results = Any[]
	first_step_count = 1
	for first_step in first_steps
		first_step_count += 1

		letter_value = div(first_step - 1, per_letter) + 1
		nuc_start = string("ACGT"[letter_value])
		pos_in = ((first_step - 1) % per_letter) + 1
		alt_aa_pos = 0
		is_x = pos_in < sample_x_len
		if is_x
			aa_pos = (Y_range[1] - 1) + (((pos_in - 1) % (ADJUSTED_Y_LEN)) + 1) #maybe -1 on the Y_range business...
			alt_aa_pos = ((pos_in - 1) % (ADJUSTED_Y_LEN)) + 1
			cur_state = div((pos_in - 1) % ((ADJUSTED_Y_LEN) * 7), (ADJUSTED_Y_LEN)) + 1
			last_codon = div(pos_in - 1, ((ADJUSTED_Y_LEN) * 7)) + 1
		else
			transformed = mod((first_step - 1), per_letter) + 1 - sample_x_len
			aa_pos = (X_range[1] - 1) + (((transformed - 1) % (ADJUSTED_X_LEN)) + 1)
			alt_aa_pos = ((transformed - 1) % (ADJUSTED_X_LEN)) + 1
			cur_state = div((transformed - 1) % ((ADJUSTED_X_LEN) * 7), (ADJUSTED_X_LEN)) + 1
			last_codon = div(transformed - 1, ((ADJUSTED_X_LEN) * 7)) + 1
		end

		cur_score = all_scores[first_step]

		if is_x
			cur_score = -Float64(bbt[nuc_start][X_len, aa_pos, cur_state, last_codon])
		else
			cur_score = -Float64(bbt[nuc_start][aa_pos, Y_len, cur_state, last_codon])
		end

		if is_x
			one_path = recurse_back(X_hmm, Y_hmm, X_len, aa_pos, cur_state, cur_score, string(trinucs[last_codon][1]), trinucs[last_codon] * nuc_start, 0.0, [(X_len, aa_pos, cur_state)], L_values, thresh, soft_temp)
			push!(results, one_path)
		else
			one_path = recurse_back(X_hmm, Y_hmm, aa_pos, Y_len, cur_state, cur_score, string(trinucs[last_codon][1]), trinucs[last_codon] * nuc_start, 0.0, [(aa_pos, Y_len, cur_state)], L_values, thresh, soft_temp)
			push!(results, one_path)
		end
	end
	return results
end

function sample_start_point(X_hmm, X_prot, X_len, X_tail_log2d, Y_hmm, Y_prot, Y_len, Y_tail_log2d, num_samples, thresh=1.0, soft_temp=1.0, first_softmax_temp=0.01, gen_bbt=false)
	#Go along all endpoints in all letters and probabilistically choose a start point based on scores.
	vecs_to_sample_x = Dict{AbstractString, Array{Float64, 1}}()
	vecs_to_sample_y = Dict{AbstractString, Array{Float64, 1}}()
	letter_legal_x = Dict{AbstractString, Array{Float64, 1}}()
	letter_legal_y = Dict{AbstractString, Array{Float64, 1}}()

	tail_offset = 0.0

	if !(gen_bbt == false)
		bbt = gen_bbt
	end

	MIN_STOP_OFFSET = 2
	#This is a bit of trickery: our model has a hard time with double stop-codons (because overlapping stop codons are impossible)
	#so we set a minimum offset to keep away from the corner case of X_len, Y_len end points that are impossible to finish.
	#Should be set to 2. That works well. Setting it to 0 or 1 wouldn't work so well (at least in my experience).

	L_values = Dict{Any, Any}()

	for letter_iterator in ["A", "C", "G", "T"]
		top_L = bbt[letter_iterator][1:end, 1:end, 1:end, 1:end]
		for integer in 1:64
			top_L[1:X_len, Y_len:Y_len, 1:7, integer:integer] += reshape(repeat(X_tail_log2d, 1, 7), X_len, 1, 7)
			top_L[X_len:X_len, 1:Y_len, 1:7, integer:integer] += reshape(repeat(Y_tail_log2d, 1, 7), 1, Y_len, 7)
		end
		L_values[letter_iterator] = top_L[1:end, 1:end, 1:end, 1:end]

		final_x_points = top_L[X_len, 1:(Y_len - MIN_STOP_OFFSET), 1:7, 1:64]
		final_y_points = top_L[1:(X_len - MIN_STOP_OFFSET), Y_len, 1:7, 1:64]

		vec_to_sample_x = vec(final_x_points)
		vec_to_sample_y = vec(final_y_points)
		vecs_to_sample_x[letter_iterator] = vec_to_sample_x
		vecs_to_sample_y[letter_iterator] = vec_to_sample_y
	end

	all_scores = Float64[]
	for letter_iterator in ["A", "C", "G", "T"]
		for i in vecs_to_sample_x[letter_iterator]
			push!(all_scores, i)
		end
		for j in vecs_to_sample_y[letter_iterator]
			push!(all_scores, j)
		end
	end

	first_steps = math.softmax(all_scores, num_samples, first_softmax_temp)
	sample_x_len = length(vecs_to_sample_x["A"])
	sample_y_len = length(vecs_to_sample_y["A"])
	per_letter = sample_x_len + sample_y_len

	results = Any[]
	first_step_count = 1
	for first_step in first_steps
		first_step_count += 1

		letter_value = div(first_step - 1, per_letter) + 1
		nuc_start = string("ACGT"[letter_value])
		pos_in = ((first_step - 1) % per_letter) + 1
		is_x = pos_in < sample_x_len
		if is_x
			aa_pos = ((pos_in - 1) % (Y_len-MIN_STOP_OFFSET)) + 1
			cur_state = div((pos_in - 1) % ((Y_len-MIN_STOP_OFFSET) * 7), (Y_len - MIN_STOP_OFFSET)) + 1
			last_codon = div(pos_in - 1, ((Y_len - MIN_STOP_OFFSET) * 7)) + 1
		else
			transformed = mod((first_step - 1), per_letter) + 1 - sample_x_len
			aa_pos = ((transformed - 1) % (X_len - MIN_STOP_OFFSET)) + 1
			cur_state = div((transformed - 1) % ((X_len - MIN_STOP_OFFSET) * 7), (X_len - MIN_STOP_OFFSET)) + 1
			last_codon = div(transformed - 1, ((X_len - MIN_STOP_OFFSET) * 7)) + 1
		end

		cur_score = all_scores[first_step]

		if is_x
			cur_score = -Float64(bbt[nuc_start][X_len, aa_pos, cur_state, last_codon])
		else
			cur_score = -Float64(bbt[nuc_start][aa_pos, Y_len, cur_state, last_codon])
		end

		if is_x
			one_path = recurse_back(X_hmm, Y_hmm, X_len, aa_pos, cur_state, cur_score, string(trinucs[last_codon][1]), trinucs[last_codon] * nuc_start, 0.0, [(X_len, aa_pos, cur_state)], L_values, thresh, soft_temp)
			push!(results, one_path)
		else
			one_path = recurse_back(X_hmm, Y_hmm, aa_pos, Y_len, cur_state, cur_score, string(trinucs[last_codon][1]), trinucs[last_codon] * nuc_start, 0.0, [(aa_pos, Y_len-1, cur_state)], bbt, thresh, soft_temp)
			push!(results, one_path)
		end
	end
	return results
end

#end


text:
# Flatten the pages of a Petri Net Markup Language

"""
$(TYPEDSIGNATURES)

Merge page content into the 1st page of each pnml net.

Note that refrence nodes are still present. They can be removed later
with [`deref!`](@ref).
"""
function flatten_pages! end
flatten_pages!(model::PNML.PnmlModel) = flatten_pages!.(nets(model))

"""
$(TYPEDSIGNATURES)

Collect keys from all pages and move to first page.
"""
function flatten_pages!(net::PnmlNet)
    # Place content of subpages of 1st page before sibling page's content.
    subpages = firstpage(net).subpages
    if subpages !== nothing
        foldl(flatten_pages!, subpages; init=firstpage(net))
        empty!(subpages)
    end
    # Sibling pages.
    if length(net.pages) > 1
        foldl(flatten_pages!, net.pages[2:end]; init=firstpage(net))
        resize!(net.pages, 1)
    end
    deref!(net) # Resolve reference nodes
    return net
end

"After appending `r` to `l`, recursivly flatten `r` into `l`, then empty `r`."
function flatten_pages!(l::Page, r::Page)
    append_page!(l, r)
    if r.subpages !== nothing
        foldl(flatten_pages!, r.subpages; init=l)
    end
    r !== nothing && empty!(r)
    return l
end

"""
Append selected fields of `r` to fields of `l`.
NB: subpages are omitted from `append_page!` See [`flatten_pages!`](@ref).
Names and xml are omitted because they are scalar values, not collections.
"""
function append_page!(l::Page, r::Page;
                      keys = [:places, :transitions, :arcs,
                              :refTransitions, :refPlaces, :declaration],
                      comk = [:tools, :labels])
    @debug "append_page!($(pid(l)), $(pid(r)))"
    foreach(keys) do key
        update_maybe!(getproperty(l, key), getproperty(r, key))
    end
    # Optional fields to append.
    foreach(comk) do key
        update_maybe!(getproperty(l.com,key), getproperty(r.com,key))
    end

    l
end

# 2 Things, each could be union of nothing and `T`.
# `T`  has field `key` that is to be appended.
function update_maybe!(l::T, r::T, key::Symbol) where {T <: Maybe{Any}}
    if !isnothing(getproperty(r, key))
        if isnothing(getproperty(l, key))
            setproperty!(l, key, getproperty(r, key))
        else
            append!(getproperty(l, key), getproperty(r, key))
        end
    end
end

#! TODO test this (how would it be used?)
function update_maybe!(l::T, r::T) where {T <: Maybe{Any}}
    if !isnothing(r)
        if isnothing(l)
            l = r
        else
            append!(l, r)
        end
    end
end

"""
$(TYPEDSIGNATURES)

Remove reference nodes from arcs. Expects [`flatten_pages!`](@ref) to have
been applied so that everything is on one page (default is first page).

# Axioms
  1) All ids in a network are unique in that they only have one instance in the XML.
  2) A chain of reference Places or Transitions always ends at a Place or Transition.
  3) All ids are valid.
  4) No cycles.
"""
function deref! end

deref!(net::PnmlNet, page_idx=1) = deref!(pages(net)[page_idx])

function deref!(page::Page)
    for arc in arcs(page)
        while arc.source ∈ refplace_ids(page)
            arc.source = deref_place(page, arc.source)
        end
        while arc.target ∈ refplace_ids(page)
            arc.target = deref_place(page, arc.target)
        end
        while arc.source ∈ reftransition_ids(page)
            arc.source = deref_transition(page, arc.source)
        end
        while arc.target ∈ reftransition_ids(page)
            arc.target = deref_transition(page, arc.target)
        end
    end
    # Remove reference nodes.
    empty!(page.refPlaces)
    empty!(page.refTransitions)
    page
end

"""
$(TYPEDSIGNATURES)

Return id of referenced place.
"""
function deref_place end

deref_place(net::PnmlNet, id::Symbol, page_idx=1) = deref_place(pages(net)[page_idx], id)
deref_place(page::Page, id::Symbol) = begin
    rp = refplace(page, id)
    isnothing(rp) ? nothing : rp.ref
end

"""
$(TYPEDSIGNATURES)

Return id of referenced transition.
"""
function deref_transition end

deref_transition(net::PnmlNet, id::Symbol, page_idx=1) =
        deref_transition(net.pages[page_idx], id)
deref_transition(page::Page, id::Symbol) = begin
    rt = reftransition(page, id)
    isnothing(rt) ? nothing : rt.ref
end


text:
########## Types

"""
    Salsa.Runtime()
    Salsa.Runtime{ContextType}(context::ContextType)
    Salsa.Runtime{ContextType,ST<:AbstractSalsaStorage}(ctx::ContextType, storage::ST)

A Salsa Runtime instance contains all the state for Salsa derived functions and inputs. The
values for inputs and computed derived values are stored in a storage backend, and the
Runtime keeps track of the dependencies of a derived function, which it then communicates to
the storage backend.

Users should create a Salsa.Runtime() to store all the incremental state for a program, and
then pass it to their various `@input`s and `@derived` functions.

Inputs and Derived functions take a Runtime() instance as their first parameter.

By providing a `context` instance of `ContextType`, you can attach extra state to a Runtime,
which you can access from within your derived functions without having any affect on derived
function invalidation. (This should be used very careful because it is essentially an escape
hatch that allows you to violate the requirement that derived functions are pure
functions. The intended use is for situations where you need this context in order to
produce some results [e.g. a factory of some kind], but where the context plays no part in
determining the _value_ of the results.)
"""
abstract type Runtime{ContextType,ST<:AbstractSalsaStorage} end

struct EmptyContext end

########## Accessors

"""
    context(runtime::Runtime{CT})::CT

Return the context object that this runtime was constructed with.
"""
function context end

"""
    storage(runtime::Runtime{CT,ST})::ST

Return the storage object that this runtime was constructed with.
"""
function storage end

########## Tracing

# All runtimes except top-level ones are responsible for tracking dependencies throughout a
# region of the dataflow. Storage implementations use the following methods to access
# dependency information from a runtime.

"""
    trace(runtime::Runtime)

Return the trace object associated with this runtime.
"""
function trace end

"""
    collect_call_stack(runtime::Runtime)::Vector{DependencyKey}

Returns the stack trace of Salsa functions called on the runtime. Stack traces are only
tracked in debug mode. When debug mode is disabled, this function always returns an empty
vector.
"""
function collect_call_stack end

"""
    collect_trace(runtime::Runtime)::Vector{DependencyKey}

Returns demand (deduplicated but in-order) incurred over the lifetime of the runtime.
"""
function collect_trace end

"""
    destruct_trace!(runtime::Runtime)

Free up resources associated with the runtime. This is called by Salsa automatically when
the demand associated with that runtime is resolved (i.e. right before a result is returned
to the user).
"""
function destruct_trace! end

########## Lookups

"""
    memoized_lookup(runtime::Runtime, key::InputKey{F,TT})
    memoized_lookup(runtime::Runtime, key::DerivedKey{F,TT})

Look up `key` in the storage associated with that `runtime`, and return the (cached) value
if it exists.

For derived keys, if no valid cached value is found, compute the value from scratch by
running the derived function `F` indicated in the `DerivedKey{F}` type with the function
arguments extracted from `key.args`. The result is then cached again.
"""
function memoized_lookup end

"""
    unmemoized_input_lookup(runtime::Runtime, input, args...)
    unmemoized_input_lookup(runtime::Runtime, input_key::InputKey)

Sometimes we need to verify a new input value against the current one in user code (for
example to gracefully handle an over-deletion. This check can inherently not be memoized
(because it would lead to a cyclic dependency).
"""
function unmemoized_input_lookup end

# Branching constructor for any runtime capable of tracing. Responsible for adding the
# provided key to the trace.
function new_trace_runtime! end

########## State

"""
    new_epoch!(runtime::Runtime)

Salsa supports running with versioned storage engines. Versioned storage works with epochs
within which the dependency order must be respected. I.e. it is not allowed to assign the
same input twice in the same epoch, or to read a derived value then assign an input on which
that derived value depends.

This function allows Salsa users to start a new epoch in such situations. There are costs
associated with creating new epochs, so this functionality should be used intentionally and
sparingly.

Salsa's DefaultStorage does not support versioning. Calls to this function have no effect
when running with DefaultStorage.
"""
function new_epoch! end

"""
    previous_output(runtime::Runtime)

Must be called from within a `@derived` function. Returns the value that is currently cached
for that function and its arguments (if any). This can be used to implement incremental
logic, that reuses the previously computed value to compute the next value more efficiently.
"""
function previous_output end


text:
# ---
# title: 1669. Merge In Between Linked Lists
# id: problem1669
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Linked List
# link: <https://leetcode.com/problems/merge-in-between-linked-lists/description/>
# hidden: true
# ---
#
# You are given two linked lists: `list1` and `list2` of sizes `n` and `m`
# respectively.
#
# Remove `list1`'s nodes from the `ath` node to the `bth` node, and put `list2`
# in their place.
#
# The blue edges and nodes in the following figure incidate the result:
#
# ![](https://assets.leetcode.com/uploads/2020/11/05/fig1.png)
#
# _Build the result list and return its head._
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/05/merge_linked_list_ex1.png)
#
#
#
#     Input: list1 = [0,1,2,3,4,5], a = 3, b = 4, list2 = [1000000,1000001,1000002]
#     Output: [0,1,2,1000000,1000001,1000002,5]
#     Explanation: We remove the nodes 3 and 4 and put the entire list2 in their place. The blue edges and nodes in the above figure indicate the result.
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/05/merge_linked_list_ex2.png)
#
#
#
#     Input: list1 = [0,1,2,3,4,5,6], a = 2, b = 5, list2 = [1000000,1000001,1000002,1000003,1000004]
#     Output: [0,1,1000000,1000001,1000002,1000003,1000004,6]
#     Explanation: The blue edges and nodes in the above figure indicate the result.
#
#
#
#
# **Constraints:**
#
#   * `3 <= list1.length <= 104`
#   * `1 <= a <= b < list1.length - 1`
#   * `1 <= list2.length <= 104`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
function createPointStandards(n::Int64)
    points,weights  = computeXYZandWeights(n)
	return points,weights
end

function computeXYZandWeights(_norder::Int64)
        n = _norder
        pointsxyz = zeros(Float64,n*n,3); # Even though we only need the x and y coordinate
                                      # we still store the corresponding z coordinate
                                      # as our quadrature point actually lives on the sphere
        weights = zeros(n*n);

        # Construct Gauss quadrature
        mu,gaussweights = gausslegendre(n)

        # around z axis equidistant
        phi = [(k+0.5)*pi/n for k=0:2*n-1]

        range = 1:Int(n/2) # we only use the upper half of the sphere as quadrature point since we do pseudo three d

        # Transform between (mu,phi) and (x,y,z)
        x = sqrt.(1.0 .- mu[range].^2).*cos.(phi)'
        y = sqrt.(1.0 .- mu[range].^2).*sin.(phi)'
        z =           mu[range]    .*ones(size(phi))'
        weights = 2.0*pi/n*repeat(gaussweights[range],1,2*n)

        # assign
        pointsxyz[:,1] = x[:]
        pointsxyz[:,2] = y[:]
        pointsxyz[:,3] = z[:]

        weights = weights[:]

        return pointsxyz, weights
    end

max_stars_repo_name: vavrines/CSN

max_issues_repo_name: vavrines/CSN

max_forks_repo_name: vavrines/CSN


text:
# for using keyword_vgh()
# change the order of parameters.

# function ComputeLL(LLs::SharedArray{Float64,1}, ratdata, ntrials::Int#, args, x::Vector{T})
function ComputeLL(LLs, ratdata, ntrials::Int#, args, x::Vector{T})
    ;kwargs...)
    LL = 0.

    @sync @parallel for i in 1:ntrials
        RightClickTimes, LeftClickTimes, maxT, rat_choice = TrialData(ratdata, i)
        Nsteps = Int(ceil(maxT/dt))

        # LLs[i] = LogLikelihood(params, RightClickTimes, LeftClickTimes, Nsteps, rat_choice)
        LLs[i] = LogLikelihood(RightClickTimes, LeftClickTimes, Nsteps, rat_choice
                ;kwargs...)#make_dict(args, x)...)
    end

    LL = -sum(LLs)
    return LL
end

function ComputeLL_bbox(LLs, ratdata, ntrials::Int#, args, x::Vector{T})
    ;kwargs...)
    LL  = 0.
    (k,v) = kwargs[1]
    LLs = SharedArray(typeof(v), ntrials)#zeros(eltype(params),ntrials)

    @sync @parallel for i in 1:ntrials
        RightClickTimes, LeftClickTimes, maxT, rat_choice = TrialData(ratdata, i)
        Nsteps = Int(ceil(maxT/dt))
        LLs[i] = LogLikelihood(RightClickTimes, LeftClickTimes, Nsteps, rat_choice
            ;kwargs...)#;make_dict(args, x)...)
    end
    LL = -sum(LLs)
    return LL
end

function ComputeGrad{T}(ratdata, ntrials::Int, args, x::Vector{T})
    LL        = 0.
    LLgrad    = zeros(T,length(x))

    # do we still need wrapper?
    function WrapperLL(;kwargs...)#(params::Vector{T})
        LL  = 0.
        (k,v) = kwargs[1]
        LLs = SharedArray(typeof(v), ntrials)#zeros(eltype(params),ntrials)

        @sync @parallel for i in 1:ntrials
            RightClickTimes, LeftClickTimes, maxT, rat_choice = TrialData(ratdata, i)
            Nsteps = Int(ceil(maxT/dt))
            LLs[i] = LogLikelihood(RightClickTimes, LeftClickTimes, Nsteps, rat_choice
                ;kwargs...)#;make_dict(args, x)...)
        end
        LL = -sum(LLs)
        return LL
    end

    do_hess = false

    # LL = WrapperLL(;make_dict(args, x)...)
    # println(LL)
    LL, LLgrad = GeneralUtils.keyword_vgh((;params...)
        -> WrapperLL(;params...), args, x, do_hess)

    # LL, LLgrad = GeneralUtils.keyword_vgh(WrapperLL, args, x ,do_hess)

    # LL, LLgrad = GeneralUtils.vgh(WrapperLL, params, do_hess)

    # result =  DiffBase.GradientResult(params)
    # ForwardDiff.gradient!(result, WrapperLL, params);
    # LL     = DiffBase.value(result)
    # LLgrad = DiffBase.gradient(result)
    return LL, LLgrad
end

function ComputeHess{T}(ratdata, ntrials::Int, args, x::Vector{T})
    LL        = 0.
    LLgrad    = zeros(T,length(x))
    LLhess    = zeros(T,length(x),length(x))

    function WrapperLL(;kwargs...)#(params::Vector{T})
        LL  = 0.
        (k,v) = kwargs[1]
        LLs = SharedArray(typeof(v), ntrials)#zeros(eltype(params),ntrials)

        @sync @parallel for i in 1:ntrials
            RightClickTimes, LeftClickTimes, maxT, rat_choice = TrialData(ratdata, i)
            Nsteps = Int(ceil(maxT/dt))
            LLs[i] = LogLikelihood(RightClickTimes, LeftClickTimes, Nsteps, rat_choice
                ;kwargs...)#;make_dict(args, x)...)

        end
        LL = -sum(LLs)
        return LL
    end

    do_hess = true

    # LL = WrapperLL(;make_dict(args, x)...)
    # println(LL)
    LL, LLgrad, LLhess = GeneralUtils.keyword_vgh((;params...)
        -> WrapperLL(;params...), args, x, do_hess)

    # LL, LLgrad, LLhess = GeneralUtils.vgh(WrapperLL, params)

    # result =  DiffBase.HessianResult(params)
    # ForwardDiff.hessian!(result, WrapperLL, params);
    # LL     = DiffBase.value(result)
    # LLgrad = DiffBase.gradient(result)
    # LLhess = DiffBase.hessian(result)
    return LL, LLgrad, LLhess
end

function TrialsLikelihood{T}(LL::AbstractArray{T,1}, ratdata, ntrials::Int
     ;kwargs...)
    for i in 1:ntrials
        RightClickTimes, LeftClickTimes, maxT, rat_choice = TrialData(ratdata, i)
        Nsteps = Int(ceil(maxT/dt))

        LL[i] = LogLikelihood(RightClickTimes, LeftClickTimes, Nsteps, rat_choice
                            ;kwargs...)
    end
end


text:
"""
## Indexing
`pa` is an instance of `SeisForwExpt`.
* `pa[:data,i]` : get data for ith supersource (defaults to first), after `update!`
* `pa[:ageom]` : get `ageom` originally input while creating `pa`
* `pa[:srcwav]`: `srcwav` used to create `pa`
* `pa[:snaps,i]` : snapshots corresponding to `tsnaps` and ith source, after `update!`
"""
function Base.getindex(pa::PFdtd, s::Symbol, iss::Int = 1)
    @assert iss ≤ length(pa.c.ageom[1])
    @assert s in [:snaps, :data, :ageom, :srcwav]
    if (s == :snaps)
        ip = findall(in.(iss, pa.sschunks))[1]
        issp = findall(x -> x == iss, pa.sschunks[ip])[]
        @sync begin
            p = procs(pa.p)[ip]
            @sync remotecall_wait(p) do
                pap = localpart(pa.p)
				return pap[1].ss[issp].snaps
            end
        end
    elseif (s == :data)
        (iss == 0) ? (return pa.c.data[1]) : (return pa.c.data[1][iss])
    elseif (s == :ageom)
        return pa.c.ageom
    elseif (s == :srcwav)
        return pa.c.srcwav
    end
end

function Base.show(io::Base.IO, pa::PFdtd)
    println(typeof(pa), "")
    println("pa[:snaps]")
    println("pa[:data] : modeled data after running `update!`")
    println("pa[:ageom]")
    println("pa[:srcwav]")
end

Base.show(io::Base.IO, pa::Vector{P_x_worker_x_pw{N}}) where {N} = nothing
Base.show(io::Base.IO, pa::P_x_worker_x_pw) = nothing
Base.show(io::Base.IO, pa::P_x_worker_x_pw_x_ss) = nothing


text:
using Test
using SpikingNeuralNetworks
using Plots
SNN = SpikingNeuralNetworks
SNN.@load_units

for tests in ["spike_current_search_test.jl"]
    include(tests)
end


text:
const es_1_a_29 =[
"un"
"dos"
"tres"
"cuatro"
"cinco"
"seis"
"siete"
"ocho"
"nueve"
"diez"
"once"
"doce"
"trece"
"catorce"
"quince"
"dieciséis"
"diecisiete"
"dieciocho"
"diecinueve"
"veinte"
"veintiun"
"veintidós"
"veintitrés"
"veinticuatro"
"veinticinco"
"veintiséis"
"veintisiete"
"veintiocho"
"veintinueve"
]

const es_decenas = [
    ""
    ""
    "treinta"
    "cuarenta"
    "cincuenta"
    "sesenta"
    "setenta"
    "ochenta"
    "noventa"
]

const es_centenas = [
"ciento"
"doscientos"
"trescientos"
"cuatrocientos"
"quinientos"
"seiscientos"
"setecientos"
"ochocientos"
"novecientos"
]

const es_multiplos_1e6_singular = [
    "millón"
    "billón"
    "trillón"
    "cuatrillón"
    "quintillón"
    "sextillón"
    "septillón"
    "octillón"
    "nonillón"
    "decillón"
    "undecillón"
    "duodecillón"
]

const es_multiplos_1e6_plural = [
    "millones"
    "billones"
    "trillones"
    "cuatrillones"
    "quintillones"
    "sextillones"
    "septillones"
    "octillones"
    "nonillones"
    "decillones"
    "undecillones"
    "duodecillones"
]

const es_decimales = [
    "décimo"
    "centécimo"
    "milésimo"
    "diezmilésimo"
    "cienmilésimo"
    "millonésimo"
]


text:
function mc2b!(mc::AbstractVector, α)
    b = mc

    for i=length(mc)-1:-1:1
        @inbounds b[i] = b[i] - α*b[i+1]
    end
    b
end

mc2b(mc::AbstractVector, α) = mc2b!(copy(mc), α)

function mc2b!{T<:MelCepstrum}(state::SpectralParamState{T})
    assert_not_ready_to_filt(state)

    α = allpass_alpha(paramdef(state))
    state.ready_to_filt = true
    α == zero(α) && return state

    data = rawdata(state)
    state.data = mc2b!(data, α)
    state
end

function mc2b!{T<:MelGeneralizedCepstrum}(state::SpectralParamState{T})
    throw(ArgumentError("""
                        use `mgc2b` instead of `mc2b` for mel-generalized cepstrums"""))
end

function mc2b{T<:MelGeneralizedCepstrum}(state::SpectralParamState{T})
    mc2b!(copy(state))
end


text:
__precompile__()

module EzExcel

@static if VERSION < v"0.7.0-DEV.2005"
    using Missings
end
using EzXML
using ZipFile
using ZipFile: ReadableFile, Reader

include("parse.jl")
include("index.jl")
include("format.jl")

include("workbook.jl")
include("worksheet.jl")
include("cell.jl")
include("show.jl")

export WorkBook, WorkSheet, Cell,
       sheetname, sheets,
       size, length

end


text: Dict{String,Any}("electron_mass" => Dict{String,Any}("value" => "9.109109383e-31","type" => "float"))

max_stars_repo_name: SmiteWindows/julia

max_issues_repo_name: SmiteWindows/julia

max_forks_repo_name: SmiteWindows/julia


text:
include("../functions/Functions.jl")
include("../evaluate/Evaluate.jl")
include("../opt/pgd.jl")
include("../opt/agd.jl")
include("../opt/GD-const.jl")
include("../data/data.jl")
include("../plot/Plotting.jl")

using Plotting

X, y = getOldData()

x0 = zeros(size(X,2),1)

lambda = 0.000001

iter = 50
step = 0.01

f = f_logreg(X,y,x0,lambda)
g = g_logreg(X,y,x0,0)

g_l2 = g_logreg(X,y,x0,lambda)

prox_operator = SqrNormL2(lambda/2)

@time (x_pgd,state,val_array_pgd,stop_array_pgd) = pgd(X,y,x0,f,g,1e-8,lambda,iter,step,size(X,1),prox_operator)
@time (x_agd,state,val_array_agd,stop_array_agd) = agd(X,y,x0,f,g_l2,1e-8,lambda,iter,step,size(X,1),prox_operator)
@time (x_gd,state,val_array_gd,stop_array_gd) = gd(X,y,x0,f,g_l2,1e-8,lambda,iter,step,size(X,1))

plotTripleArray(val_array_pgd,val_array_agd,val_array_gd, "GD-FVAL-COMPARE.svg")
plotTripleArray(stop_array_pgd,stop_array_agd,stop_array_gd, "GD-STOP-COMPARE.svg")

max_stars_repo_name: zagazao/num-opt

max_issues_repo_name: zagazao/num-opt

max_forks_repo_name: zagazao/num-opt


text:
export TwoFactorRepresentation
"""
    X = U*Z'

    Two factor factorization of a low rank matrix into factors U ∈ ℝⁿˣʳ and Z ∈ ℝᵐˣʳ.
    U should span the range of X while Z spans the co-range. The factorization is non-unique and
    U need not be orthogonal! However, U is often chosen orthogonal for cheap (pseudo)inversion. In that case,
    note that orthogonality of U is not guaranteed to and in fact will rarely be preserved under the operations that
    are supported (multiplication, addition, etc.). In order to reorthonormalize U, simply call `orthonormalize!(X, alg)`
    where alg refers to the algorithm used to compute the orthonormalization:
        * GradientDescent()
        * QRFact()
        * SVDFact()
"""
mutable struct TwoFactorRepresentation{uType, zType} <: AbstractLowRankRepresentation
    U::uType
    Z::zType
end

rank(LRA::TwoFactorRepresentation) = size(LRA.U, 2)
size(LRA::TwoFactorRepresentation) = (size(LRA.U,1), size(LRA.Z,1))
size(LRA::TwoFactorRepresentation, ::Val{1}) = size(LRA.U,1)
size(LRA::TwoFactorRepresentation, ::Val{2}) = size(LRA.Z,1)
size(LRA::TwoFactorRepresentation, i::Int) = size(LRA, Val(i))
Matrix(LRA::TwoFactorRepresentation) = LRA.U*LRA.Z'
getindex(LRA::TwoFactorRepresentation, i::Int, j::Int) = sum(LRA.U[i,k]*LRA.Z[j,k] for k in 1:rank(LRA)) # good enough for now
getindex(LRA::TwoFactorRepresentation, i, j::Int) = TwoFactorRepresentation(LRA.U[i,:], LRA.Z[[j],:])
getindex(LRA::TwoFactorRepresentation, i::Int, j) = TwoFactorRepresentation(LRA.U[[i],:], LRA.Z[j,:])
getindex(LRA::TwoFactorRepresentation, i, j) = TwoFactorRepresentation(LRA.U[i,:], LRA.Z[j,:])
getindex(LRA::TwoFactorRepresentation, ::Colon, j::AbstractVector) = TwoFactorRepresentation(LRA.U, LRA.Z[j,:])
getindex(LRA::TwoFactorRepresentation, i::AbstractVector, ::Colon) = TwoFactorRepresentation(LRA.U[i,:], LRA.Z)
getindex(LRA::TwoFactorRepresentation, ::Colon, j::Int) = TwoFactorRepresentation(LRA.U, LRA.Z[[j],:])
getindex(LRA::TwoFactorRepresentation, i::Int, ::Colon) = TwoFactorRepresentation(LRA.U[[i],:], LRA.Z)
hcat(A::TwoFactorRepresentation, B::TwoFactorRepresentation) = TwoFactorRepresentation(hcat(A.U, B.U), blockdiagonal(A.Z, B.Z))
vcat(A::TwoFactorRepresentation, B::TwoFactorRepresentation) = TwoFactorRepresentation(blockdiagonal(A.U, B.U), hcat(A.Z, B.Z))

# simple support of adjoints, probably not ideal though
adjoint(LRA::TwoFactorRepresentation) = TwoFactorRepresentation(conj(LRA.Z),conj(LRA.U))

## *
*(A::AbstractMatrix, B::TwoFactorRepresentation) = TwoFactorRepresentation(A*B.U, B.Z)
*(A::TwoFactorRepresentation, B::AbstractMatrix) = TwoFactorRepresentation(A.U, B'*A.Z)
*(A::TwoFactorRepresentation, ::UniformScaling) = A
*(::UniformScaling, A::TwoFactorRepresentation) = A
function *(A::TwoFactorRepresentation,B::TwoFactorRepresentation)
    if rank(A) ≤ rank(B) # minimize upper bound on rank
        return TwoFactorRepresentation(A.U, B.Z*(B.U'*A.Z))
    else
        return TwoFactorRepresentation(A.U*(A.Z'*B.U), B.Z)
    end
end
*(A::TwoFactorRepresentation, α::Number) = TwoFactorRepresentation(A.U, α*A.Z)
*(α::Number, A::TwoFactorRepresentation) = TwoFactorRepresentation(A.U, α*A.Z)
*(A::TwoFactorRepresentation, v::AbstractVector) = A.U*(A.Z'*v)
*(v::AbstractVector, A::TwoFactorRepresentation) = (v*A.U)*A.Z'

## +
+(A::TwoFactorRepresentation, B::TwoFactorRepresentation) = TwoFactorRepresentation(hcat(A.U, B.U), hcat(A.Z, B.Z))
+(A::AbstractMatrix, B::TwoFactorRepresentation) = A + Matrix(B)
+(A::TwoFactorRepresentation, B::AbstractMatrix) = Matrix(A) + B
-(A::TwoFactorRepresentation, B::TwoFactorRepresentation) = TwoFactorRepresentation(hcat(A.U, B.U), hcat(A.Z, -B.Z))
-(A::TwoFactorRepresentation, B::AbstractMatrix) = Matrix(A) - B
-(A::AbstractMatrix, B::TwoFactorRepresentation) = A - Matrix(B)

## .*
function hadamard(A::TwoFactorRepresentation, B::TwoFactorRepresentation)
    @assert size(A) == size(B) "elementwise product is only defined between matrices of equal dimension"
    rA, rB = rank(A), rank(B)
    r_new = rA*rB
    U = ones(eltype(A.U), size(A,1), r_new)
    Z = ones(eltype(A.Z), size(A,2), r_new)
    AUcols = [@view A.U[:,i] for i in 1:rA]
    AZcols = [@view A.Z[:,i] for i in 1:rA]
    BUcols = [@view B.U[:,i] for i in 1:rB]
    BZcols = [@view B.Z[:,i] for i in 1:rB]
    k = 0
    for r1 in 1:rA, r2 in 1:rB
        k += 1
        U[:,k] = AUcols[r1] .* BUcols[r2]
        Z[:,k] = AZcols[r1] .* BZcols[r2]
    end
    return TwoFactorRepresentation(U,Z)
end

## .^
function elpow(A::TwoFactorRepresentation, d::Int)
    @assert d >= 1 "elementwise power operation 'elpow' only defined for positive powers"
    r = rank(A)
    r_new = binomial(r+d-1, d)
    Ucols = [@view A.U[:,i] for i in 1:r]
    Zcols = [@view A.Z[:,i] for i in 1:r]
    U = ones(eltype(A.U), size(A,1), r_new)
    Z = ones(eltype(A.Z), size(A,2), r_new)
    k = 0
    for exps in multiexponents(r,d)
        k += 1
        Z[:, k] .*= multinomial(exps...)
        for j in 1:r
            if exps[j] > 0
                U[:, k] .*= Ucols[j].^exps[j]
                Z[:, k] .*= Zcols[j].^exps[j]
            end
        end
    end
    return TwoFactorRepresentation(U, Z)
end

function elpow(A,d::Int)
    return A.^d
end

## special cases
function add_to_cols(LRA::TwoFactorRepresentation, v::AbstractVector)
    return LRA + TwoFactorRepresentation(v, ones(eltype(v), size(LRA, 2)))
end

function multiply_cols(LRA::TwoFactorRepresentation, v::AbstractVector)
    return hadamard(LRA, TwoFactorRepresentation(v, ones(eltype(v), size(LRA, 2))))
end

function add_to_cols(A, v::AbstractVector)
    return A .+ v
end

function multiply_cols(A, v::AbstractVector)
    return v .* A
end

function add_to_rows(LRA::TwoFactorRepresentation, v::AbstractVector)
    return LRA + TwoFactorRepresentation(ones(eltype(v), size(LRA, 1)), v)
end

function multiply_rows(LRA::TwoFactorRepresentation, v::AbstractVector)
    return hadamard(LRA, TwoFactorRepresentation(ones(eltype(v), size(LRA, 1)), v))
end

function add_to_rows(A, v)
    return A .+ v'
end

function multiply_rows(A, v)
    return A .* v'
end

function add_scalar(LRA::TwoFactorRepresentation, α::Number)
    return LRA + TwoFactorRepresentation(ones(eltype(α), size(LRA, 1)), α*ones(eltype(α), size(LRA,2)))
end

function add_scalar(A, α::Number)
    return A .+ α
end

## rounding
function svd(A::TwoFactorRepresentation, alg=QRFact())
    orthonormalize!(A, alg)
    U_, S_, V_ = svd(A.Z)
    return SVDLikeRepresentation(A.U*U_, S_, V_)
end

function qr(A::TwoFactorRepresentation, alg=QRFact())
    orthonormalize!(A, alg)
    Q_, R_ = qr(A.Z')
    return TwoFactorRepresentation(A.U*Q_, R_')
end

function round(A::TwoFactorRepresentation, alg = SVDFact(); tol = sqrt(eps(eltype(A.U))), rmin::Int = 1, rmax::Int = rank(A), alg_orthonormalize = QRFact())
    orthonormalize!(A, alg_orthonormalize)
    Z_lr = truncated_svd(A.Z, alg; tol=tol, rmin = rmin, rmax = rmax)
    return TwoFactorRepresentation(A.U*Z_lr.V, Z_lr.U*Z_lr.S)
end

function round(A::TwoFactorRepresentation, rank::Int, alg = SVDFact(); alg_orthonormalize = QRFact())
    orthonormalize!(A, alg_orthonormalize)
    Z_lr = truncated_svd(A.Z, rank, alg)
    return TwoFactorRepresentation(A.U*Z_lr.V, Z_lr.U*Z_lr.S)
end

## orthonormalization
function orthonormalize!(LRA::TwoFactorRepresentation, alg)
    orthonormalize!(LRA.U, LRA.Z, alg)
end


text:
mtcars = dataset("datasets", "mtcars")

priors = [
    Prior_Ridge(),
    Prior_Laplace(),
    Prior_Cauchy(),
    Prior_TDist(),
    Prior_Uniform(),
]

model = @fitmodel((MPG ~ HP + WT + Gear), mtcars, LinearRegression())
@test sizeof(model) > 0

for prior in priors
    model = @fitmodel((MPG ~ HP + WT + Gear), mtcars, LinearRegression(), prior)
    @test sizeof(model) > 0
end


text:
struct EconomicDispatch <: AbstractOperationsModel end
struct SCEconomicDispatch <: AbstractOperationsModel end

function EconomicDispatch(system::PSY.System, transmission::Type{S}; optimizer::Union{Nothing,JuMP.OptimizerFactory}=nothing, kwargs...) where {S <: PM.AbstractPowerFormulation}

    devices = Dict{Symbol, PSI.DeviceModel}(:ThermalGenerators => PSI.DeviceModel(PSY.ThermalGen, PSI.ThermalDispatch),
                                            :RenewableGenerators => PSI.DeviceModel(PSY.RenewableGen, PSI.RenewableFullDispatch),
                                            :Loads => PSI.DeviceModel(PSY.PowerLoad, PSI.StaticPowerLoad))

    branches = Dict{Symbol, PSI.DeviceModel}(:Lines => PSI.DeviceModel(PSY.Branch, PSI.SeriesLine))
    services = Dict{Symbol, PSI.ServiceModel}(:Reserves => PSI.ServiceModel(PSY.Reserve, PSI.AbstractReservesForm))

    return PowerOperationModel(EconomicDispatch,
                                   transmission,
                                    devices,
                                    branches,
                                    services,
                                    system,
                                    optimizer = optimizer; kwargs...)

end

function SCEconomicDispatch(system::PSY.System; optimizer::Union{Nothing,JuMP.OptimizerFactory}=nothing, kwargs...)

    if :PTDF in keys(kwargs)

        PTDF = kwargs[:PTDF]

    else
        @info "PTDF matrix not provided. It will be constructed using PowerSystems.PTDF"
        PTDF, A = PowerSystems.buildptdf(system.branches, system.buses);
    end

    devices = Dict{Symbol, PSI.DeviceModel}(:ThermalGenerators => PSI.DeviceModel(PSY.ThermalGen, PSI.ThermalDispatch),
    :RenewableGenerators => PSI.DeviceModel(PSY.RenewableGen, PSI.RenewableFullDispatch),
    :Loads => PSI.DeviceModel(PSY.PowerLoad, PSI.StaticPowerLoad))

    branches = Dict{Symbol, PSI.DeviceModel}(:Lines => PSI.DeviceModel(PSY.Branch, PSI.SeriesLine))
    services = Dict{Symbol, PSI.ServiceModel}(:Reserves => PSI.ServiceModel(PSY.Reserve, PSI.AbstractReservesForm))

    return PowerOperationModel(EconomicDispatch,
                                    StandardPTDFForm,
                                    devices,
                                    branches,
                                    services,
                                    system,
                                    optimizer = optimizer; PTDF = PTDF, kwargs...)

end


text:
using Base.Iterators: flatten, repeated
using PrettyTables: pretty_table, simple

struct Event{T, U, V}
    time::T
    tag::U
    val::V
end

Base.show(io::IO, e::Event) = print(io, "time: $(e.time),\ttag:$(e.tag),\tvalue:$(e.val)")

struct SkipValidation end # indicates no validation of data in EventTS constructor

struct EventTS{T, U, V} <: AbstractVector{Event{T}}
    timestamps::Vector{T}
    tag::U
    values::V
    function EventTS(timestamps::Vector{T}, tag::U, values::V, ::SkipValidation) where {T, U, V}
        new{T,U,V}(timestamps, tag, values)
    end
end

function EventTS(;timestamps, tag, values)
    @assert issorted(timestamps)
    @assert length(timestamps) == length(values)
    _validate_tag(timestamps, tag)
    EventTS(timestamps, tag, values, SkipValidation())
end

_validate_tag(timestamps, tag) = _validate_tag(timestamps, tag, tagtype(tag))
_validate_tag(timestamps, tag, ::EventTag) = nothing
function _validate_tag(timestamps, tag, ::SeriesTag)
    @assert length(timestamps) == length(tag)
    @assert tagtype(first(tag)) == EventTag()
end

timestamps(ts::EventTS) = (e.time for e in ts)
timestamps(ts::EventTS, tag) = (e.time for e in ts if e.tag==tag)

Base.values(ts::EventTS) = (e.val for e in ts)
Base.values(ts::EventTS, tag) = (e.val for e in ts if e.tag==tag)

tag(ts::EventTS) = ts.tag
tags(ts::EventTS) = _tags(ts::EventTS, tagtype(ts))
_tags(ts::EventTS, ::EventTag) = repeated(ts.tag, length(ts))
_tags(ts::EventTS, ::SeriesTag) = ts.tag

tagtype(s::EventTS) = tagtype(s.tag)
duration(ts::EventTS) =  ts.timestamps[end] - ts.timestamps[1]

function _2matrix(ts::EventTS)
    mat = Matrix{Any}(undef, length(ts), 3)
    mat[:,1] .= timestamps(ts)
    mat[:,2] .= tags(ts)
    mat[:,3] .= ts.values
    mat
end

pretty_print(ts::EventTS) = pretty_table(_2matrix(ts), [:time, :tag, :value]; tf=simple)

###
### Abstract array Interface
###

Base.IndexStyle(::EventTS) = IndexLinear()
Base.size(ts::EventTS) = size(ts.timestamps)
Base.setindex!(::EventTS, v, i) = @error "EventTS does not support set_index!"
Base.getindex(ts::EventTS, I) = _getindex(ts, I, tagtype(ts))

function _getindex(ts::EventTS, s, ::EventTag)
    EventTS(ts.timestamps[s], ts.tag, ts.values[s], SkipValidation())
end
function _getindex(ts::EventTS, s, ::SeriesTag)
    EventTS(ts.timestamps[s], ts.tag[s], ts.values[s], SkipValidation())
end
function _getindex(ts::EventTS, i::Integer, ::EventTag)
    Event(ts.timestamps[i], ts.tag, ts.values[i])
end
function _getindex(ts::EventTS, i::Integer, ::SeriesTag)
    Event(ts.timestamps[i], ts.tag[i], ts.values[i])
end

splice(ts::EventTS{T}...) where {T} = foldl(_splice, ts)

function _timestamps_values_sortperm(ts1, ts2)
    timestamps = flatten((ts1.timestamps, ts2.timestamps)) |> collect
    values = flatten((ts1.values, ts2.values)) |> collect
    perm = sortperm(timestamps)
    timestamps[perm], values[perm], perm
end

_splice(ts1::EventTS, ts2::EventTS) = _splice(ts1, tagtype(ts1), ts2, tagtype(ts2))
function _splice(ts1, ::EventTag, ts2, ::EventTag)
    timestamps, values, perm = _timestamps_values_sortperm(ts1, ts2)
    if ts1.tag == ts2.tag
        return EventTS(timestamps, ts1.tag, values, SkipValidation())
    else
        tag = flatten((repeated(ts1.tag, length(ts1)), repeated(ts2.tag, length(ts2)))) |> collect
        return EventTS(timestamps, tag[perm], values, SkipValidation())
    end
end
_splice(ts1, ::EventTag, ts2, ::SeriesTag) = _splice(ts2, tagtype(ts2), ts1, tagtype(ts1))
function _splice(ts1, ::SeriesTag, ts2, ::EventTag)
    timestamps, values, perm = _timestamps_values_sortperm(ts1, ts2)
    tag2 = repeated(ts2.tag, length(ts2)) |> collect
    tag = [ts1.tag; tag2]
    EventTS(timestamps, tag[perm], values, SkipValidation())
end
function _splice(ts1, ::SeriesTag, ts2, ::SeriesTag)
    timestamps, values, perm = _timestamps_values_sortperm(ts1, ts2)
    tag = vcat(ts1.tag, ts2.tag)
    EventTS(timestamps, tag[perm], values, SkipValidation())
end

Base.split(ts::EventTS) = split(ts, tagtype(ts))
Base.split(ts::EventTS, ::EventTag) = [ts,]
function Base.split(ts::EventTS, ::SeriesTag)
    tags_ =  tags(ts) |> unique |> sort
    [EventTS(timestamps(ts, t) |> collect, t,
             values(ts,t) |> collect, SkipValidation()) for t in tags_]
end

drop_repeated(ts::EventTS; keep_end=true) = drop_repeated(ts, tagtype(ts), keep_end=keep_end)
function drop_repeated(ts::EventTS, ::SeriesTag; keep_end)
    splice([drop_repeated(ts_, keep_end=keep_end && tag(ts_)==ts.tag[end])
           for ts_ in split(ts)]...)
end
function drop_repeated(ts::EventTS, ::EventTag; keep_end)
    select = [true; [a!=b for (a,b) in neighbors(ts.values)]]
    keep_end && (select[end] = true)
    EventTS(ts.timestamps[select], ts.tag, ts.values[select], SkipValidation())
end

function _combined_val(tag, val, utags, vals)
    idx = searchsortedfirst(utags, tag)
    out = Tuple(i==idx ? val : v for (i,v) in enumerate(vals))
    vals[idx] = val
    out
end

merge_tags(ts::EventTS) = merge_tags(ts, tagtype(ts))

merge_tags(ts::EventTS, ::EventTag) = ts
function merge_tags(ts::EventTS, ::SeriesTag)
    utags =  tags(ts) |> unique |> sort
    vals = Any[nothing for i in  1:length(utags)]
    values = [_combined_val(t, v, utags, vals) for (t,v) in zip(ts.tag, ts.values)]
    EventTS(ts.timestamps, Tuple(utags), values, SkipValidation())
end


text:
"""
```
struct MMADualObj{TApprox <: AbstractMMAApprox} <: AbstractFunction
    model::MMAApproxModel{TApprox}
    x::AbstractVector # primal solution
    approxfg::AbstractVector # approximate [f(x); g(x)]
end
```

The dual objective function of the MMA model, `model`. The approximate objective and constraint values are cached in `approxfg` when the dual objective is called. For every input `λ`, the optimal primal solution is computed and cached in the field `x`.
"""
@params struct MMADualObj{TApprox <: AbstractMMAApprox} <: AbstractFunction
    model::MMAApproxModel{TApprox}
    x::AbstractVector # primal solution
    approxfg::AbstractVector # approximate [f(x); g(x)]
end

"""
    MMADualObj(model::MMAApproxModel)

Constructs the dual objective function from the approximate MMA model, `model`.
"""
function MMADualObj(model::MMAApproxModel{<:MMAApprox})
    xk = getxk(model.approx_objective_ineq_constraints)
    x = copy(xk)
    approxfg = copy(getfk(model.approx_objective_ineq_constraints))
    return MMADualObj(model, x, approxfg)
end
function MMADualObj(model::MMAApproxModel{<:XMMAApprox})
    xk = getxk(model.approx_objective_ineq_constraints)
    yk = model.approx_objective_ineq_constraints.y
    zk = model.approx_objective_ineq_constraints.z[]
    x = [xk; yk; zk]
    approxfg = copy(getfk(model.approx_objective_ineq_constraints))
    return MMADualObj(model, x, approxfg)
end

getdim(::MMADualObj) = 1

# To be maximized
"""
    (f::MMADualObj)(λ)

Evaluates the dual objective value at the dual solution `λ`. The optimal primal solution is computed and used automatically.
"""
function (f::MMADualObj)(λ)
    x = optimizeprimal!(f, λ)
    return f(x, λ)
end

"""
    (f::MMADualObj)(x, λ)

Evaluates the dual objective value at the primal solution `x` and the dual solution `λ`.
"""
function (f::MMADualObj)(x, λ)
    temp = f.model.approx_objective_ineq_constraints(x)
    return temp[1] + dot(temp[2:end], λ)
end

"""
    optimizeprimal!(f::MMADualObj, λ::AbstractVector)

Finds and returns the primal minimizer of the Lagrangian function for a given dual solution `λ`.
"""
function optimizeprimal!(f::MMADualObj{TA}, λ::AbstractVector) where {TA}
    @unpack model, x = f
    @unpack p, q, l, u, r, ρ = getmmaapprox(model.approx_objective_ineq_constraints)
    min = getmin(model)
    max = getmax(model)
    T = eltype(x)

    map!(view(x, 1:size(p, 2)), 1:size(p, 2)) do j
        @views t1 = p[1, j] + λ' * p[2:end, j]
        @views t2 = q[1, j] + λ' * q[2:end, j]
        temp1 = t1 / (u[j] - min[j])^2 - t2 / (min[j] - l[j])^2
        temp2 = t1 / (u[j] - max[j])^2 - t2 / (max[j] - l[j])^2
        if isfinite(min[j]) && temp1 >= 0
            return min[j]
        elseif isfinite(max[j]) && temp2 <= 0
            return max[j]
        else
            return (sqrt(t1) * l[j] + sqrt(t2) * u[j]) / (sqrt(t1) + sqrt(t2))
        end
    end

    if TA <: XMMAApprox
        @unpack a, c, d = model.approx_objective_ineq_constraints
        map!(view(x, size(p, 2) + 1 : size(p, 2) + length(c)), 1:length(c)) do i
            temp = c[i] - λ[i]
            if temp >= 0
                return zero(T)
            else
                return -temp / d[i]
            end
        end
        temp = a[1] - dot(a[2:end], λ)
        ε = 1e-4
        if temp >= 0
            x[end] = zero(T)
        else
            x[end] = -temp / (2 * ε)
        end
    end

    return x
end
ChainRulesCore.@non_differentiable optimizeprimal!(f::MMADualObj, λ::AbstractVector)

"""
    getoptimalx(f::MMADualObj)

Returns the cached optimal primal solution.
"""
getoptimalx(f::MMADualObj) = f.x[1:getnvars(getparent(f.model))]

"""
    getapproxfg(dualobj::MMADualObj)

Returns the cached approximate objective and constraint values.
"""
getapproxfg(dualobj::MMADualObj) = getapproxfg(dualobj.model)

"""
```
struct MMADualModel <: AbstractModel
    parent::MMAApproxModel
    obj::MMADualObj
end
```

The dual model with a primal model `parent` and dual objective `obj`.
"""
@params struct MMADualModel <: AbstractModel
    parent::MMAApproxModel
    obj::MMADualObj
end

"""
    MMADualModel(model::MMAApproxModel)

Constructs an instance of `MMADualModel` for the approximate MMA model, `model`.
"""
function MMADualModel(model::MMAApproxModel)
    obj = MMADualObj(model)
    return MMADualModel(model, obj)
end

"""
    getoptimalx(m::MMADualModel)

Returns the last cached optimal primal solution.
"""
getoptimalx(m::MMADualModel) = getoptimalx(m.obj)

getparent(m::MMADualModel) = m.parent

@doc docσ
getσ(m::MMADualModel) = getσ(getparent(m))

@doc docσ
setσ!(m::MMADualModel, σ) = setσ!(getparent(m), σ)

@doc docρ
getρ(m::MMADualModel) = getρ(getparent(m))

@doc docρ
setρ!(m::MMADualModel, ρ) = setρ!(getparent(m), ρ)

@doc docupdateapprox!
function updateapprox!(m::MMADualModel, args...)
    updateapprox!(getparent(m), args...)
    return m
end

getobjective(m::MMADualModel) = m.obj

@doc docxk
getxk(m::MMADualModel) = getxk(getparent(m))

@doc docfk
getfk(m::MMADualModel) = getfk(getparent(m))

@doc doc∇fk
get∇fk(m::MMADualModel) = get∇fk(getparent(m))


text:
using dummyMindReader
using Documenter

DocMeta.setdocmeta!(dummyMindReader, :DocTestSetup, :(using dummyMindReader); recursive=true)

makedocs(;
    modules=[dummyMindReader],
    authors="DanielRivasMD <danielrivasmd@gmail.com> and contributors",
    repo="https://github.com/DanielRivasMD/dummyMindReader.jl/blob/{commit}{path}#{line}",
    sitename="dummyMindReader.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://DanielRivasMD.github.io/dummyMindReader.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/DanielRivasMD/dummyMindReader.jl",
)

max_stars_repo_name: DanielRivasMD/dummyMindReader

max_issues_repo_name: DanielRivasMD/dummyMindReader

max_forks_repo_name: DanielRivasMD/dummyMindReader


text:
#  This module replaces the JTAGManager module.
#  This implements the write function.
#  The write function is then used by the send function in the module
#  ImageSerializer.
#  This version uses Julia to create and manage Avalon Packets, instead
#  of the Altera C code.

module SPIManager

export SPI_GetNumChannels, SPI_OpenChannel, SPI_InitChannel, SPI_CloseChannel
export write, writeLEDs

using ProgressMeter
using Printf

include("AvalonPackets.jl")

#  Use ccall to use MPSSE C library functions.
#  Need to wrap all FTDI MPSSE functions in Julia functions.
#  Check for number of channels and report.  This will only work
#  if the FTDI device is plugged into USB and two kernal modules have
#  been removed:
#  sudo rmmod ftdi_sio usbserial
#  Note that this removal process is done by the udev system and a script.
#  Look in the directory /etc/udev/rules.d
#  The module removal is done a boot by a systemd service.

channel = Ref{Cuint}(0)  # pointer to 0

function SPI_GetNumChannels()
    ccall((:SPI_GetNumChannels,"libMPSSE"),UInt32, (Ref{Cuint},), channel)
#  The full path is required if the command ldconfig has not been used to add the path.
#ccall((:SPI_GetNumChannels,"/usr/local/lib/libMPSSE.so"),UInt32, (Ref{Cuint},), channel)
println("Number channels = $channel")
end

#  FTDI MPSSE struct which is used to configure the SPI channel.
struct ChannelConfig_t
    ClockRate::UInt32
    LatencyTimer::UInt8
    configOptions::UInt32
    Pin::UInt32
    reserved::UInt16
end

#  Configure the channel:
SPI_CONFIG_OPTION_MODE0= 0x00000000
SPI_CONFIG_OPTION_CS_DBUS3 =	0x00000000
SPI_CONFIG_OPTION_CS_ACTIVELOW =	0x00000020
configOptions = SPI_CONFIG_OPTION_MODE0 | SPI_CONFIG_OPTION_CS_DBUS3 | SPI_CONFIG_OPTION_CS_ACTIVELOW
Pin = 0x00000000
#  Initialize the struct.  First parameter is clock rate in Hz.  Maximum 30 MHz.
#  The 2nd parameter is the latency timer.  Set this value low for fast ReadWrite.
channelConf = ChannelConfig_t(30000000, 1, configOptions, Pin,0)

#  OPen and Initialize the SPI device.
function SPI_InitChannel()
handle = Ref{UInt32}()  #  This is written to by SPI_OpenChannel function.
opencall =  ccall((:SPI_OpenChannel, "/usr/local/lib/libMPSSE.so"), UInt32, (Cuint, Ref{UInt32}), 0, handle)
println("opencall = $opencall")
Printf.@printf("handle = %x\n", handle[])
initcall = ccall((:SPI_InitChannel, "/usr/local/lib/libMPSSE.so"), UInt32, (UInt32, Ref{ChannelConfig_t}), handle[], channelConf)
println("initcall = $initcall")
return handle[]
end

#  Close the SPI Channel and report.
function SPI_CloseChannel(handle)
spi_close = ccall((:SPI_CloseChannel, "/usr/local/lib/libMPSSE.so"), UInt32, (UInt32,), handle)
println("spiclose = $spi_close")
end

#  Need a "write" function which is similar to that in JTAGManager.
#  This function needs to safely open and close the SPI device channel.
#  This is done using the collection of FTDI MPSSE functions above.

#  This is the main Read/Write/ReadWrite wrapper.  Uncomment the 3 C function variants as required.
function spi_readwrite(handle, read_buf::Array{UInt8,1}, write_buf::Array{UInt8,1}, size)
   sizeTransfered = Ref{UInt32}(0)
   #ftdi_error = ccall((:SPI_Write, "/usr/local/lib/libMPSSE.so"), UInt32, (UInt32, Ptr{UInt8}, UInt32, Ref{UInt32}, UInt32), handle, write_buf, size, sizeTransfered, 0)
   #ftdi_error = @time ccall((:SPI_Read, "/usr/local/lib/libMPSSE.so"), UInt32, (UInt32, Ptr{UInt8}, UInt32, Ref{UInt32}, UInt32), handle, readbuffer, size, sizeTransfered, 0)
   ftdi_error = ccall((:SPI_ReadWrite, "/usr/local/lib/libMPSSE.so"), UInt32, (UInt32, Ptr{UInt8}, Ptr{UInt8}, UInt32, Ref{UInt32}, UInt32), handle, read_buf, write_buf, size, sizeTransfered, 0x00)
   if ftdi_error != 0
    println("FTDI Error!")
   end
   return ftdi_error
end

#  Transaction op codes:
SEQUENTIAL_WRITE = 0x04
SEQUENTIAL_READ  = 0x14
NON_SEQUENTIAL_WRITE = 0x00
NON_SEQUENTIAL_READ = 0x10

#  This is the primary write function which uses FTDI SPI device to write to FPGA Avalon bus.
function write(
            address,
            data
       )

    # Open the channel to the SPI device.
    SPI_GetNumChannels()
    thisHandle = SPI_InitChannel()

    # Compute the number of transfers
    #  This should be the number of individual frames in the animated GIF.
    bytes_per_transfer = 62000  #  Actual will be higher than this due to Avalon packet overhead.
    data_length = length(data)
    ntransfers = ceil(Int, data_length / bytes_per_transfer)
    println("The number ntransfers = $ntransfers")
    println("The image data length is $data_length")
    avalonbuffer2 = zeros(UInt8, 65536)
    readbuffer = zeros(UInt8, 65536)
    #  Alternate ways of creating arrays.
    #    readbuffer = Vector{UInt8}(undef, 66000)
    #    readbuffer = Base.Libc.malloc(166000)
    buffersize = 0
    transfer_data_length = 0

@showprogress for i in 1:ntransfers

        # Compute start address for this transfer
        start = address + bytes_per_transfer * (i-1) # Starts at address
        # Get the index range of data for this transfer
        idx_low = bytes_per_transfer * (i-1) + 1
        idx_high = min(bytes_per_transfer * i, length(data))

        #  This if statement is so the function can handle a single byte.
        if length(data) > 1  #  Required because view will not work on single byte.
        transfer_data = view(data, idx_low:idx_high)
        else
        transfer_data = data
        end

        transfer_data_length = length(transfer_data)

        # Transform to Avalon packet.  See include file AvalonPackets.jl.
        header, header_and_datalength = AvalonPackets.do_transaction(SEQUENTIAL_WRITE, transfer_data_length, start, transfer_data)
        send_data, send_data_length = AvalonPackets.tx_packet(header, header_and_datalength)
        avalonbuffer, buffersize = AvalonPackets.byte_to_core(send_data, send_data_length)

#  Write the data to the FTDI SPI device and thus to the FPGA:
spiwrite = spi_readwrite(thisHandle, readbuffer, avalonbuffer, buffersize)

        #  The data returned by byte_to_core to the avalonbuffer must be analyzed.
        #  It will include a response packet which contains the number of bytes written/read.
        #  The response packet is tacked onto the end.  So back up 12 indices from the end, and that should
        #  be the response packet.  There may be a few idle characters which can be ignored.
        #  Interpret the response.
response_packet = AvalonPackets.clean_response(readbuffer[(buffersize-12):buffersize])
#  This function returns a 16 bit unsigned integer, which should be the number of bytes written.
response_data = AvalonPackets.get_response_data(response_packet)

#  Verify the number of bytes written is what was intended.
if response_data != transfer_data_length
    println("Bad data transfer!!!")
    println("Response data = $response_data")
    println("i = $i")
    println("Buffersize = $buffersize")
    #  Look for possible problems caused by special characters.
    for k in 1:bytes_per_transfer
        if transfer_data[k] == 0x7a
            Printf.@printf("transfer_data[%d] = %x\n", k, transfer_data[k])
        elseif transfer_data[k] == 0x7b
            Printf.@printf("transfer_data[%d] = %x\n", k, transfer_data[k])
        elseif transfer_data[k] == 0x7c
            Printf.@printf("transfer_data[%d] = %x\n", k, transfer_data[k])
        elseif transfer_data[k] == 0x7d
            Printf.@printf("transfer_data[%d] = %x\n", k, transfer_data[k])
        elseif transfer_data[k] == 0x4a
            Printf.@printf("transfer_data[%d] = %x\n", k, transfer_data[k])
        elseif transfer_data[k] == 0x4d
            Printf.@printf("transfer_data[%d] = %x\n", k, transfer_data[k])
        end
    end
end

end # ntransfers
#    for k in (buffersize - 12):buffersize
#            Printf.@printf("readbuffer[%d] = %x\n", k, readbuffer[k])
#        end

#Printf.@printf("buffersize = %d\n", buffersize)

#println("Transfer data length = $transfer_data_length")
SPI_CloseChannel(thisHandle)
end

#  This function was used to debug the write function above.
function writeLEDs(
            address,
            data
       )

       avalonbuffer = zeros(UInt8, 100)
       readbuffer = zeros(UInt8, 100)

    # Open the channel to the SPI device.
    SPI_GetNumChannels()
    thisHandle = SPI_InitChannel()

    transfer_data_length = length(data)
    #    header = AvalonPackets.do_transaction(0x04, transfer_data_length, address, data)
    #    avalonbuffer = AvalonPackets.tx_packet(header)
    #    buffersize = length(avalonbuffer)

    header, header_and_datalength = AvalonPackets.do_transaction(SEQUENTIAL_WRITE, transfer_data_length, address, data)
    send_data, send_data_length = AvalonPackets.tx_packet(header, header_and_datalength)
    avalonbuffer, buffersize = AvalonPackets.byte_to_core(send_data, send_data_length)

        for i in 1:buffersize
        Printf.@printf("avalonbuffer[%d] = %x\n", i, avalonbuffer[i])
        end

        spiwrite = spi_readwrite(thisHandle, readbuffer, avalonbuffer, buffersize)
        println("spiwrite = $spiwrite")
        response_packet = AvalonPackets.clean_response(readbuffer[(buffersize-12):buffersize])
        response_data = AvalonPackets.get_response_data(response_packet)
        println("Response data = $response_data")

        for i in (buffersize - 12):buffersize
        Printf.@printf("readbuffer[%d] = %x\n", i, readbuffer[i])
        end

        for i in 1:length(response_packet)
        Printf.@printf("response_packet[%d] = %x\n", i, response_packet[i])
        end

        println("buffersize = $buffersize")
SPI_CloseChannel(thisHandle)
end

include("spi_read_avalon.jl")

end # module

max_stars_repo_name: Greg-R/spiavalonfpga

max_issues_repo_name: Greg-R/spiavalonfpga

max_forks_repo_name: Greg-R/spiavalonfpga


text:
# ---
# title: 1130. Minimum Cost Tree From Leaf Values
# id: problem1130
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming, Stack, Tree
# link: <https://leetcode.com/problems/minimum-cost-tree-from-leaf-values/description/>
# hidden: true
# ---
#
# Given an array `arr` of positive integers, consider all binary trees such
# that:
#
#   * Each node has either 0 or 2 children;
#   * The values of `arr` correspond to the values of each  **leaf** in an in-order traversal of the tree.  _(Recall that a node is a leaf if and only if it has 0 children.)_
#   * The value of each non-leaf node is equal to the product of the largest leaf value in its left and right subtree respectively.
#
# Among all possible binary trees considered, return the smallest possible sum
# of the values of each non-leaf node.  It is guaranteed this sum fits into a
# 32-bit integer.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [6,2,4]
#     Output: 32
#     Explanation:
#     There are two possible trees.  The first has non-leaf node sum 36, and the second has non-leaf node sum 32.
#
#         24            24
#        /  \          /  \
#       12   4        6    8
#      /  \               / \
#     6    2             2   4
#
#
#
#
# **Constraints:**
#
#   * `2 <= arr.length <= 40`
#   * `1 <= arr[i] <= 15`
#   * It is guaranteed that the answer fits into a 32-bit signed integer (ie. it is less than `2^31`).
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using SEML, Test

make_model("../example/testcase/case3tests.txt", model="FBA")
#sleep(1)
rm("autogeneratedmodel", recursive=true)

max_stars_repo_name: varnerlab/SEML

max_issues_repo_name: varnerlab/SEML

max_forks_repo_name: varnerlab/SEML


text:
using MyExample
using Documenter

DocMeta.setdocmeta!(MyExample, :DocTestSetup, :(using MyExample); recursive=true)

makedocs(;
    modules=[MyExample],
    authors="Edouard <edouard.oudet@imag.fr> and contributors",
    repo="https://github.com/edljk/MyExample.jl/blob/{commit}{path}#{line}",
    sitename="MyExample.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://edljk.github.io/MyExample.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/edljk/MyExample.jl",
    devbranch="main",
)