text:
using .Locales
using .LC

# Construct LocaleId from language tags
@test LocaleId("") === DEFAULT
@test LocaleId("C") === ROOT
@test LocaleId("en") === ENGLISH
@test LocaleId("en-US") === US
@test_throws ArgumentError("missing language prefix") LocaleId("c")
@test string(LocaleId("de")) == "de"
@test string(LocaleId("de_latn")) == "de_Latn"
@test string(LocaleId("de-de")) == "de_DE"
@test string(LocaleId("de-111")) == "de_111"
@test string(LocaleId("de-variAnt1")) == "de_variant1"
@test string(LocaleId("de-variAnt1_varia2")) == "de_variant1_varia2"
@test string(LocaleId("de-1erw-varia2")) == "de_1erw_varia2"
@test string(LocaleId("de-a_aB")) == "de_a_ab"
@test string(LocaleId("DE_X_1_UZ_A_B")) == "de_x_1_uz_a_b"
@test string(LocaleId("DE_latn_de")) == "de_Latn_DE"
@test string(LocaleId("DE_latn_de-variant1")) == "de_Latn_DE_variant1"
@test string(LocaleId("DE_latn_de_a-uu-11")) == "de_Latn_DE_a_uu_11"
@test string(LocaleId("DE_latn_de-varia-a-uu-11")) == "de_Latn_DE_varia_a_uu_11"
@test string(LocaleId("DE_latn_de-varia-a-uu-11")) == "de_Latn_DE_varia_a_uu_11"
@test string(LocaleId("DE_latn_de-1var-a-uu-11")) == "de_Latn_DE_1var_a_uu_11"
@test string(LocaleId("DE_latn_de-1va2-a-uu-11")) == "de_Latn_DE_1va2_a_uu_11"
@test string(LocaleId("DE_latn_de-x-11-a-uu-11")) == "de_Latn_DE_x_11_a_uu_11"
@test string(LocaleId("DE_latn_de-b-22-a-uu-11")) == "de_Latn_DE_a_uu_11_b_22"
@test string(LocaleId("DE_latn_de-z-22-a-uu-11")) == "de_Latn_DE_a_uu_11_z_22"

@test string(LocaleId("zh-han")) == "han"
@test string(LocaleId("deu")) == "de"

@test LocaleId("en") ⊆ LocaleId("C")
@test BOTTOM ⊆ LocaleId("en-Latn-ab-valencia-a-bc-x-1-2")
@test LocaleId("en-a-aa-b-bb") ⊆ LocaleId("en-b-bb")
@test hash(LocaleId("is-a-aa-b-bb")) == hash(LocaleId("is-b-bb-a-aa"))
@test LocaleId("is-a-aa-b-bb") === LocaleId("is-b-bb-a-aa")
@test LocaleId("C") != LocaleId("")

# exceptions
create2(l="", lex=String[], s="", r="", v=String[], d=Dict{Char,Vector{Symbol}}()) =
    create_locid(l, lex, s, r, v, d)

@test_throws ArgumentError("missing language prefix") LocaleId("a-ab")
@test_throws ArgumentError("no language prefix 'a'") LocaleId("a", "")
@test_throws ArgumentError("only one language extension allowed 'de_abc_def'") LocaleId("de-abc-def")
@test_throws ArgumentError("no language exensions allowed 'abcd_def'") create2("abcd", ["def"])
@test_throws ArgumentError("no script 'Abc'") LocaleId("ab", "abc", "de")
@test_throws ArgumentError("no region 'DEF'") LocaleId("ab", "abcd", "def")
@test_throws ArgumentError("no variants 'vari'") create2("ab", String[], "", "", ["vari"])
@test_throws ArgumentError("no variants '1va'") create2("ab", String[], "", "", ["1va"])
@test_throws ArgumentError("no variants 'asdfghjkl'") create2("ab", String[], "", "", ["asdfghjkl"])
@test_throws ArgumentError("no variants '1990_asdfghjkl'") create2("ab", String[], "", "", ["1990", "asdfghjkl"])
@test_throws ArgumentError("language tag contains invalid characters: 'Ä'") LocaleId("Ä")
@test_throws ArgumentError("multiple occurrence of singleton 'u'") LocaleId("en-u-u1-u-u2")
@test_throws ArgumentError("no language tag: 'en_asdfghjkl' after 1") LocaleId("en_asdfghjkl")
@test_throws ArgumentError("no language tag: 'asdfghjkl' after 0") LocaleId("asdfghjkl")
@test_throws ArgumentError("no language tag: 'x_asdfghjkl' after 1") LocaleId("x-asdfghjkl")
@test_throws ArgumentError("too many language extensions 'en_abc_def_ghi_jkl'") LocaleId("en-abc-def-ghi-jkl")

#various constructors
@test LocaleId("") === DEFAULT
@test LocaleId("C") === ROOT
@test LocaleId() === BOTTOM
@test LocaleId("de", "de") == LocaleId("de-de")
@test LocaleId("de", "latn", "de") == LocaleId("de-Latn-de")
@test LocaleId("de", "latn", "de", "bavarian") == LocaleId("de-Latn-de-bavarian")
@test LocaleId("de", "de") == LocaleId("de-de")
@test LocaleId("de", "de") == LocaleId("de-de")
@test LocaleId("de", "de") == LocaleId("de-de")
@test LocaleId("de", "de") == LocaleId("de-de")

# predefined locales
@test ENGLISH === LocaleId("en")
@test FRENCH === LocaleId("fr")
@test GERMAN === LocaleId("de")
@test ITALIAN === LocaleId("it")
@test JAPANESE === LocaleId("ja")
@test KOREAN === LocaleId("ko")
@test CHINESE === LocaleId("zh")
@test SIMPLIFIED_CHINESE === LocaleId("zh_CN")
@test TRADITIONAL_CHINESE === LocaleId("zh_TW")
@test FRANCE === LocaleId("fr_FR")
@test GERMANY === LocaleId("de_DE")
@test ITALY === LocaleId("it_IT")
@test JAPAN === LocaleId("ja_JP")
@test KOREA === LocaleId("ko_KR")
@test CHINA === LocaleId("zh_CN")
@test PRC === LocaleId("zh_CN")
@test TAIWAN === LocaleId("zh_TW")
@test UK === LocaleId("en_GB")
@test US === LocaleId("en_US")
@test CANADA === LocaleId("en_CA")

# inclusion and ordering
@test LocaleId() ⊆ LocaleId("en-Latf-gb-variant-variant2-a-123-26-x-1-2-23a-bv")
@test LocaleId("en-Latf-gb-variant-variant2-a-123-26-x-1-2-23a-bv") ⊆ LocaleId("")
@test !(LocaleId("fr-Latn") ⊆ LocaleId("fr-CA"))
@test !(LocaleId("fr-CA") ⊆ LocaleId("fr-Latn"))
@test !(LocaleId("fr-Latn") < LocaleId("fr-CA"))
@test LocaleId("fr-Latn") ≥ LocaleId("fr-CA")

# loading initial values from environment variables
delete!(ENV, "LANG")
delete!(ENV, "LC_MESSAGES")
delete!(ENV, "LC_COLLATE")
delete!(ENV, "LC_TIME")
delete!(ENV, "LC_NUMERIC")
delete!(ENV, "LC_MONETARY")

ENV["LC_ALL"] = "en_GB"
@test default_locale(LC.MESSAGES) === LocaleId("en-GB")
ENV["LC_ALL"] = "en_GB.utf8"
@test default_locale(LC.COLLATE) === LocaleId("en-GB")
ENV["LC_ALL"] = "en_GB.utf8@oed"
@test default_locale(LC.TIME) === LocaleId("en-GB-x-posix-oed")
ENV["LC_ALL"] = "en_GB@oed"
@test default_locale(LC.NUMERIC) === LocaleId("en-GB-x-posix-oed")

delete!(ENV, "LC_ALL")
ENV["LC_MONETARY"] = "en_US.UTF-8"
@test default_locale(LC.MONETARY) === LocaleId("en-US")
ENV["LC_TIME"] = "en_CA"
@test default_locale(LC.TIME) === LocaleId("en-ca")

ENV["LC_TIME"] = ""
delete!(ENV, "LC_MONETARY")
ENV["LANG"] = "fr_FR@guadelo"
@test default_locale(LC.TIME) === LocaleId("fr-FR-x-posix-guadelo")
ENV["LANG"] = "C"
@test default_locale(LC.TIME) === LocaleId("C")


text:
"""
$(SIGNATURES)
Solve the dynamic programming problem, starting from the terminal time step
Computes the gain matrices K and d by applying the principle of optimality at
each time step, solving for the gradient (s) and Hessian (S) of the cost-to-go
function. Also returns parameters Δv for line search (see Synthesis and Stabilization of Complex Behaviors through
Online Trajectory Optimization)
"""
function backwardpass!(results::SolverVectorResults,solver::Solver)
    if solver.opts.square_root
        Δv = _backwardpass_sqrt!(results, solver)
    else
        Δv = _backwardpass!(results, solver)
    end
    return Δv
end

function _backwardpass!(res::SolverVectorResults,solver::Solver)
    # Get problem sizes
    n,m,N = get_sizes(solver)
    m̄,mm = get_num_controls(solver)
    n̄,nn = get_num_states(solver)

    # Objective
    costfun = solver.obj.cost

    # Minimum time and infeasible options
    solver.state.minimum_time ? R_minimum_time = solver.opts.R_minimum_time : nothing
    solver.state.infeasible ? R_infeasible = solver.opts.R_infeasible*Matrix(I,n,n) : nothing

    dt = solver.dt

    X = res.X; U = res.U; K = res.K; d = res.d; S = res.S; s = res.s

    Qx = res.bp.Qx; Qu = res.bp.Qu; Qxx = res.bp.Qxx; Quu = res.bp.Quu; Qux = res.bp.Qux
    Quu_reg = res.bp.Quu_reg; Qux_reg = res.bp.Qux_reg

    # TEMP resets values for now - this will get fixed
    for k = 1:N-1
        Qx[k] = zeros(nn); Qu[k] = zeros(mm); Qxx[k] = zeros(nn,nn); Quu[k] = zeros(mm,mm); Qux[k] = zeros(mm,nn)
        Quu_reg[k] = zeros(mm,mm); Qux_reg[k] = zeros(mm,nn)
    end

    # Boundary Conditions
    S[N], s[N] = quaternion_expansion(costfun, X[N])

    if solver.state.minimum_time
        s[N][n̄] = R_minimum_time*X[N][n̄]
        S[N][n̄,n̄] = R_minimum_time
    end

    # Initialize expected change in cost-to-go
    Δv = zeros(2)

    # Terminal constraints
    if res isa ConstrainedIterResults
        C = res.C; Iμ = res.Iμ; λ = res.λ
        Cx = res.Cx; Cu = res.Cu

        S[N] += Cx[N]'*Iμ[N]*Cx[N]
        s[N] += Cx[N]'*(Iμ[N]*C[N] + λ[N])
    end

    # Backward pass
    k = N-1
    while k >= 1
        solver.state.minimum_time ? dt = U[k][m̄]^2 : dt = solver.dt

        x = X[k][1:n]
        u = U[k][1:m]

        if solver.opts.sat_att == true
            expansion = quaternion_expansion(costfun,x,u)
            Qxx[k],Quu[k],Qux[k],Qx[k],Qu[k] = expansion .* dt
        else
            expansion = taylor_expansion(costfun,x,u)
            Qxx[k][1:n,1:n],Quu[k][1:m,1:m],Qux[k][1:m,1:n],Qx[k][1:n],Qu[k][1:m] = expansion .* dt
        end

        # Minimum time expansion components
        if solver.state.minimum_time
            ℓ1 = stage_cost(costfun,x,u)
            h = U[k][m̄]
            tmp = 2*h*expansion[5]

            Qu[k][m̄] = h*(2*ℓ1 + R_minimum_time)
            Quu[k][1:m,m̄] = tmp
            Quu[k][m̄,1:m] = tmp'
            Quu[k][m̄,m̄] = (2*ℓ1 + R_minimum_time)
            Qux[k][m̄,1:n] = 2*h*expansion[4]'

            Qx[k][n̄] = R_minimum_time*X[k][n̄]
            Qxx[k][n̄,n̄] = R_minimum_time
        end

        # Infeasible expansion components
        if solver.state.infeasible
            Qu[k][m̄+1:mm] = R_infeasible*U[k][m̄+1:m̄+n]
            Quu[k][m̄+1:mm,m̄+1:mm] = R_infeasible
        end

        # Compute gradients of the dynamics
        fdx, fdu = res.fdx[k], res.fdu[k]

        if solver.opts.sat_att == true
            #renormalize for quaternion
            qk = X[k]
            sk = qk[4]
            vk = qk[5:7]

            qn = X[k+1]
            sn = qn[4]
            vn = qn[5:7]

            Gk = [-vk'; sk*Array(1.0*Diagonal(I,3)) + hat(vk)]
            Gn = [-vn'; sn*Array(1.0*Diagonal(I,3)) + hat(vn)]

            #recompute A and B matrices through permutation matrices
            perm_Gn = zeros(7,8)
            perm_Gk = zeros(8,7)
            perm_Gn[1:3,1:3]= Array(1.0*Diagonal(I,3))
            perm_Gn[4:6,4:7] = Gn'
            perm_Gk[1:3,1:3]= Array(1.0*Diagonal(I,3))
            perm_Gk[4:7,4:6] = Gk

            fdx = perm_Gn*fdx*perm_Gk
            fdu = perm_Gn*fdu
        end

        # Gradients and Hessians of Taylor Series Expansion of Q
        Qx[k] += fdx'*s[k+1]
        Qu[k] += fdu'*s[k+1]
        Qxx[k] += fdx'*S[k+1]*fdx
        Quu[k] += fdu'*S[k+1]*fdu
        Qux[k] += fdu'*S[k+1]*fdx

        # # Constraints
        # if res isa ConstrainedIterResults
        #     Qx[k] += Cx[k]'*(Iμ[k]*C[k] + λ[k])
        #     Qu[k] += Cu[k]'*(Iμ[k]*C[k] + λ[k])
        #     Qxx[k] += Cx[k]'*Iμ[k]*Cx[k]
        #     Quu[k] += Cu[k]'*Iμ[k]*Cu[k]
        #     Qux[k] += Cu[k]'*Iμ[k]*Cx[k]
        # end

        if solver.opts.bp_reg_type == :state
            Quu_reg[k] = Quu[k] + res.ρ[1]*fdu'*fdu
            Qux_reg[k] = Qux[k] + res.ρ[1]*fdu'*fdx
        elseif solver.opts.bp_reg_type == :control
            Quu_reg[k] = Quu[k] + res.ρ[1]*I
            Qux_reg[k] = Qux[k]
        end

        # Regularization
        if !isposdef(Hermitian(Array(Quu_reg[k])))  # need to wrap Array since isposdef doesn't work for static arrays
            # increase regularization
            @logmsg InnerIters "Regularizing Quu "
            regularization_update!(res,solver,:increase)

            # reset backward pass
            k = N-1
            Δv[1] = 0.
            Δv[2] = 0.
            continue
        end

        # Compute gains
        K[k] = -Quu_reg[k]\Qux_reg[k]
        d[k] = -Quu_reg[k]\Qu[k]

        # Calculate cost-to-go (using unregularized Quu and Qux)
        s[k] = (Qx[k] + K[k]'*Quu[k]*d[k] + K[k]'*Qu[k] + Qux[k]'*d[k])
        S[k] = (Qxx[k] + K[k]'*Quu[k]*K[k] +
         K[k]'*Qux[k] + Qux[k]'*K[k])
        S[k] = 0.5*(S[k] + S[k]')

        # calculated change is cost-to-go over entire trajectory
        Δv[1] += d[k]'*Qu[k]
        Δv[2] += 0.5*d[k]'*Quu[k]*d[k]

        k = k - 1;
    end

    # decrease regularization after backward pass
    regularization_update!(res,solver,:decrease)

    return Δv
end

function al_expansion!(res::SolverIterResults, solver::Solver, k::Int)
    X = res.X; U = res.U
    costfun = solver.obj.cost
    n,m,N = get_sizes(solver)
    m̄,mm = get_num_controls(solver)
    n̄,nn = get_num_states(solver)
    R_infeasible = solver.opts.R_infeasible

    Qx = res.bp.Qx; Qu = res.bp.Qu; Qxx = res.bp.Qxx; Quu = res.bp.Quu; Qux = res.bp.Qux

    x = X[k][1:n]

    if k == N
        Qxx[k][1:n,1:n],Qx[k][1:n] = taylor_expansion(costfun,x)
    else
        u = U[k][1:m]
        expansion = taylor_expansion(costfun,x,u)
        Qxx[k][1:n,1:n],Quu[k][1:m,1:m],Qux[k][1:m,1:n],Qx[k][1:n],Qu[k][1:m] = expansion .* solver.dt
    end

    if solver.state.constrained
        C = res.C; Iμ = res.Iμ; λ = res.λ; μ = res.μ
        Cx = res.Cx; Cu = res.Cu

        if solver.opts.al_type == :default
            Qx[k] += Cx[k]'*(Iμ[k]*C[k] + λ[k])
            Qxx[k] += Cx[k]'*Iμ[k]*Cx[k]
            if k < solver.N
                Qu[k] += Cu[k]'*(Iμ[k]*C[k] + λ[k])
                Quu[k] += Cu[k]'*Iμ[k]*Cu[k]
                Qux[k] += Cu[k]'*Iμ[k]*Cx[k]
            end
        elseif solver.opts.al_type == :algencan
            a = res.active_set[k]
            Ia = Diagonal(a)
            Iμk = Diagonal(μ[k])

            # Declare some temp variables that are re-used
            gg = Ia*Iμk
            g = Ia*(λ[k] + Iμk*C[k])

            # Expansion
            Qx[k] += Cx[k]'g
            Qxx[k] += Cx[k]'gg*Cx[k]
            if k < solver.N
                Qu[k] += Cu[k]'g
                Quu[k] += Cu[k]'gg*Cu[k]
                Qux[k] += Cu[k]'gg*Cx[k]
            end

            # Qx[k] += Cx[k]'*(Iμ[k]*C[k] + λ[k])
            # Qxx[k] += Cx[k]'*Iμ[k]*Cx[k]
            # if k < solver.N
            #     Qu[k] += Cu[k]'*(Iμ[k]*C[k] + λ[k])
            #     Quu[k] += Cu[k]'*Iμ[k]*Cu[k]
            #     Qux[k] += Cu[k]'*Iμ[k]*Cx[k]
            # end
        end

        # Minimum time expansion components
        if solver.state.minimum_time
            R_minimum_time = solver.opts.R_minimum_time
            Qx[k][n̄] = R_minimum_time*X[k][n̄]
            Qxx[k][n̄,n̄] = R_minimum_time
            if k < N
                ℓ1 = stage_cost(costfun,x,u)
                h = U[k][m̄]
                tmp = 2*h*expansion[5]

                Qu[k][m̄] = h*(2*ℓ1 + R_minimum_time)
                Quu[k][1:m,m̄] = tmp
                Quu[k][m̄,1:m] = tmp'
                Quu[k][m̄,m̄] = (2*ℓ1 + R_minimum_time)
                Qux[k][m̄,1:n] = 2*h*expansion[4]'
            end
        end

        # Infeasible expansion components
        if solver.state.infeasible && k < N
            Qu[k][m̄+1:mm] = R_infeasible*U[k][m̄+1:m̄+n]
            Quu[k][m̄+1:mm,m̄+1:mm] = R_infeasible*Diagonal(I,n)
        end
    end
end

"""
$(SIGNATURES)
Perform a backwards pass with Cholesky Factorizations of the Cost-to-Go to
avoid ill-conditioning.
"""
function _backwardpass_sqrt!(res::SolverVectorResults,solver::Solver)
    # Get problem sizes
    n,m,N = get_sizes(solver)
    m̄,mm = get_num_controls(solver)
    n̄,nn = get_num_states(solver)

    # Objective
    costfun = solver.obj.cost

    # Minimum time and infeasible options
    solver.state.minimum_time ? R_minimum_time = solver.opts.R_minimum_time : nothing
    solver.state.infeasible ? R_infeasible = solver.opts.R_infeasible*Matrix(I,n,n) : nothing

    dt = solver.dt

    X = res.X; U = res.U; K = res.K; d = res.d; Su = res.S; s = res.s

    Qx = res.bp.Qx; Qu = res.bp.Qu; Qxx = res.bp.Qxx; Quu = res.bp.Quu; Qux = res.bp.Qux
    Quu_reg = res.bp.Quu_reg; Qux_reg = res.bp.Qux_reg

    # TEMP resets values for now - this will get fixed
    for k = 1:N-1
        Qx[k] = zeros(nn); Qu[k] = zeros(mm); Qxx[k] = zeros(nn,nn); Quu[k] = zeros(mm,mm); Qux[k] = zeros(mm,nn)
        Quu_reg[k] = zeros(mm,mm); Qux_reg[k] = zeros(mm,nn)
    end

    # Boundary Conditions
    Su[N][1:n,1:n], s[N][1:n] = taylor_expansion(costfun, X[N][1:n])

    if solver.state.minimum_time
        s[N][n̄] = R_minimum_time*X[N][n̄]
        Su[N][n̄,n̄] = R_minimum_time
    end

    # Take square root (via cholesky)
    try
        Su[N][1:nn,1:nn] = cholesky(Su[N][1:nn,1:nn]).U # if no terminal cost is provided cholesky will fail gracefully
    catch PosDefException
        if sum([Su[N][i,i] for i = 1:n]) != 0. # TODO there may be something faster here, but lu fails for this case
            # tmp = svd(Su[N])
            # Su[N] = Diagonal(sqrt.(tmp.S))*tmp.V'
            tmp = eigen(Su[N])
            Su[N] = Diagonal(sqrt.(tmp.values))*tmp.vectors'
        elseif tr(Su[N][1:n,1:n]) == 0. && n̄ > n
            Su[N][n̄,n̄] = sqrt(Su[N][n̄,n̄])
        end
    end

    # Terminal constraints
    if res isa ConstrainedIterResults
        C = res.C; Iμ = res.Iμ; λ = res.λ
        Cx = res.Cx; Cu = res.Cu
        Iμ_sqrt = sqrt.(Iμ[N])

        Su[N][1:nn,1:nn] = chol_plus(Su[N],Iμ_sqrt*Cx[N])
        s[N] += Cx[N]'*(Iμ[N]*C[N] + λ[N])
    end

    # Backward pass
    Δv = zeros(2)
    tmp1 = []
    tmp2 = []

    k = N-1
    while k >= 1
        solver.state.minimum_time ? dt = U[k][m̄]^2 : dt = solver.dt

        x = X[k][1:n]
        u = U[k][1:m]

        expansion = taylor_expansion(costfun,x,u)
        Qxx[k][1:n,1:n],Quu[k][1:m,1:m],Qux[k][1:m,1:n],Qx[k][1:n],Qu[k][1:m] = expansion .* dt

        # Minimum time expansion components
        if solver.state.minimum_time
            ℓ1 = stage_cost(costfun,x,u)
            h = U[k][m̄]
            tmp = 2*h*expansion[5]

            Qu[k][m̄] = h*(2*ℓ1 + R_minimum_time)
            Quu[k][1:m,m̄] = tmp
            Quu[k][m̄,1:m] = tmp'
            Quu[k][m̄,m̄] = (2*ℓ1 + R_minimum_time)
            Qux[k][m̄,1:n] = 2*h*expansion[4]'

            Qx[k][n̄] = R_minimum_time*X[k][n̄]
            Qxx[k][n̄,n̄] = R_minimum_time
        end

        # Infeasible expansion components
        if solver.state.infeasible
            Qu[k][m̄+1:mm] = R_infeasible*U[k][m̄+1:m̄+n]
            Quu[k][m̄+1:mm,m̄+1:mm] = R_infeasible
        end

        # Compute gradients of the dynamics
        fdx, fdu = res.fdx[k], res.fdu[k]

        # Gradients and Hessians of Taylor Series Expansion of Q
        Qx[k] += fdx'*s[k+1]
        Qu[k] += fdu'*s[k+1]
        Qux[k] += (fdu'*Su[k+1]')*(Su[k+1]*fdx)
        try
            Qxx[k][1:nn,1:nn] = cholesky(Qxx[k][1:nn,1:nn]).U
        catch
            if sum([Qxx[k][i,i] for i = 1:n]) != 0. #TODO faster
                tmp = eigen(Qxx[k])
                Qxx[k] = Diagonal(sqrt.(tmp.values))*tmp.vectors'
            elseif tr(Qxx[k][1:n,1:n]) == 0. && n̄ > n
                Qxx[k][n̄,n̄] = sqrt(Qxx[k][n̄,n̄])
            end
        end
        try
            Quu[k] = cholesky(Quu[k]).U
        catch #TODO fix this...
            error("problem with sqrt bp Quu...")
            # tmp = svd(Quu[k])
            # Quu[k] = Diagonal(sqrt.(tmp.S))*tmp.V'

            # tmp = eigen(Quu[k] + fdu'*Su[k+1]'*Su[k+1]*fdu)
            # Wuu = Diagonal(sqrt.(tmp.values))*tmp.vectors'
        end
        Wxx = chol_plus(Qxx[k], Su[k+1]*fdx)
        Wuu = chol_plus(Quu[k], Su[k+1]*fdu)

        # Constraints
        if res isa ConstrainedIterResults
            Iμ_sqrt = sqrt.(Iμ[k])

            Qx[k] += Cx[k]'*(Iμ[k]*C[k] + λ[k])
            Qu[k] += Cu[k]'*(Iμ[k]*C[k] + λ[k])
            Wxx = chol_plus(Wxx,Iμ_sqrt*Cx[k])
            Wuu = chol_plus(Wuu,Iμ_sqrt*Cu[k])
            Qux[k] += Cu[k]'*Iμ[k]*Cx[k]
        end

        if solver.opts.bp_reg_type == :state
            Wuu_reg = chol_plus(Wuu,sqrt(res.ρ[1])*fdu)
            Qux_reg[k] = Qux[k] + res.ρ[1]*fdu'*fdx
        elseif solver.opts.bp_reg_type == :control
            Wuu_reg = chol_plus(Wuu,sqrt(res.ρ[1])*Matrix(I,mm,mm))
            Qux_reg[k] = Qux[k]
        end

        #TODO find better PD check for Wuu_reg
        # # Regularization
        # if !isposdef(Hermitian(Array(Wuu_reg)))  # need to wrap Array since isposdef doesn't work for static arrays
        #     # increase regularization
        #     regularization_update!(res,solver,:increase)
        #
        #     # reset backward pass
        #     k = N-1
        #     Δv[1] = 0.
        #     Δv[2] = 0.
        #     continue
        # end

        # Compute gains
        K[k] = -Wuu_reg\(Wuu_reg'\Qux_reg[k])
        d[k] = -Wuu_reg\(Wuu_reg'\Qu[k])

        # Calculate cost-to-go
        s[k] = Qx[k] + (K[k]'*Wuu')*(Wuu*d[k]) + K[k]'*Qu[k] + Qux[k]'*d[k]

        try
            tmp1 = (Wxx')\Qux[k]'
        catch SingularException
            # if solver.opts.bp_sqrt_inv_type == :reg
            #     reg = solver.opts.bp_reg_sqrt_initial
            #     while minimum(eigvals(Wxx)) < 1e-6
            #         Wxx += reg*Matrix(I,nn,nn)
            #         try
            #             tmp1 = (Wxx')\Qux[k]'
            #             break
            #         catch
            #             reg *= solver.opts.bp_reg_sqrt_increase_factor
            #             if reg >= solver.opts.bp_reg_max
            #                 error("Square root regularization exceded")
            #             end
            #         end
            #     end
            # elseif solver.opts.bp_sqrt_inv_type == :pseudo
            tmp1 = pinv(Array(Wxx'))*Qux[k]'
            # end
        end

        try
            tmp2 = cholesky(Wuu'*Wuu - tmp1'*tmp1).U
        catch
            tmp = eigen(Wuu'*Wuu - tmp1'*tmp1)
            tmp2 = Diagonal(sqrt.(tmp.values))*tmp.vectors'
            # tmp = svd(Wuu'*Wuu - tmp1'*tmp1)
            # tmp2 = Diagonal(sqrt.(tmp.S))*tmp.V'
        end

        Su[k][1:nn,1:nn] = Wxx + tmp1*K[k]
        Su[k][nn+1:nn+mm,1:nn] = tmp2*K[k]

        # calculated change is cost-to-go over entire trajectory
        Δv[1] += d[k]'*Qu[k]
        Δv[2] += 0.5*d[k]'*Wuu'*Wuu*d[k]

        Quu_reg[k] = Array(Wuu_reg)

        Quu[k] = Array(Wuu)
        Qxx[k] = Array(Wxx)

        k = k - 1;
    end

    # decrease regularization after backward pass
    regularization_update!(res,solver,:decrease)

    return Δv
end

function chol_plus(A,B)
    n1,m = size(A)
    n2 = size(B,1)
    P = zeros(n1+n2,m)
    P[1:n1,:] = A
    P[n1+1:end,:] = B
    return qr(P).R
end

function backwardpass_max_condition_number(bp::TrajectoryOptimization.BackwardPass)
    N = length(bp.Quu)
    max_cn = 0.
    for k = 1:N-1
        cn = cond(bp.Quu_reg[k])
        if cn > max_cn
            max_cn = cn
        end
    end
    return max_cn
end

function backwardpass_max_condition_number(results::TrajectoryOptimization.SolverVectorResults)
    N = length(results.S)
    max_cn = 0.
    for k = 1:N
        cn = cond(results.S[k])
        if cn > max_cn && cn < Inf
            max_cn = cn
        end
    end
    return max_cn
end

function L_decomp(q)
    #this function takes in a desired quaternion computes the L decomposition
    L = [q[1] -q[2] -q[3] -q[4];
         -q[2] q[1] q[4] -q[3];
         -q[3] -q[4] q[1] q[2];
         -q[4] q[3] -q[2] q[1]]

         return L

end

ext: jl


text:
module NVRTC

using ..CUJulia

include("../lib/$(CUJulia.major)$(CUJulia.minor)/libnvrtc.jl")
include("../lib/$(CUJulia.major)$(CUJulia.minor)/libnvrtc_types.jl")

if is_windows()
    const libnvrtc = Libdl.find_library(["nvrtc64_80","nvrtc64_75"])
else
    const libnvrtc = Libdl.find_library(["libnvrtc"])
end
isempty(libnvrtc) && error("NVRTC library cannot be found.")

function check_nvrtcresult(status)
    status == NVRTC_SUCCESS && return
    throw(unsafe_string(nvrtcGetErrorString(status)))
end

function getlog(prog::Ptr{Void})
    logsize = Csize_t[0]
    nvrtcGetProgramLogSize(prog, logsize)
    log = Array(UInt8, logsize[1])
    nvrtcGetProgramLog(prog, log)
    unsafe_string(pointer(log))
end

function compile(code::String; headers=(), include_names=())
    p = Ptr{Void}[0]
    headers = Ptr{Void}[pointer(h) for h in headers]
    include_names = Ptr{Void}[pointer(n) for n in include_names]
    nvrtcCreateProgram(p, code, C_NULL, length(headers), headers, include_names)
    prog = p[1]
    options = ["--gpu-architecture=compute_30"]
    try
        nvrtcCompileProgram(prog, length(options), options)
    catch
        throw(getlog(prog))
    end
    ptxsize = Csize_t[0]
    nvrtcGetPTXSize(prog, ptxsize)
    ptx = Array(UInt8, ptxsize[1])
    nvrtcGetPTX(prog, ptx)
    nvrtcDestroyProgram(p)
    String(ptx)
end

end


text:
export_plot(z, all_data, ϕ_all, filename, ylabel) = nothing
export_plot_snapshot(z, all_data, ϕ_all, filename, ylabel) = nothing

# using Requires
# @init @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
#   using .Plots

"""
    plot_friendly_name(ϕ)

Get plot-friendly string, since many Unicode
characters do not render in plot labels.
"""
function plot_friendly_name(ϕ)
    s = ϕ
    s = replace(s, "ρ" => "rho")
    s = replace(s, "α" => "alpha")
    s = replace(s, "∂" => "partial")
    s = replace(s, "∇" => "nabla")
    return s
end

"""
    export_plot(z, all_data, ϕ_all, filename, ylabel)

Export plot of all variables, or all
available time-steps in `all_data`.
"""
function export_plot(z, all_data, ϕ_all, filename, ylabel)
    ϕ_all isa Tuple || (ϕ_all = (ϕ_all,))
    p = plot()
    for n in 0:(length(keys(all_data)) - 1)
        for ϕ in ϕ_all
            ϕ_string = String(ϕ)
            ϕ_name = plot_friendly_name(ϕ_string)
            ϕ_data = all_data[n][ϕ_string][:]
            plot!(ϕ_data, z, xlabel = ϕ_name, ylabel = ylabel)
        end
    end
    savefig(filename)
end

"""
    export_plot_snapshot(z, all_data, ϕ_all, filename, ylabel)

Export plot of all variables in `all_data`
"""
function export_plot_snapshot(z, all_data, ϕ_all, filename, ylabel)
    ϕ_all isa Tuple || (ϕ_all = (ϕ_all,))
    p = plot()
    for ϕ in ϕ_all
        ϕ_string = String(ϕ)
        ϕ_name = plot_friendly_name(ϕ_string)
        ϕ_data = all_data[ϕ_string][:]
        plot!(ϕ_data, z, xlabel = ϕ_name, ylabel = ylabel)
    end
    savefig(filename)
end

# end


text:
module FTCTests

using Reexport
@reexport using FaultTolerantControl
const FTC = FaultTolerantControl
using UnPack
using Plots
using Transducers
using LinearAlgebra
using JLD2, FileIO
using Printf
using NumericalIntegration
using ReferenceFrameRotations
using StaticArrays: SMatrix
using DataFrames
using Random
using DifferentialEquations

# privileged name, see #7
TRAJ_DATA_NAME = "traj.jld2"

# simulation
export run_sim, run_multiple_sim
# utils
export evaluate, extract_fault_property, get_fault_kind
# data manipulation
export partitionTrainTest

include("data_manipulation.jl")
include("run_sim.jl")
include("run_multiple_sim.jl")
include("evaluate.jl")
include("fault_info.jl")

end


text:
using DynamicTimeWarp
using Test

@test DynamicTimeWarp.trackback([
     0    0    0    0  16   65  146  227  276  292  292
    16   16   16   16   0    9   34   59   68   68   84
    65   65   65   65   9    0    4    8    8   17   66
   146  146  146  146  34    4    0    0    4   29   98
   195  195  195  195  43    4    4    4    0    9   58
   211  211  211  211  43   13   29   29    9    0   16
   211  211  211  211  59   62   94  110   58   16    0
   211  211  211  211  75  108  143  175  107   32    0
]) == (
   [1, 1, 1, 1, 2, 3, 4, 4, 5, 6, 7, 8],
   [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 11]
)

@test dtw(
   [0,4,7,9,10,9,7,4,0,0],
   [0,0,0,0,4,7,9,10,10,9,7,4,0]
) == (
   [1, 1, 1, 1, 2, 3, 4, 5, 5, 6, 7, 8, 9, 10],
   [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 13],
   0.0
)

@test linkseq([0,4,7,9,10,9,7,4,0,0], [0,0,0,0,4,7,9,10,10,9,7,4,0]) == (
   [[0.0 0.0]; [0.0 1.0]; [0.0 2.0]; [0.0 3.0];
   [1.0 4.0]; [2.0 5.0]; [3.0 6.0]; [4.0 7.0];
   [4.0 8.0]; [5.0 9.0]; [6.0 10.0]; [7.0 11.0];
   [8.0 12.0]; [9.0 12.0]]',
   [[0 0]; [0 0]; [0 0]; [0 0];
   [4 4]; [7 7]; [9 9]; [10 10];
   [10 10]; [9 9]; [7 7]; [4 4];
   [0 0]; [0 0]]'
)

max_stars_repo_name: gabriel-ss/DynamicTimeWarp

max_issues_repo_name: gabriel-ss/DynamicTimeWarp

max_forks_repo_name: gabriel-ss/DynamicTimeWarp


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    EstimationSolution

A solution to a spatial estimation problem.
"""
struct EstimationSolution{D,M,V}
  domain::D
  mean::M
  variance::V
end

# -------------
# VARIABLE API
# -------------

function Base.getindex(solution::EstimationSolution, var::Symbol)
  vars = (var, Symbol(var,"var"))
  vals = (solution.mean[var], solution.variance[var])
  georef((; zip(vars, vals)...), solution.domain)
end

# ------------
# IO methods
# ------------
function Base.show(io::IO, solution::EstimationSolution)
  N = ncoords(solution.domain)
  print(io, "$(N)D EstimationSolution")
end

function Base.show(io::IO, ::MIME"text/plain", solution::EstimationSolution)
  println(io, solution)
  println(io, "  domain: ", solution.domain)
  print(  io, "  variables: ", join(keys(solution.mean), ", ", " and "))
end


text:
#============================================

Benchmarks for insertion into Skiplist

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

using BenchmarkTools
using Random
using Skiplists

function benchmark_insertion(::Type{Skiplist}; N = 10_000)
end

function benchmark_insertion(::Type{SkiplistSet})
end


text:
# using Routing
include("../src/Routing_include.jl")
using Test
using DataStructures

solomon_dataset = OrderedDict(
    # "R101_025" => 617.1, # 1 sec
    # "R102_025" => 547.1, # 2 sec
    # "R103_025" => 454.6, # 2 sec
    # "R104_025" => 416.9, # 3 sec
    # "R105_025" => 530.5, # 1 sec
    # "R106_025" => 465.4, # 8 sec
    # "R107_025" => 424.3, # 2 sec
    # "R108_025" => 397.3, # 5 sec
    # "R109_025" => 441.3, # 1 sec
    # "R110_025" => 444.1, # 30 sec
    # "R111_025" => 428.8, # 13 sec
    # "R112_025" => 393.0, # 15 sec

    # "C101_025" => 191.3, # 1 sec
    # "C102_025" => 190.3, # 29 sec
    # "C103_025" => 190.3, # 51 sec
    # "C104_025" => 186.9, # 150 sec
    # "C105_025" => 191.3, # 1 sec
    # "C106_025" => 191.3, # 1 sec
    # "C107_025" => 191.3, # 3 sec
    # "C108_025" => 191.3, # 5 sec
    # "C109_025" => 191.3, # 13 sec

    # "RC101_025" => 461.1, # 30 sec with limit=10000   --> 4113 sec with limit 400 / 1e-7
    # # "RC102_025" => 351.8,  # oscillation? columns are added forever
    # "RC103_025" => 332.8, # 22 sec
    # "RC104_025" => 306.6, # 88 sec
    # "RC105_025" => 411.3, # 2 sec
    # "RC106_025" => 345.5, # 1 sec
    # "RC107_025" => 298.3, # 11 sec
    # "RC108_025" => 294.5, # 225 sec

    # # "RC101_050" => 944.0, # osciilation?
    # # "RC102_050" => 822.5, # osciilation?
    # # "RC103_050" => 710.9,
    # # "RC104_050" => 545.8,
    # # "RC105_050" => 855.3,
    # # "RC106_050" => 723.2,
    # # "RC107_050" => 642.7,
    # # "RC108_050" => 598.1,

    # "R101_050" => 1044.0, # 12 sec
    # "R102_050" => 909.0, # 10 sec
    # "R103_050" => 772.9, # complete B&B takes very long time
    # "R104_050" => 625.4,
    # "R105_050" => 899.3,
    # "R106_050" => 793.0,
    # "R107_050" => 711.1,
    # "R108_050" => ???,
    # "R109_050" => 786.8,
    # "R110_050" => 697.0,
    # "R111_050" => ???,
    # "R112_050" => ???,

)

comp_time = OrderedDict()

# (name, val) = pop!(solomon_dataset)
for (name, val) in solomon_dataset
    @info("Solving the Solomon Instance $name")
    solomon = load_solomon(name)
    vrptw = generate_solomon_vrptw_instance(solomon)

    @testset "$name" begin
        println("-- Solomon Instance $name --")
        start_time = time()
        @time routes, obj_val = solveVRP(vrptw)
        comp_time[name] = time() - start_time
        @show obj_val
        @test isapprox(obj_val, val, atol=1e-7)
    end
end


text:
# TODO: figure out if there's a way to get around this hacky file-splitting thing
include("./cli.jl")

Base.@ccallable function julia_main(ARGS::Vector{String})::Cint
    main()
    return 0
end


text:
## Show a quick example of how seeding changes the PageRank values

##
using SparseArrays
using DelimitedFiles
using LinearAlgebra
using Printf
##
ei,ej = copy.(eachcol(Int.(readdlm("wiki-simple.edges"))))
xy = readdlm("wiki-simple.xy")'
pages = readlines("wiki-simple.nodes")
n = length(pages)
A = sparse(ei,ej,1,n,n)
## Simple PageRank
function simplepagerank(A,α,v)
  @assert(0 ≤ α < 1, "needs probably α")
  @assert(all(vi -> vi ≥ 0, v), "needs non-negative v")
  v = v ./ sum(v) # we can normalize for them.
  d = vec(sum(A,dims=2)) # compute the degrees
  x = copy(v) # start of with v
  nsteps = 2*ceil(Int,log(eps(1.0))/log(α)) # upper bound on steps
  for i=1:nsteps
    x = α*(A'*(x./d)) .+ (1-α).*v
  end
  return x/sum(x) # renormalize to probability
end
pr = simplepagerank(A,0.85,ones(n)./n)
v = zeros(n); v[1] = 1
spr = simplepagerank(A,0.85,v)

for i=1:10
  @printf("%20s  %.4f  %.4f\n", pages[i], pr[i], spr[i])
end


text:
module independent

using QuantumOpticsBase
using ..interaction, ..system

import ..integrate

# Define Spin 1/2 operators
spinbasis = SpinBasis(1//2)
sigmax_ = sigmax(spinbasis)
sigmay_ = sigmay(spinbasis)
sigmaz_ = sigmaz(spinbasis)
sigmap_ = sigmap(spinbasis)
sigmam_ = sigmam(spinbasis)
I_spin = identityoperator(spinbasis)

"""
    independent.blochstate(phi, theta[, N=1])

Product state of `N` single spin Bloch states.

All spins have the same azimuthal angle `phi` and polar angle `theta`.
"""
function blochstate(phi::Vector{T1}, theta::Vector{T2}) where {T1<:Real, T2<:Real}
    N = length(phi)
    @assert length(theta)==N
    state = zeros(Float64, 3*N)
    state[0*N+1:1*N] = cos(phi).*sin(theta)
    state[1*N+1:2*N] = sin(phi).*sin(theta)
    state[2*N+1:3*N] = cos(theta)
    return state
end

function blochstate(phi::Real, theta::Real, N::Int=1)
    state = zeros(Float64, 3*N)
    state[0*N+1:1*N] = ones(Float64, N)*cos(phi)*sin(theta)
    state[1*N+1:2*N] = ones(Float64, N)*sin(phi)*sin(theta)
    state[2*N+1:3*N] = ones(Float64, N)*cos(theta)
    return state
end

"""
    independent.dim(state)

Number of spins described by this state.
"""
function dim(state::Vector{Float64})
    N, rem = divrem(length(state), 3)
    @assert rem==0
    return N
end

"""
    independent.splitstate(state)

Split state into sx, sy and sz parts.
"""
function splitstate(state::Vector{Float64})
    N = dim(state)
    return view(state, 0*N+1:1*N), view(state, 1*N+1:2*N), view(state, 2*N+1:3*N)
end

"""
    independent.densityoperator(sx, sy, sz)
    independent.densityoperator(state)

Create density operator from independent sigma expectation values.
"""
function densityoperator(sx::Number, sy::Number, sz::Number)
    return 0.5*(identityoperator(spinbasis) + sx*sigmax_ + sy*sigmay_ + sz*sigmaz_)
end
function densityoperator(state::Vector{Float64})
    N = dim(state)
    sx, sy, sz = splitstate(state)
    if N>1
        return DenseOperator(reduce(tensor, [densityoperator(sx[i], sy[i], sz[i]) for i=1:N]))
    else
        return DenseOperator(densityoperator(sx[i], sy[i], sz[i]))
    end
end

"""
    independent.sx(state)

Sigma x expectation values of state.
"""
sx(state::Vector{Float64}) = view(state, 1:dim(state))

"""
    independent.sy(state)

Sigma y expectation values of state.
"""
sy(state::Vector{Float64}) = view(state, dim(state)+1:2*dim(state))

"""
    independent.sz(state)

Sigma z expectation values of state.
"""
sz(state::Vector{Float64}) = view(state, 2*dim(state)+1:3*dim(state))

"""
    independent.timeevolution(T, gamma, state0)

Independent time evolution.

# Arguments
* `T`: Points of time for which output will be generated.
* `gamma`: Single spin decay rate.
* `state0`: Initial state.
"""
function timeevolution(T, gamma::Number, state0::Vector{Float64}; kwargs...)
    N = dim(state0)
    γ = gamma
    function f(ds::Vector{Float64}, s::Vector{Float64}, p, t)
        sx, sy, sz = splitstate(s)
        dsx, dsy, dsz = splitstate(ds)
        @inbounds for k=1:N
            dsx[k] = -0.5*γ*sx[k]
            dsy[k] = -0.5*γ*sy[k]
            dsz[k] = -γ*(1+sz[k])
        end
    end

    fout_(t::Float64, u::Vector{Float64}) = deepcopy(u)

    return integrate(T, f, state0, fout_; kwargs...)
end

"""
    independent.timeevolution(T, S::SpinCollection, state0)

Independent time evolution.

# Arguments
* `T`: Points of time for which output will be generated.
* `S`: SpinCollection describing the system.
* `state0`: Initial state.
"""
timeevolution(T, S::system.SpinCollection, state0::Vector{Float64}) = timeevolution(T, S.gamma, state0)

end # module


text:
## LEGACY METHODS FOR EXPORTING LEARNING NETWORKS BY HAND

# needed?
# function fitresults!(; kwargs...)
#     mach =  machine(; kwargs...) |> fit!
#     return mach.fitresult, mach.cache, mach.report
# end

# legacy method:
function fitresults(yhat::AbstractNode)
    Base.depwarn("`fitresults(::Node)` is deprecated. "*
                 "See \"Composing Models\" section of MLJ manual "*
                 "on preferred way to export learning networks "*
                 "by hand. ", Base.Core.Typeof(fitresults).name.mt.name)
    inputs = sources(yhat, kind=:input)
    targets = sources(yhat, kind=:target)
    weights = sources(yhat, kind=:weights)

    length(inputs) == 1 ||
        error("Improperly exported supervised network does "*
              "not have a unique :input source. ")
    length(targets) < 2 ||
        error("Improperly exported network has multiple :target sources. ")
    length(weights) < 2 ||
        error("Improperly exported network has multiple :weights sources. ")

    cache = anonymize!(vcat(inputs, targets, weights))

    rep = report(yhat)

    return yhat, cache, rep
end

## FALL-BACK METHODS FOR COMPOSITE MODELS (EXPORTED LEARNING NETWORKS)

# *Note.* Be sure to read Note 4 in src/operations.jl to see see how
# fallbacks are provided for operations acting on Composite models.

fitted_params(::Union{Composite,Surrogate},
              fitresult::NamedTuple) =
                  fitted_params(glb(values(fitresult)...))

# legacy code (fitresult a Node):
fitted_params(::Union{Composite,Surrogate},
              fitresult::Node) =
                  fitted_params(fitresult)

function update(model::Union{Composite,Surrogate},
                verbosity::Integer,
                fitresult::NamedTuple,
                cache,
                args...)

    # If any `model` field has been replaced (and not just mutated)
    # then we actually need to fit rather than update (which will
    # force build of a new learning network). If `model` has been
    # created using a learning network export macro, the test used
    # below is perfect. In any other case it is at least conservative,
    # ie it might force a cold restart when a warm restart suffices.

    # greatest lower bound of all nodes delivering predictions:
    glb_node = glb(values(fitresult)...)

    network_model_ids = objectid.(models(glb_node))
    field_values =
        [getproperty(model, name) for name in fieldnames(typeof(model))]
    submodels    = filter(f->f isa Model, field_values)
    submodel_ids = objectid.(submodels)
    if !issubset(submodel_ids, network_model_ids)
        return fit(model, verbosity, args...)
    end

    is_anonymized = cache isa NamedTuple{(:sources, :data)}

    if is_anonymized
        sources, data = cache.sources, cache.data
        for k in eachindex(sources)
            rebind!(sources[k], data[k])
        end
    end

    fit!(glb_node; verbosity=verbosity)
    if is_anonymized
        for s in sources
            rebind!(s, nothing)
        end
    end

    return fitresult, cache, report(glb_node)
end

# legacy version of above (fitresult a Node):
function update(model::Union{Composite,Surrogate},
                verbosity::Integer,
                yhat::Node,
                cache,
                args...)

    # If any `model` field has been replaced (and not just mutated)
    # then we actually need to fit rather than update (which will
    # force build of a new learning network). If `model` has been
    # created using a learning network export macro, the test used
    # below is perfect. In any other case it is at least conservative:
    network_model_ids = objectid.(models(yhat))
    fields = [getproperty(model, name) for name in fieldnames(typeof(model))]
    submodels    = filter(f->f isa Model, fields)
    submodel_ids = objectid.(submodels)
    if !issubset(submodel_ids, network_model_ids)
        return fit(model, verbosity, args...)
    end

    is_anonymized = cache isa NamedTuple{(:sources, :data)}

    if is_anonymized
        sources, data = cache.sources, cache.data
        for k in eachindex(sources)
            rebind!(sources[k], data[k])
        end
    end

    fit!(yhat; verbosity=verbosity)
    if is_anonymized
        for s in sources
            rebind!(s, nothing)
        end
    end

    return yhat, cache, report(yhat)
end

# legacy method (replacements defined in operations.jl):
predict(::SupervisedComposite, fitresult::Node, Xnew)     = fitresult(Xnew)

# legacy method (replacements defined in operations.jl):
transform(::UnsupervisedComposite, fitresult::Node, Xnew) = fitresult(Xnew)


text:
value(x) = x

function __init__()
  @require Distributions="31c24e10-a181-5473-b8eb-7969acd0382f" begin
    handle_distribution_u0(_u0::Distributions.Sampleable) = rand(_u0)
  end

  @require ForwardDiff="f6369f11-7733-5829-9624-2563aa707210" begin

    value(x::ForwardDiff.Dual) = ForwardDiff.value(x)

    # Support adaptive with non-dual time
    @inline function ODE_DEFAULT_NORM(u::AbstractArray{<:ForwardDiff.Dual,N},t) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((ForwardDiff.value(x) for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline function ODE_DEFAULT_NORM(u::Array{<:ForwardDiff.Dual,N},t) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((ForwardDiff.value(x) for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline ODE_DEFAULT_NORM(u::ForwardDiff.Dual,t) = abs(ForwardDiff.value(u))

    # When time is dual, it shouldn't drop the duals for adaptivity
    @inline function ODE_DEFAULT_NORM(u::AbstractArray{<:ForwardDiff.Dual,N},t::ForwardDiff.Dual) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((x for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline function ODE_DEFAULT_NORM(u::Array{<:ForwardDiff.Dual,N},t::ForwardDiff.Dual) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((x for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline ODE_DEFAULT_NORM(u::ForwardDiff.Dual,t::ForwardDiff.Dual) = abs(u)

    # Type piracy. Should upstream
    Base.nextfloat(d::ForwardDiff.Dual{T,V,N}) where {T,V,N} = ForwardDiff.Dual{T}(nextfloat(d.value), d.partials)
    Base.prevfloat(d::ForwardDiff.Dual{T,V,N}) where {T,V,N} = ForwardDiff.Dual{T}(prevfloat(d.value), d.partials)
  end

  @require Measurements="eff96d63-e80a-5855-80a2-b1b0885c5ab7" begin

    value(x::Measurements.Measurement) = Measurements.value(x)

    # Support adaptive steps should be errorless
    @inline function ODE_DEFAULT_NORM(u::AbstractArray{<:Measurements.Measurement,N},t) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((Measurements.value(x) for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline function ODE_DEFAULT_NORM(u::Array{<:Measurements.Measurement,N},t) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((Measurements.value(x) for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline ODE_DEFAULT_NORM(u::Measurements.Measurement,t) = abs(Measurements.value(u))
  end

  @require Unitful="1986cc42-f94f-5a68-af5c-568840ba703d" begin
    # Support adaptive errors should be errorless for exponentiation
    value(x::Unitful.Quantity) = x.val
    @inline function ODE_DEFAULT_NORM(u::AbstractArray{<:Unitful.Quantity,N},t) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline function ODE_DEFAULT_NORM(u::Array{<:Unitful.Quantity,N},t) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline ODE_DEFAULT_NORM(u::Unitful.Quantity,t) = abs(value(u))
  end

  @require Flux="587475ba-b771-5e3f-ad9e-33799f191a9c" begin
    value(x::Flux.Tracker.TrackedReal)  = x.data
    value(x::Flux.Tracker.TrackedArray) = x.data

    # Support adaptive with non-tracked time
    @inline function ODE_DEFAULT_NORM(u::Flux.Tracker.TrackedArray,t) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline function ODE_DEFAULT_NORM(u::AbstractArray{<:Flux.Tracker.TrackedReal,N},t) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline function ODE_DEFAULT_NORM(u::Array{<:Flux.Tracker.TrackedReal,N},t) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline ODE_DEFAULT_NORM(u::Flux.Tracker.TrackedReal,t) = abs(value(u))

    # Support TrackedReal time, don't drop tracking on the adaptivity there
    @inline function ODE_DEFAULT_NORM(u::Flux.Tracker.TrackedArray,t::Flux.Tracker.TrackedReal) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip(u,Iterators.repeated(t))) / length(u))
    end
    @inline function ODE_DEFAULT_NORM(u::AbstractArray{<:Flux.Tracker.TrackedReal,N},t::Flux.Tracker.TrackedReal) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip(u,Iterators.repeated(t))) / length(u))
    end
    @inline function ODE_DEFAULT_NORM(u::Array{<:Flux.Tracker.TrackedReal,N},t::Flux.Tracker.TrackedReal) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip(u,Iterators.repeated(t))) / length(u))
    end
    @inline ODE_DEFAULT_NORM(u::Flux.Tracker.TrackedReal,t::Flux.Tracker.TrackedReal) = abs(u)
  end

end


text:
### A Pluto.jl notebook ###
# v0.16.4

using Markdown
using InteractiveUtils

# ╔═╡ dab179ae-4a5a-11eb-317c-c7fa9d9ccf8f
begin
	import Pkg
	Pkg.activate(Base.current_project())
	using ReinforcementLearning
end

# ╔═╡ 109c4fb2-4a5b-11eb-08d5-bd6b1eb0ebe9
using Plots

# ╔═╡ 1fcd93f0-4a5c-11eb-252d-9da5bc78b08b
using StatsPlots

# ╔═╡ 1fbc2952-4b1b-11eb-3b65-75c1058a9537
using Flux

# ╔═╡ db64341a-4b1b-11eb-3f7b-f11b26f442f4
using Statistics

# ╔═╡ 538080e4-4a5a-11eb-0570-65614c5797f0
md"""
# Ten Armed Bandits Environment

In this chapter, we'll use the `MultiArmBanditsEnv` to study two main concepts in reinforcement learning: **exploration** and **exploitation**.

Let's take a look at the environment first.
"""

# ╔═╡ dfe484d4-4a5a-11eb-0224-573d091b3d08
env = MultiArmBanditsEnv()

# ╔═╡ 1427b132-4a5b-11eb-3506-e744a9a5595c
violin(
	[
		[
			begin
				reset!(env)
				env(a)
				reward(env)
			end
			for _ in 1:100
		]
		for a in action_space(env)
	],
	leg=false
)

# ╔═╡ 69bc9e66-4a5c-11eb-0288-1930cdb31d9d
md"""
The above figure shows the reward distribution of each action. （Figure 2.1）
"""

# ╔═╡ c0ca4172-4aac-11eb-255d-8b0005441fb0
md"""
Now we create a testbed to calulate the average reward and perfect action percentage.
"""

# ╔═╡ 4bf0f782-4aad-11eb-291c-afa853f150a3
"""
A customized hook to record whether the action to take is the best action or not.
"""
begin
	Base.@kwdef struct CollectBestActions <: AbstractHook
		best_action::Int
		isbest::Vector{Bool} = []
	end

	function (h::CollectBestActions)(::PreActStage, agent, env, action)
		push!(h.isbest, h.best_action==action)
	end
end

# ╔═╡ 1ff8d726-4aad-11eb-0d88-c7f6080c4072
md"""
Writing a customized hook is easy.

1. Define your `struct` and make it inherit from `AbstractHook` (optional).
1. Write your customized runtime logic by overwriting some of the following functions. By default, they will do nothing if your hook inherits from `AbstractHook`.
    - `(h::YourHook)(::PreActStage, agent, env, action)`
    - `(h::YourHook)(::PostActStage, agent, env)`
    - `(h::YourHook)(::PreEpisodeStage, agent, env)`
    - `(h::YourHook)(::PostEpisodeStage, agent, env)`
"""

# ╔═╡ 1dfcd040-4b1a-11eb-248e-990c0e029f43
function bandit_testbed(
	;explorer=EpsilonGreedyExplorer(0.1),
	true_reward=0.0,
	init=0.,
	opt=InvDecay(1.0)
)
   env = MultiArmBanditsEnv(;true_reward=true_reward)
   agent = Agent(
	   policy=QBasedPolicy(
		   learner = TDLearner(
			   approximator = TabularQApproximator(
				   n_state=length(state_space(env)),
				   n_action=length(action_space(env)),
				   init=init,
				   opt = opt
			   ),
			   γ = 1.0,
			   method=:SARSA,
			   n = 0,
		   ),
		   explorer = explorer
	   ),
	   trajectory=VectorSARTTrajectory()
	)
	h1 = CollectBestActions(;best_action=findmax(env.true_values)[2])
	h2 = TotalRewardPerEpisode(;is_display_on_exit=false)
	run(agent, env, StopAfterStep(1000), ComposedHook(h1, h2))
    h1.isbest, h2.rewards
end

# ╔═╡ 15e65f7e-4b1b-11eb-26b7-ef85168b5112
begin
	p = plot(layout=(2, 1))
	for ϵ in [0.1, 0.01, 0.0]
    	stats = [
			bandit_testbed(;explorer=EpsilonGreedyExplorer(ϵ))
			for _ in 1:2000
		]
    	plot!(p, mean(x[2] for x in stats);
			subplot=1, legend=:bottomright, label="epsilon=$ϵ")
    	plot!(p, mean(x[1] for x in stats);
			subplot=2, legend=:bottomright, label="epsilon=$ϵ")
	end
	p
end

# ╔═╡ b302e1e6-4b1c-11eb-1cb8-91be2a24d5a7
begin
	# now compare the effect of setting init value
	p_2_3 = plot(legend=:bottomright)

	v1 = mean(
		bandit_testbed(
			;explorer=EpsilonGreedyExplorer(0.),
			init=5.,
			opt=Descent(0.1)
		)[1]
		for _ in 1:2000
	)
	plot!(p_2_3, v1, label="Q_1=5, epsilon=0.")

	v2 = mean(
		bandit_testbed(
			;explorer=EpsilonGreedyExplorer(0.1),
			init=0.,
			opt=Descent(0.1)
		)[1]
		for _ in 1:2000
	)
	plot!(p_2_3, v2, label="Q_1=0, epsilon=0.1")
	p_2_3
end

# ╔═╡ 21c988d8-4b80-11eb-15aa-550e091315f5
begin
	p_2_4 = plot(legend=:bottomright)
	plot!(p_2_4, mean(bandit_testbed(;explorer=UCBExplorer(10), opt=Descent(0.1))[2] for _ in 1:5000), label="UpperConfidenceBound, c=2")
	plot!(p_2_4, mean(bandit_testbed(;explorer=EpsilonGreedyExplorer(0.1), opt=Descent(0.1))[2] for _ in 1:5000), label="epsilon-greedy, epsilon=0.1")

	p_2_4
end

# ╔═╡ 5b3f7b36-4b1e-11eb-0035-2bced9aead26
md"""
Similar to the `bandit_testbed` function, we'll create a new function to test the performance of `GradientBanditLearner`.
"""

# ╔═╡ 6d93c3d0-4b1e-11eb-2b41-af6689ba72f4
function gb_bandit_testbed(
	;baseline=0.,
	explorer=WeightedExplorer(is_normalized=true),
	true_reward=0.0,
	init=0.,
	opt=InvDecay(1.0)
)
    env = MultiArmBanditsEnv(;true_reward=true_reward)
	agent = Agent(
	   policy=QBasedPolicy(
		   learner = GradientBanditLearner(
			   approximator = TabularQApproximator(
				   n_state=length(state_space(env)),
				   n_action=length(action_space(env)),
				   init=init,
				   opt = opt
			   ),
			   baseline=baseline
		   ),
		   explorer = explorer
	   ),
	   trajectory=VectorSARTTrajectory()
	)

	h1 = CollectBestActions(;best_action=findmax(env.true_values)[2])
	h2 = TotalRewardPerEpisode(;is_display_on_exit=false)
	run(agent, env, StopAfterStep(1000), ComposedHook(h1, h2))
    h1.isbest, h2.rewards
end

# ╔═╡ 04e8320c-4b1f-11eb-3340-47f7392a8282
md"""
Note that there's a keyword argument named `baseline` in the `GradientBanditLearner`. It can be either a number or a callable function (`reward -> value`). One of such functions mentioned in the book is to calculate the average of seen rewards.
"""

# ╔═╡ b291cb0c-4b1f-11eb-3ee5-cfcdfdcae00b
begin
	Base.@kwdef mutable struct SampleAvg
		t::Int = 0
		avg::Float64 = 0.0
	end

	function (s::SampleAvg)(x)
		s.t += 1
		s.avg += (x - s.avg) / s.t
		s.avg
	end
end

# ╔═╡ e0e72a60-4b1f-11eb-1001-89777fd3d0f7
begin
	baseline = SampleAvg()
	[baseline(x) for x in 1:10]
end

# ╔═╡ 42525d24-4b20-11eb-099c-b10c90af166e
begin
	true_reward = 4.0

	p_2_5 = plot(legend=:bottomright)

	plot!(p_2_5, mean(gb_bandit_testbed(;opt=Descent(0.1), baseline=SampleAvg(), true_reward=true_reward)[1] for _ in 1:2000), label="alpha = 0.1, with baseline")
	plot!(p_2_5, mean(gb_bandit_testbed(;opt=Descent(0.4), baseline=SampleAvg(), true_reward=true_reward)[1] for _ in 1:2000), label="alpha = 0.4, with baseline")
	plot!(p_2_5, mean(gb_bandit_testbed(;opt=Descent(0.1), true_reward=true_reward)[1] for _ in 1:2000), label="alpha = 0.1, without baseline")
	plot!(p_2_5, mean(gb_bandit_testbed(;opt=Descent(0.4), true_reward=true_reward)[1] for _ in 1:2000), label="alpha = 0.4, without baseline")

	p_2_5
end

# ╔═╡ aad675d2-4b80-11eb-3d78-a1ef731d7d8b
begin

	p_2_6 = plot(legend=:topleft)

	plot!(p_2_6, -7:-2, [mean(mean(bandit_testbed(;explorer=EpsilonGreedyExplorer(2.0^i))[2] for _ in 1:2000)) for i in -7:-2], label="epsilon greedy")
	plot!(p_2_6, -5:1, [mean(mean(gb_bandit_testbed(;explorer=WeightedExplorer(is_normalized=true), opt=Descent(2.0^i))[2] for _ in 1:2000)) for i in -5:1], label="gradient")
	plot!(p_2_6, -4:2, [mean(mean(bandit_testbed(;explorer=UCBExplorer(10; c=2.0^i))[2] for _ in 1:2000)) for i in -4:2], label="UCB")
	plot!(p_2_6, -2:2, [mean(mean(bandit_testbed(;explorer=EpsilonGreedyExplorer(0.), init=(2.0^i))[2] for _ in 1:2000)) for i in -2:2], label="greedy with initialization")

	p_2_6
end

# ╔═╡ Cell order:
# ╟─538080e4-4a5a-11eb-0570-65614c5797f0
# ╠═dab179ae-4a5a-11eb-317c-c7fa9d9ccf8f
# ╠═dfe484d4-4a5a-11eb-0224-573d091b3d08
# ╠═109c4fb2-4a5b-11eb-08d5-bd6b1eb0ebe9
# ╠═1fcd93f0-4a5c-11eb-252d-9da5bc78b08b
# ╠═1427b132-4a5b-11eb-3506-e744a9a5595c
# ╟─69bc9e66-4a5c-11eb-0288-1930cdb31d9d
# ╟─c0ca4172-4aac-11eb-255d-8b0005441fb0
# ╠═4bf0f782-4aad-11eb-291c-afa853f150a3
# ╟─1ff8d726-4aad-11eb-0d88-c7f6080c4072
# ╠═1fbc2952-4b1b-11eb-3b65-75c1058a9537
# ╠═db64341a-4b1b-11eb-3f7b-f11b26f442f4
# ╠═1dfcd040-4b1a-11eb-248e-990c0e029f43
# ╠═15e65f7e-4b1b-11eb-26b7-ef85168b5112
# ╠═b302e1e6-4b1c-11eb-1cb8-91be2a24d5a7
# ╠═21c988d8-4b80-11eb-15aa-550e091315f5
# ╠═5b3f7b36-4b1e-11eb-0035-2bced9aead26
# ╠═6d93c3d0-4b1e-11eb-2b41-af6689ba72f4
# ╟─04e8320c-4b1f-11eb-3340-47f7392a8282
# ╠═b291cb0c-4b1f-11eb-3ee5-cfcdfdcae00b
# ╠═e0e72a60-4b1f-11eb-1001-89777fd3d0f7
# ╠═42525d24-4b20-11eb-099c-b10c90af166e
# ╠═aad675d2-4b80-11eb-3d78-a1ef731d7d8b


text:
# 32-bit fixed point; parameter `f` is the number of fraction bits
struct Fixed{T <: Signed,f} <: FixedPoint{T,  f}
    i::T

    # constructor for manipulating the representation;
    # selected by passing an extra dummy argument
    Fixed{T, f}(i::Integer, _) where {T,f} = new{T, f}(i % T)
    Fixed{T, f}(x) where {T,f} = convert(Fixed{T,f}, x)
    Fixed{T, f}(x::Fixed{T,f}) where {T,f} = x
    Fixed{T, f}(x::Char) where {T,f} = throw(ArgumentError("Fixed cannot be constructed from a Char"))
    Fixed{T, f}(x::Complex) where {T,f} = Fixed{T, f}(convert(real(typeof(x)), x))
    Fixed{T, f}(x::Base.TwicePrecision) where {T,f} = Fixed{T, f}(convert(Float64, x))
end

reinterpret(::Type{Fixed{T,f}}, x::T) where {T <: Signed,f} = Fixed{T,f}(x, 0)

typechar(::Type{X}) where {X <: Fixed} = 'Q'
signbits(::Type{X}) where {X <: Fixed} = 1

for T in (Int8, Int16, Int32, Int64)
    for f in 0:sizeof(T)*8-1
        sym = Symbol(String(take!(showtype(_iotypealias, Fixed{T,f}))))
        @eval begin
            const $sym = Fixed{$T,$f}
            export $sym
        end
    end
end

# basic operators
-(x::Fixed{T,f}) where {T,f} = Fixed{T,f}(-x.i,0)
abs(x::Fixed{T,f}) where {T,f} = Fixed{T,f}(abs(x.i),0)

+(x::Fixed{T,f}, y::Fixed{T,f}) where {T,f} = Fixed{T,f}(x.i+y.i,0)
-(x::Fixed{T,f}, y::Fixed{T,f}) where {T,f} = Fixed{T,f}(x.i-y.i,0)

# with truncation:
#*{f}(x::Fixed32{f}, y::Fixed32{f}) = Fixed32{f}(Base.widemul(x.i,y.i)>>f,0)
# with rounding up:
*(x::Fixed{T,f}, y::Fixed{T,f}) where {T,f} = Fixed{T,f}((Base.widemul(x.i,y.i) + (one(widen(T)) << (f-1)))>>f,0)

/(x::Fixed{T,f}, y::Fixed{T,f}) where {T,f} = Fixed{T,f}(div(convert(widen(T), x.i) << f, y.i), 0)

# # conversions and promotions
convert(::Type{Fixed{T,f}}, x::Integer) where {T,f} = Fixed{T,f}(round(T, convert(widen1(T),x)<<f),0)
convert(::Type{Fixed{T,f}}, x::AbstractFloat) where {T,f} = Fixed{T,f}(round(T, trunc(widen1(T),x)<<f + rem(x,1)*(one(widen1(T))<<f)),0)
convert(::Type{Fixed{T,f}}, x::Rational) where {T,f} = Fixed{T,f}(x.num)/Fixed{T,f}(x.den)

rem(x::Integer, ::Type{Fixed{T,f}}) where {T,f} = Fixed{T,f}(rem(x,T)<<f,0)
rem(x::Real,    ::Type{Fixed{T,f}}) where {T,f} = Fixed{T,f}(rem(Integer(trunc(x)),T)<<f + rem(Integer(round(rem(x,1)*(one(widen1(T))<<f))),T),0)

# convert{T,f}(::Type{AbstractFloat}, x::Fixed{T,f}) = convert(floattype(x), x)
float(x::Fixed) = convert(floattype(x), x)

convert(::Type{BigFloat}, x::Fixed{T,f}) where {T,f} =
    convert(BigFloat,x.i>>f) + convert(BigFloat,x.i&(one(widen1(T))<<f - 1))/convert(BigFloat,one(widen1(T))<<f)
convert(::Type{TF}, x::Fixed{T,f}) where {TF <: AbstractFloat,T,f} =
    convert(TF,x.i>>f) + convert(TF,x.i&(one(widen1(T))<<f - 1))/convert(TF,one(widen1(T))<<f)

convert(::Type{Bool}, x::Fixed{T,f}) where {T,f} = x.i!=0
function convert(::Type{Integer}, x::Fixed{T,f}) where {T,f}
    isinteger(x) || throw(InexactError())
    convert(Integer, x.i>>f)
end
function convert(::Type{TI}, x::Fixed{T,f}) where {TI <: Integer,T,f}
    isinteger(x) || throw(InexactError())
    convert(TI, x.i>>f)
end

convert(::Type{TR}, x::Fixed{T,f}) where {TR <: Rational,T,f} =
    convert(TR, x.i>>f + (x.i&(1<<f-1))//(one(widen1(T))<<f))

promote_rule(ft::Type{Fixed{T,f}}, ::Type{TI}) where {T,f,TI <: Integer} = Fixed{T,f}
promote_rule(::Type{Fixed{T,f}}, ::Type{TF}) where {T,f,TF <: AbstractFloat} = TF
promote_rule(::Type{Fixed{T,f}}, ::Type{Rational{TR}}) where {T,f,TR} = Rational{TR}

# TODO: Document and check that it still does the right thing.
decompose(x::Fixed{T,f}) where {T,f} = x.i, -f, 1


text:
using Metaheuristics, Seaborn
pygui(true)

# Define the foxhole function
function  foxhole(x)
    a = [-32 -16 0 16 32 -32 -16 0 16 32 -32 -16 0 16 32 -32 -16 0 16 32 -32 -16 0 16 32;
         -32 -32 -32 -32 -32 -16 -16 -16 -16 -16 0 0 0 0 0 16 16 16 16 16 32 32 32 32 32]
    return 500. .- 1. ./ (0.002 .+ sum( 1. ./ (j .+ 1. .+ (x[1] .- a[1,j]) .^ 6. .+ (x[2] .- a[2,j]) .^ 6.) for j in 1:size(a,2)))
end

# Bounds
lb, ub = [-50, -50], [50, 50]

# Optimize
@elapsed results = Metaheuristics.optimize(foxhole, lb, ub,
                                           Metaheuristics.Clearing(),
                                           options = Metaheuristics.Options(multithreads=false))

# Plots
X, Y = -50:1:50, permutedims(-50:1:50)
plot_surface(repeat(Y,length(Y),1), repeat(X,1,length(X)), foxhole([X,Y]), rstride=2,cstride=1, color="gray")
scatter3D(results.method.population[:,1], results.method.population[:,2], zs=results.method.population[:,5], c="red")


text:
function make_API_call(url::String)
    try
        response = HTTP.get(url)
        return String(response.body)
    catch e
        #return "Error occurred: $(e)"
        return e
    end
end

function read_CSV_url(url::String)
    try
        response = HTTP.get(url)
        df = CSV.read(response.body, DataFrame)
        return df
    catch e
        return e
    end
end


text:
# prototypes for space transforms and cutoffs
function transform end
function transform_d end
function inv_transform end

# ------------------------------------------------------------
#  Abstract polynomial degree business

"""
`AbstractDegree` : object specifying a degree can be called via
`degree(D, arg)` or via `D(arg)`
"""
abstract type AbstractDegree end

(D::AbstractDegree)(args...) = degree(D, args...)

"""
interface functions for `OneParticleBasis`
"""
function add_into_A! end

"""
interface functions for `OneParticleBasis`
"""
function add_into_A_dA! end

"""
`function degree(D::AbstractDegree, arg)` : compute some notion of degree of
the `arg` argument.
"""
function degree end

"""
`function scaling(b, p)`:

a scaling factor for a basis functions ϕ, which gives a rought estimate on
   the magnitude of ∇ᵖϕ e.g.,
```
ϕ = r^n Ylm
```
has scaling factor `n^p + l^p`, though sharper estimates are also possible.
"""
function scaling end

using LinearAlgebra: Diagonal

diagonal_regulariser(basis; diff = 0) = Diagonal(scaling(basis, diff))

"""
every scalar basis must implement this
"""
function rand_radial end

struct NaiveTotalDegree end

# for a one-particle basis function
degree(b::NamedTuple, Deg::NaiveTotalDegree, basis::OneParticleBasis) =
      degree(b, basis)

# for an ν-correlation basis function
# in this case `bb` should be a Vector of NamedTuples
degree(bb, Deg::NaiveTotalDegree, basis::OneParticleBasis) =
      length(bb) == 0 ? 0 : sum( degree(b, basis) for b in bb )


text:
module GaussianProcessODEs

using DifferentialEquations
using Zygote
using Distances
using LinearAlgebra
using Optim
using QuadGK
using KernelFunctions
using Measurements
using ThreadsX

include("kernelfuns.jl")
include("kernelutils.jl")
include("gp_de.jl")
include("kernel_opt.jl")

include("base.jl")
include("problems.jl")
include("utils.jl")
include("hamiltonian.jl")

end # module


text:
module ResampleAndFit

using DataFrames
using Dates
using Statistics
import FileTools
import LsqFit
import DSP

include("fitdata.jl")
include("resampledata.jl")
include("interpdata.jl")
include("filtdata.jl")
include("fillnans.jl")
include("correctdata.jl")
include("mergedata.jl")

export aggregate2, time2regular, isregular, getresolution, cut2interval!
export interpdf, interp1, meshgrid, mesh2vec, interp2
export fitexp, evalexp, fitpoly, evalpoly
export mmconv, findblocks, filtblocks, demean, findnanblocks, detrend, defirst
export fillnans, replacenans!, na2nan!, missing2nan
export correctinterval, correctinterval!, prepcorrpar
export mergetimeseries

end #module


text:
using Pivot
import Pivot: Static, Filter, Security, Logger, Rendering
import Security: setin_session!, getin_session, setin_cookie!, getin_cookie

# Setting up a relative directory for a home path
appdir = @__FILE__() |> abspath |> dirname

app = Engine()

# middleware -------------------------------------------------------------------
use!(app,
     Logger.simple,
     Rendering.template(appdir),
     Static.public(appdir;
                   ignore=[]),
     Static.favicon,
     Filter.cookie_todict,
     Filter.body_todict,
     Filter.query_todict,
     Security.session(),
     Security.secret("ap!@#\$%^&*:""ofij3209f20930239jafa)()(@#\$#@)"),
     Security.csrf)
# custom middleware
use!(app) do app, ctx
    # lets do something to the context
    ctx.data[:awesome] = "Totally"
    app(ctx)
end
# ------------------------------------------------------------------------------

handle!(app, GET, "/") do ctx
    Static.render(ctx, "index.html",
                  Dict("name"       => "alexander",
                       "csrf_token" => getin_session(ctx, "csrf_token")))
end

# SESSIONS ---------------------------------------------------------------------
handle!(app, GET, "/:name") do ctx
    name = ctx.data[:params]["name"]
    try
        getin_session(ctx, name)
    catch
        "No key in session with id $name"
    end
end

handle!(app, POST, "/:name") do ctx
    name = ctx.data[:params]["name"]
    value = ctx.data[:query]["value"]
    setin_session!(ctx, name, value)
    getin_session(ctx, name)
end
# ------------------------------------------------------------------------------

# COOKIES ----------------------------------------------------------------------
handle!(app, GET, "/cookie/:name") do ctx
    name = ctx.data[:params]["name"]
    try
        getin_cookie(ctx, name)
    catch
        "No key in session with id $name"
    end
end

handle!(app, POST, "/cookie/:name") do ctx
    name = ctx.data[:params]["name"]
    value = ctx.data[:query]["value"]
    setin_cookie!(ctx, name, value)
    getin_cookie(ctx, name)
end
# ------------------------------------------------------------------------------

# FORM -------------------------------------------------------------------------
handle!(app, POST, "/person") do ctx
    Static.render(ctx, "index.html",
                  Dict("name" => "alex"))
end
# ------------------------------------------------------------------------------

handle!(app, GET, "/awesome") do ctx
    ctx.data[:awesome]
end

Pivot.run(app, 8080)


text:
@testset "GapRange" begin
    gr = GapRange(1:4, 5:10)
    @test getindex(gr, 2) == 2
    @test getindex(gr, 6) == 6
    @test getindex(gr, 1:3) == 1:3
    @test getindex(gr, 5:6) == 5:6
    @test getindex(gr, 2:8) == 2:8
    @test getindex(1:20, gr) == 1:10

    @test first(gr) == 1
    @test last(gr) == 10

    @test GapRange(1:5, 6:10) == GapRange(6:10, 1:5)
    @test GapRange(5:-1:1, 10:-1:6) == GapRange(10:-1:6, 5:-1:1)

    for (isf, gr) in (
            (true, GapRange(1, 2:9)),
            (true, GapRange(2:9, 1)),
            (true, GapRange(9, 1:8)),
            (true, GapRange(1:8, 9)),
            (false, GapRange(8:-1:1, 9)),
            (false, GapRange(9:-1:2, 1)),
            (false, GapRange(1, 9:-1:2)),
            (false, GapRange(9, 8:-1:1)))
        @test length(gr) == 9
        if isf
            @test gr == 1:9
        end
    end

    @test length(GapRange(1, 3)) == 2

    @test length(GapRange(3, 1)) == 2

    @test StaticRanges.first_length(GapRange(1, 3)) == 1
    @test StaticRanges.last_length(GapRange(1, 3)) == 1

    @testset "iterate" begin
        for (gr_i, r_i) in zip(GapRange(1:5, 6:10), 1:10)
            @test gr_i == r_i
        end
    end

    @test length(GapRange(1, 2:9)) == 9
    @test length(GapRange(1:8, 9)) == 9
    @testset "errors" begin
        # TODO get gaprange tests working with exceptions
        @test_throws ErrorException GapRange(3, 2:9)
        @test_throws ErrorException GapRange(1:8, 3)
        @test_throws ErrorException GapRange(1:8, 1:8)
        @test_throws ErrorException GapRange(8:-1:1, 3)
        @test_throws ErrorException GapRange(3, 8:-1:1)
        @test_throws ErrorException GapRange(8:-1:1, 8:-1:1)
        @test_throws ErrorException GapRange(8:-1:1, 1:8)
        @test_throws ErrorException GapRange(1:8, 8:-1:1)
        @test_throws ErrorException GapRange(1:5, 10:-1:6)
        @test_throws ErrorException GapRange(5:-1:1, 6:10)
        @test_throws ErrorException GapRange{Int,UnitRange{Int},UnitRange{Float64}}(1:3,UnitRange(4.0, 6.0))
    end

    @test size(GapRange(1:3, 7:10)) == (7,)
    @test size(GapRange(1:5, 6:10)) == (10,)
end


text:
using GraphDatasets: cat_tuple_types

@testset "cat_tuple_types" begin

    @test cat_tuple_types(Tuple{}, Tuple{}) == Tuple{}
    @test cat_tuple_types(Tuple{}, Tuple{Int, String}) == Tuple{Int, String}
    @test cat_tuple_types(Tuple{Int, String}, Tuple{}) == Tuple{Int, String}

    @test cat_tuple_types(Tuple{}, @NamedTuple{a::String, b::Int}) == @NamedTuple{a::String, b::Int}
    @test cat_tuple_types(@NamedTuple{a::String, b::Int}, Tuple{}) == @NamedTuple{a::String, b::Int}

    @test cat_tuple_types(Tuple{Int, String}, Tuple{Char, Vector, Tuple}) == Tuple{Int, String, Char, Vector, Tuple}

    @test cat_tuple_types(@NamedTuple{a::Int, b::String}, @NamedTuple{c::Char, d::Vector, e::Tuple}) ==
        @NamedTuple{a::Int, b::String, c::Char, d::Vector, e::Tuple}

    @test cat_tuple_types(Tuple{Int}, Tuple{}, Tuple{Char, Int}, Tuple{String}) ==
            Tuple{Int, Char, Int, String}

    @test cat_tuple_types(@NamedTuple{a::String}, Tuple{}, @NamedTuple{b::Int}, @NamedTuple{c::String, e::Char}) ==
            @NamedTuple{a::String, b::Int, c::String, e::Char}
end


text:
@testset "Shortest paths" begin
  @testset "Constructor with $i nodes" for i in [2, 5, 10]
    # Edge weights: high gain likelihood for the path 1->2->3...n, low for the other edges.
    n = i
    ε = 1 / (n + 1)
    graph = complete_digraph(n)
    reward = Dict{Tuple{Int64, Int64}, Distribution}(
      (i, j) => Bernoulli((i + 1 == j) ? (1 - ε) : (ε))
      for i in 1:n, j in 1:n if i != j)
    instance = ElementaryPath(graph, reward, 1, n, ElementaryPathNoSolver())

    @test instance.n_arms == ne(graph)
    @test instance.reward == reward

    # Error: source makes no sense.
    @test_throws ErrorException ElementaryPath(graph, reward, - 1, n, ElementaryPathNoSolver())
    @test_throws ErrorException ElementaryPath(graph, reward, n + 1, n, ElementaryPathNoSolver())

    # Error: destination makes no sense.
    @test_throws ErrorException ElementaryPath(graph, reward, 1, - 1, ElementaryPathNoSolver())
    @test_throws ErrorException ElementaryPath(graph, reward, 1, n + 1, ElementaryPathNoSolver())
  end

  @testset "State with $i nodes" for i in [2, 5, 10]
    n = i
    ε = 1 / (n + 1)
    graph = complete_digraph(n)
    reward = Dict{Tuple{Int64, Int64}, Distribution}(
      (i, j) => Bernoulli((i + 1 == j) ? (1 - ε) : (ε))
      for i in 1:n, j in 1:n if i != j)
    instance = ElementaryPath(graph, reward, 1, n, ElementaryPathNoSolver())

    state = initial_state(instance)

    @test state.round == 0
    @test state.regret == 0.0
    @test state.reward == 0.0
    @test length(state.arm_counts) == ne(graph)
    @test length(state.arm_reward) == ne(graph)
    @test length(state.arm_average_reward) == ne(graph)

    for i in 1:n
      for j in 1:n
        if i == j
          continue
        end

        @test state.arm_counts[(i, j)] == 0
        @test state.arm_reward[(i, j)] == 0.0
        @test state.arm_average_reward[(i, j)] == 0.0
      end
    end
  end

  @testset "Trace with $i nodes" for i in [2, 5, 10]
    n = i
    ε = 1 / (n + 1)
    reward = Dict{Tuple{Int64, Int64}, Distribution}(
      (i, j) => Bernoulli((i + 1 == j) ? (1 - ε) : (ε))
      for i in 1:n, j in 1:n if i != j)
    instance = ElementaryPath(complete_digraph(n), reward, 1, n, ElementaryPathNoSolver())

    trace = initial_trace(instance)

    @test length(trace.states) == 0
    @test length(trace.arms) == 0
    @test length(trace.reward) == 0
    @test length(trace.policy_details) == 0
    @test length(trace.time_choose_action) == 0

    @test eltype(trace.states) == State{Tuple{Int, Int}}
    @test eltype(trace.arms) == Vector{Tuple{Int, Int}}
    @test eltype(trace.reward) == Vector{Float64}
    @test eltype(trace.time_choose_action) == Int
  end

  @testset "Pull with $i nodes" for i in [2, 5, 10]
    reward = Dict{Tuple{Int64, Int64}, Distribution}(
      (i, j) => Bernoulli((i == 1 && j == i) ? 1 : 1 / i)
      for i in 1:i, j in 1:i if i != j)
    instance = ElementaryPath(complete_digraph(i), reward, 1, i, ElementaryPathLightGraphsDijkstraSolver())

    Random.seed!(2)
    reward, regret = pull(instance, [(1, i)])
    @test reward ≈ [1.0]
    @test regret ≈ 0.0
  end

  @testset "Check feasibility with 3 nodes" begin
    n = 3
    ε = 1 / (n + 1)
    reward = Dict{Tuple{Int64, Int64}, Distribution}(
      (i, j) => Bernoulli((i + 1 == j) ? (1 - ε) : (ε))
      for i in 1:n, j in 1:n if i != j)
    instance = ElementaryPath(complete_digraph(n), reward, 1, n, ElementaryPathNoSolver())

    @test is_feasible(instance, [(1, 3)])
    @test is_feasible(instance, [(1, 2), (2, 3)])
    @test ! is_feasible(instance, Tuple{Int, Int}[])
    @test ! is_feasible(instance, [(1, 2), (2, 3), (3, 1)])
    @test ! is_feasible(instance, [(2, 3), (3, 2), (1, 2)])
    @test ! is_feasible(instance, [(2, 3), (3, 2), (2, 1), (1, 2)])
  end

  @testset "Check partial acceptability with 3 nodes" begin
    n = 3
    ε = 1 / (n + 1)
    reward = Dict{Tuple{Int64, Int64}, Distribution}(
      (i, j) => Bernoulli((i + 1 == j) ? (1 - ε) : (ε))
      for i in 1:n, j in 1:n if i != j)
    instance = ElementaryPath(complete_digraph(n), reward, 1, n, ElementaryPathNoSolver())

    @test is_partially_acceptable(instance, [(1, 3)])
    @test is_partially_acceptable(instance, [(1, 2)])
    @test is_partially_acceptable(instance, [(1, 2), (2, 3)])
    @test is_partially_acceptable(instance, Tuple{Int, Int}[])
    @test ! is_partially_acceptable(instance, [(1, 2), (2, 3), (3, 1)])
    @test ! is_partially_acceptable(instance, [(2, 3), (3, 2), (1, 2)])
    @test ! is_partially_acceptable(instance, [(2, 3), (3, 2), (2, 1), (1, 2)])
  end

  @testset "LightGraphs.jl Dijkstra solver" begin
    @testset "Constructor with $i nodes" for i in [2, 5, 10]
      n = i
      ε = 1 / (n + 1)
      reward = Dict{Tuple{Int64, Int64}, Distribution}(
        (i, j) => Bernoulli((i + 1 == j) ? (1 - ε) : (ε))
        for i in 1:n, j in 1:n if i != j)
      instance = ElementaryPath(complete_digraph(n), reward, 1, n, ElementaryPathLightGraphsDijkstraSolver())

      @test instance.solver != nothing
      @test instance.solver.graph != nothing
      @test size(instance.solver.weights_matrix, 1) == n
      @test size(instance.solver.weights_matrix, 2) == n
      @test instance.solver.source == 1
      @test instance.solver.destination == n
    end

    @testset "Solve with $i nodes" for i in [2, 5, 10]
      n = i
      ε = 1 / (n + 1)
      reward = Dict{Tuple{Int64, Int64}, Distribution}(
        (i, j) => Bernoulli((i + 1 == j) ? (1 - ε) : (ε))
        for i in 1:n, j in 1:n if i != j)
      instance = ElementaryPath(complete_digraph(n), reward, 1, n, ElementaryPathLightGraphsDijkstraSolver())

      Random.seed!(i)
      drawn = Dict((i, j) => rand() for i in 1:n, j in 1:n)
      solution = solve_linear(instance, drawn)
      @test is_feasible(instance, solution)
    end
  end

  if ! is_travis
    @testset "LP solver" begin
      @testset "Constructor with $i nodes" for i in [5]#[2, 5, 10]
        n = i
        ε = 1 / (n + 1)
        reward = Dict{Tuple{Int64, Int64}, Distribution}(
          (i, j) => Bernoulli((i + 1 == j) ? (1 - ε) : (ε))
          for i in 1:n, j in 1:n if i != j)
        instance = ElementaryPath(complete_digraph(n), reward, 1, n, ElementaryPathLPSolver(Gurobi.Optimizer))

        @test instance.solver != nothing
        @test instance.solver.graph != nothing
        @test instance.solver.source == 1
        @test instance.solver.destination == n
      end

      @testset "Solve with $i nodes" for i in [2, 5, 10]
        n = i
        ε = 1 / (n + 1)
        reward = Dict{Tuple{Int64, Int64}, Distribution}(
          (i, j) => Bernoulli((i + 1 == j) ? (1 - ε) : (ε))
          for i in 1:n, j in 1:n if i != j)
        instance = ElementaryPath(complete_digraph(n), reward, 1, n, ElementaryPathLPSolver(Gurobi.Optimizer))

        Random.seed!(i)
        drawn = Dict((i, j) => rand() for i in 1:n, j in 1:n if i != j)
        solution = solve_linear(instance, drawn)
        @test is_feasible(instance, solution)
      end
    end
  end

  @testset "Algos solver" begin
    @testset "Constructor with $i nodes" for i in [2, 5, 10]
      n = i
      ε = 1 / (n + 1)
      reward = Dict{Tuple{Int64, Int64}, Distribution}(
        (i, j) => Bernoulli((i + 1 == j) ? (1 - ε) : (ε))
        for i in 1:n, j in 1:n if i != j)
      instance = ElementaryPath(complete_digraph(n), reward, 1, n, ElementaryPathAlgosSolver())

      @test instance.solver != nothing
      @test instance.solver.graph != nothing
      @test instance.solver.source == 1
      @test instance.solver.destination == n
    end

    @testset "Solve with $i nodes" for i in [2, 5, 10]
      n = i
      ε = 1 / (n + 1)
      reward = Dict{Tuple{Int64, Int64}, Distribution}(
        (i, j) => Bernoulli((i + 1 == j) ? (1 - ε) : (ε))
        for i in 1:n, j in 1:n if i != j)
      instance = ElementaryPath(complete_digraph(n), reward, 1, n, ElementaryPathAlgosSolver())

      Random.seed!(i)
      drawn = Dict((i, j) => rand() for i in 1:n, j in 1:n if i != j)
      solution = solve_linear(instance, drawn)
      @test is_feasible(instance, solution)
    end

    @testset "Solve with $i nodes and a budget" for i in [3]# [2, 5, 10]
      # Test case: all edges have a 100% probability of reward and a zero weight,
      # except the direct edge from the source to the destination (no reward,
      # some weight). The algorithm is looking for a path with a nonzero minimum
      # weight, i.e. there is only one possible solution (1 -> n) with zero reward.
      n = i
      ε = 1 / (n + 1)
      reward = Dict{Tuple{Int64, Int64}, Distribution}(
        (i, j) => Bernoulli((i == 1 && j == n) ? 0 : 1)
        for i in 1:n, j in 1:n if i != j)
      instance = ElementaryPath(complete_digraph(n), reward, 1, n, ElementaryPathAlgosSolver())

      # There is no real randomness in this test case, hence no need to set the seed.
      rewards = Dict((i, j) => Float64(rand(reward[i, j])) for i in 1:n, j in 1:n if i != j)
      weights = Dict((i, j) => (i == 1 && j == n) ? 2 : 0 for i in 1:n, j in 1:n if i != j)
      solution = solve_budgeted_linear(instance, rewards, weights, 2)
      @test solution == [(1, n)]
    end
  end

  if ! is_travis
    @testset "Solver equivalence" begin
      function compare_lp_algos(graph, s, d, rewards, weights, ε, budgets)
        ε = 0.003
        rewards_discrete = Dict(k => round(Int, v / ε, RoundUp) for (k, v) in rewards)
        distr = Dict{Tuple{Int, Int}, Distribution}((src(e), dst(e)) => Bernoulli(0.0) for e in edges(graph))

        @testset "Budget $budget" for budget in budgets
          i_lp = ElementaryPath(graph, distr, s, d, ElementaryPathLPSolver(Gurobi.Optimizer))
          i_al = ElementaryPath(graph, distr, s, d, ElementaryPathAlgosSolver())
          s_lp = ! is_travis && solve_budgeted_linear(i_lp.solver, weights, rewards_discrete, budget)
          s_al = solve_budgeted_linear(i_al.solver, weights, rewards_discrete, budget)

          if length(s_lp) == 0
            @test length(s_lp) == length(s_al)
            continue
          end
          @test length(s_lp) >= 1
          @test length(s_al) >= 1

          @test sum(rewards_discrete[arm] for arm in s_lp) >= budget
          @test sum(rewards_discrete[arm] for arm in s_al) >= budget

          o_lp = sum(rewards[arm] for arm in s_lp) + sqrt(sum(weights[arm] for arm in s_lp))
          o_al = sum(rewards[arm] for arm in s_al) + sqrt(sum(weights[arm] for arm in s_al))
          @test o_al ≈ o_lp atol=ε

          # Compare directly against Bellman-Ford.
          if budget == 0
            weights_edge = Dict(Edge(k...) => v for (k, v) in weights)
            isp = CombinatorialBandits.ElementaryPathInstance(graph, weights_edge, s, d)
            s_bf = CombinatorialBandits.lp_dp(isp).path
            o_bf = sum(rewards[src(arm), dst(arm)] for arm in s_bf) + sqrt(sum(weights[src(arm), dst(arm)] for arm in s_bf))

            @test o_al ≈ o_bf atol=ε
          end
        end

        @testset "Overall solution" begin
          i_lp = ElementaryPath(graph, distr, s, d, ElementaryPathLPSolver(Gurobi.Optimizer))
          i_al = ElementaryPath(graph, distr, s, d, ElementaryPathAlgosSolver())
          s_lp = CombinatorialBandits.optimise_linear_sqrtlinear(i_lp, ESCB2Exact(), rewards, weights)
          s_al = CombinatorialBandits.optimise_linear_sqrtlinear(i_al, ESCB2Budgeted(ε, true), rewards, weights)
          o_lp = sum(rewards[arm] for arm in s_lp) + sqrt(sum(weights[arm] for arm in s_lp))
          o_al = sum(rewards[arm] for arm in s_al) + sqrt(sum(weights[arm] for arm in s_al))

          @test o_lp ≈ o_al atol=ε
        end
      end

      function get_graph(n)
        graph = DiGraph(n)
        for i in 1:(n - 1)
          for j in (i + 1):n
            add_edge!(graph, i, j)
          end
        end
        return graph, 1, n
      end

      @testset "First case" begin
        graph, s, d = get_graph(5)
        ε = 0.003
        rewards = Dict((1, 2) => 0.0,(1, 3) => 0.0,(1, 4) => 1.0,(1, 5) => 1.0,(2, 3) => 0.3333333333333333,(2, 4) => 0.0,(2, 5) => 0.0,(3, 4) => 0.3333333333333333,(3, 5) => 0.0,(4, 5) => 0.0)
        weights = Dict((1, 3) => 15.681708132568822,(1, 2) => 3.1363416265137642,(1, 4) => 15.681708132568822,(1, 5) => 15.681708132568822,(2, 3) => 5.227236044189607,(2, 4) => 15.681708132568822,(2, 5) => 15.681708132568822,(3, 4) => 5.227236044189607,(3, 5) => 15.681708132568822,(4, 5) => 3.1363416265137642)

        compare_lp_algos(graph, s, d, rewards, weights, ε, [0, 1, 112, 113, 224, 225, 334, 335])
         # No more solution after 334
      end

      @testset "Second case" begin
        graph, s, d = get_graph(5)
        ε = 0.003
        rewards = Dict((1, 2) => 0.18126888217522658,(1, 3) => 0.1527777777777778,(1, 4) => 0.23404255319148937,(1, 5) => 0.8253968253968254,(2, 3) => 0.13100436681222707,(2, 4) => 0.12307692307692308,(2, 5) => 0.10810810810810811,(3, 4) => 0.14847161572052403,(3, 5) => 0.16666666666666666,(4, 5) => 0.16129032258064516)
        weights = Dict((1, 2) => 0.12134857271805984,(1, 3) => 0.557866355134414,(1, 4) => 0.8546037780782512,(1, 5) => 0.3187807743625223,(2, 3) => 0.17539902868854937,(2, 4) => 0.6179442703027355,(2, 5) => 1.0855777721534543,(3, 4) => 0.17539902868854937,(3, 5) => 0.557866355134414,(4, 5) => 0.11778996354744224)

        compare_lp_algos(graph, s, d, rewards, weights, ε, [0, 1, 37, 38, 42, 43, 44, 45, 50, 51, 52, 54, 55, 56, 61, 62, 79, 80, 82, 83, 87, 88, 89, 91, 92, 93, 94, 95, 98, 99, 276, 277, 552, 553, 828, 829])
      end

      @testset "Third case" begin
        graph, s, d = get_graph(10)
        ε = 1/63
        rewards = Dict((1, 2) => 0.08148148148148149,(1, 3) => 0.0958904109589041,(1, 4) => 0.09090909090909091,(1, 5) => 0.05,(1, 6) => 0.125,(1, 7) => 0.0,(1, 8) => 0.09090909090909091,(1, 9) => 0.2,(1, 10) => 0.0, (2, 3) => 0.07746478873239436,(2, 4) => 0.02631578947368421,(2, 5) => 0.11538461538461539,(2, 6) => 0.1111111111111111,(2, 7) => 0.2,(2, 8) => 0.0,(2, 9) => 0.2727272727272727,(2, 10) => 0.0,(3, 4) => 0.14953271028037382,(3, 5) => 0.0625,(3, 6) => 0.2,(3, 7) => 0.0,(3, 8) => 0.07692307692307693,(3, 9) => 0.25,(3, 10) => 0.09090909090909091,(4, 5) => 0.11764705882352941,(4, 6) => 0.0625,(4, 7) => 0.0,(4, 8) => 0.0625,(4, 9) => 0.07142857142857142,(4, 10) => 0.0,(5, 6) => 0.06666666666666667,(5, 7) => 0.09090909090909091,(5, 8) => 0.0, (5, 9) => 0.05555555555555555,(5, 10) => 0.125,(6, 7) => 0.14772727272727273,(6, 8) => 0.06060606060606061,(6, 9) => 0.041666666666666664,(6, 10) => 0.047619047619047616,(7, 8) => 0.12149532710280374,(7, 9) => 0.047619047619047616,(7, 10) => 0.08823529411764706,(8, 9) => 0.06716417910447761,(8, 10) => 0.11538461538461539,(9, 10) => 0.0830188679245283)
        weights = Dict((1, 2) => 0.6156966918752366,(1, 3) => 2.277234339812519,(1, 4) => 5.037518388070118,(1, 5) => 8.311905340315693,(1, 6) => 10.389881675394617,(1, 7) => 12.78754667740876,(1, 8) => 15.112555164210352,(1, 9) => 16.623810680631387,(1, 10) => 12.78754667740876,(2, 3) => 1.1706908930022104,(2, 4) => 4.3746870212187865,(2, 5) => 6.39377333870438,(2, 6) => 9.235450378128549,(2, 7) => 11.082540453754259,(2, 8) => 15.112555164210352,(2, 9) => 15.112555164210352,(2, 10) => 18.470900756257098,(3, 4) => 1.5536271664141483,(3, 5) => 5.194940837697309,(3, 6) => 6.649524272252555,(3, 7) => 11.082540453754259,(3, 8) => 12.78754667740876,(3, 9) => 13.853175567192823,(3, 10) => 15.112555164210352,(4, 5) => 1.9557424330154574,(4, 6) => 5.194940837697309,(4, 7) => 8.749374042437573,(4, 8) => 10.389881675394617,(4, 9) => 11.874150486165277,(4, 10) => 13.853175567192823,(5, 6) => 2.216508090750852,(5, 7) => 5.037518388070118,(5, 8) => 7.916100324110185,(5, 9) => 9.235450378128549,(5, 10) => 10.389881675394617,(6, 7) => 1.889069395526294,(6, 8) => 5.037518388070118,(6, 9) => 6.926587783596411,(6, 10) => 7.916100324110185,(7, 8) => 1.5536271664141483,(7, 9) => 3.9580501620550925,(7, 10) => 4.889356082538644,(8, 9) => 1.2405828866142827,(8, 10) => 2.1312577795681267,(9, 10) => 0.6273136105898637)

        compare_lp_algos(graph, s, d, rewards, weights, ε, collect(0:20))
      end
    end
  end
end


text:
using ProtoBuf
using ProtoBuf.Gen
gen()


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 5fe9a49a-fdc2-11ea-3ff9-01e000224b3a
using Pkg, DrWatson

# ╔═╡ 5fe9db1c-fdc2-11ea-0c8c-7db1e606cbcb
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanQuap
	using StatisticalRethinking
end

# ╔═╡ 5b3cfe6e-fdc2-11ea-2462-cb4109cb77ee
md"## Clip-05-49.1s.jl"

# ╔═╡ 5fea720a-fdc2-11ea-2dde-493a3df79ff2
begin
	df = CSV.read(sr_datadir("milk.csv"), DataFrame; delim=';')
	df[!, :clade_id] = Int.(indexin(df[:, :clade], unique(df[:, :clade])))
	scale!(df, [:kcal_per_g])
	PRECIS(df[:, [:clade_id, :kcal_per_g]])
end

# ╔═╡ 5ffa4cfc-fdc2-11ea-0dd3-659a65ad4949
stan5_9 = "
data{
  int <lower=1> N;              // Sample size
  int <lower=1> k;
  vector[N] K;
  int clade_id[N];
}
parameters{
  vector[k] a;
  real<lower=0> sigma;
}
model{
  vector[N] mu;
  sigma ~ exponential( 1 );
  a ~ normal( 0 , 0.5 );
  for ( i in 1:N ) {
      mu[i] = a[clade_id[i]];
  }
  K ~ normal( mu , sigma );
}
";

# ╔═╡ 0233ffa6-fdc1-11ea-202f-d18f84a47bb1
md"##### Define the SampleModel."

# ╔═╡ 90cf5430-fdc2-11ea-279a-19fbadfdf16b
begin
	data = (N = size(df, 1), clade_id = df.clade_id,
    	K = df.kcal_per_g_s, k = length(unique(df[:, :clade_id])))
	init = (a = [-1.0, 0,0, 1.0, 1.0], sigma = 1.0,)
	q5_9s, m5_9s, _ = stan_quap("m5.9", stan5_9; data, init)

	if !isnothing(m5_9s)
		post5_9s_df = read_samples(m5_9s, :dataframe)
		PRECIS(post5_9s_df)
	end
end

# ╔═╡ 7e7479dc-fdc2-11ea-2644-a75204ec3851
if !isnothing(q5_9s)
	quap5_9s_df = sample(q5_9s)
	PRECIS(quap5_9s_df)
end

# ╔═╡ 7e7528a0-fdc2-11ea-1031-973876c110eb
rethinking = "
       mean   sd  5.5% 94.5%
a[1]  -0.48 0.22 -0.83 -0.14
a[2]   0.37 0.22  0.02  0.71
a[3]   0.68 0.26  0.26  1.09
a[4]  -0.59 0.27 -1.02 -0.15
sigma  0.72 0.10  0.57  0.87
";

# ╔═╡ 7e89657a-fdc2-11ea-33a9-779fdfd53519
md"## End of clip-05-49.1s.jl"

# ╔═╡ Cell order:
# ╟─5b3cfe6e-fdc2-11ea-2462-cb4109cb77ee
# ╠═5fe9a49a-fdc2-11ea-3ff9-01e000224b3a
# ╠═5fe9db1c-fdc2-11ea-0c8c-7db1e606cbcb
# ╠═5fea720a-fdc2-11ea-2dde-493a3df79ff2
# ╠═5ffa4cfc-fdc2-11ea-0dd3-659a65ad4949
# ╟─0233ffa6-fdc1-11ea-202f-d18f84a47bb1
# ╠═90cf5430-fdc2-11ea-279a-19fbadfdf16b
# ╠═7e7479dc-fdc2-11ea-2644-a75204ec3851
# ╠═7e7528a0-fdc2-11ea-1031-973876c110eb
# ╟─7e89657a-fdc2-11ea-33a9-779fdfd53519


text:
using BenchmarkTools
using Plots
using ProgressBars
using OpenQuantumSystems
# using OpenQuantumSystemsPrivate
using LinearAlgebra
using Random
using QuadGK
using ProgressMeter
using Expokit
# using DifferentialEquations
using DataFrames
using JSON
using Git

import OrdinaryDiffEq, DiffEqCallbacks, DelayDiffEq
import SparseArrays: sparse
import QuantumOpticsBase

BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

# dic_commits = JSON.parsefile("test/benchmark_commits.json")
open("benchmark/benchmark_commits.json", "r") do f
    global dic_commits = JSON.parse(f)
end

commit_hash = readlines(`git log -n 1 --no-merges --format=%h`)[1]
commit_date = readlines(`git log -n 1 --no-merges --format=%cs`)[1]
dic_commit = Dict()
dic_commit["commit_date"] = commit_date

println

## Aggregate building
### agg_dimer_small
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 5.0
t1_1 = @benchmark begin
    HR = 0.01
    shift = (2.0 * HR)
    modes = [Mode(180.0, shift)]
    mols = [
        Molecule([Mode(300.0, shift)], 3, [12500.0, 12710.0]),
        Molecule([Mode(300.0, shift)], 3, [12500.0, 12690.0]),
    ]

    agg = Aggregate(mols)
    for mol_i = 2:length(agg.molecules)
        agg.coupling[mol_i, mol_i+1] = 200
        agg.coupling[mol_i+1, mol_i] = 200
    end

    aggIndices = getIndices(agg)
    vibindices = getVibIndices(agg, aggIndices)
    aggIndLen = length(aggIndices)
    base = GenericBasis([aggIndLen])
    FCFact = getFranckCondonFactors(agg, aggIndices)
    FCProd = getFCProd(agg, FCFact, aggIndices, vibindices)
    Ham = getAggHamiltonian(agg, aggIndices, FCFact)

    Ham_bath = getAggHamiltonianBath(agg)
    Ham_sys = getAggHamiltonianSystem(agg; groundEnergy = false)
    b = GenericBasis([aggIndLen])
    b_sys = GenericBasis([size(Ham_sys, 1)])
    b_bath = GenericBasis([size(Ham_bath, 1)])

    Ham_int = getAggHamiltonianInteraction(agg, aggIndices, FCFact)
    Ham_S = Ham - Ham_int
    b_sys = GenericBasis([size(Ham_sys, 1)])
    b_bath = GenericBasis([size(Ham_bath, 1)])
    Ham_B = tensor(Ham_bath, OneDenseOperator(b_sys))
    E0 = Ham_B.data[1, 1]
    for I = 1:aggIndLen
        Ham_B.data[I, I] -= E0
    end
    Ham_B = DenseOperator(b, b, Ham_B.data)
end

### agg_dimer_medium
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 20.0
t1_2 = @benchmark begin
    HR = 0.01
    shift = (2.0 * HR)
    modes = [Mode(180.0, shift)]
    mols = [
        Molecule([Mode(300.0, shift)], 7, [12500.0, 12710.0]),
        Molecule([Mode(300.0, shift)], 7, [12500.0, 12690.0]),
    ]

    agg = Aggregate(mols)
    for mol_i = 2:length(agg.molecules)
        agg.coupling[mol_i, mol_i+1] = 200
        agg.coupling[mol_i+1, mol_i] = 200
    end

    aggIndices = getIndices(agg)
    vibindices = getVibIndices(agg, aggIndices)
    aggIndLen = length(aggIndices)
    base = GenericBasis([aggIndLen])
    FCFact = getFranckCondonFactors(agg, aggIndices)
    FCProd = getFCProd(agg, FCFact, aggIndices, vibindices)
    Ham = getAggHamiltonian(agg, aggIndices, FCFact)

    Ham_bath = getAggHamiltonianBath(agg)
    Ham_sys = getAggHamiltonianSystem(agg; groundEnergy = false)
    b = GenericBasis([aggIndLen])
    b_sys = GenericBasis([size(Ham_sys, 1)])
    b_bath = GenericBasis([size(Ham_bath, 1)])

    Ham_int = getAggHamiltonianInteraction(agg, aggIndices, FCFact)
    Ham_S = Ham - Ham_int
    b_sys = GenericBasis([size(Ham_sys, 1)])
    b_bath = GenericBasis([size(Ham_bath, 1)])
    Ham_B = tensor(Ham_bath, OneDenseOperator(b_sys))
    E0 = Ham_B.data[1, 1]
    for I = 1:aggIndLen
        Ham_B.data[I, I] -= E0
    end
    Ham_B = DenseOperator(b, b, Ham_B.data)
end

### agg_building
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 1.0
t1_4 = @benchmark begin
    HR = 0.01
    shift = (2.0 * HR)
    modes = [Mode(180.0, shift)]
    mols = [
        Molecule([Mode(300.0, shift)], 7, [12500.0, 12710.0]),
        Molecule([Mode(300.0, shift)], 7, [12500.0, 12690.0]),
    ]

    agg = Aggregate(mols)
end
HR = 0.01
shift = (2.0 * HR)
modes = [Mode(180.0, shift)]
mols = [
    Molecule([Mode(300.0, shift)], 7, [12500.0, 12710.0]),
    Molecule([Mode(300.0, shift)], 7, [12500.0, 12690.0]),
]

agg = Aggregate(mols)
for mol_i = 2:length(agg.molecules)
    agg.coupling[mol_i, mol_i+1] = 200
    agg.coupling[mol_i+1, mol_i] = 200
end

aggIndices = getIndices(agg)
vibindices = getVibIndices(agg, aggIndices)
aggIndLen = length(aggIndices)
base = GenericBasis([aggIndLen])

### getFranckCondonFactors
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 10.0
t1_5 = @benchmark begin
    FCFact = getFranckCondonFactors(agg, aggIndices)
end
FCFact = getFranckCondonFactors(agg, aggIndices)

### getFCProd
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 100.0
t1_6 = @benchmark begin
    FCProd = getFCProd(agg, FCFact, aggIndices, vibindices)
end
FCProd = getFCProd(agg, FCFact, aggIndices, vibindices)

### getFCProd
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 10.0
t1_7 = @benchmark begin
    Ham = getAggHamiltonian(agg, aggIndices, FCFact)
end
Ham = getAggHamiltonian(agg, aggIndices, FCFact)
Ham_bath = getAggHamiltonianBath(agg)
Ham_sys = getAggHamiltonianSystem(agg; groundEnergy = false)
b = GenericBasis([aggIndLen])
b_sys = GenericBasis([size(Ham_sys, 1)])
b_bath = GenericBasis([size(Ham_bath, 1)])

### getAggHamiltonianInteraction
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 10.0
t1_8 = @benchmark begin
    Ham_int = getAggHamiltonianInteraction(agg, aggIndices, FCFact)
end

dic_commit["agg_dimer_small"] = t1_1.times # 27
dic_commit["agg_dimer_medium"] = t1_2.times # 147
dic_commit["agg_building"] = t1_4.times
dic_commit["getFranckCondonFactors"] = t1_5.times
dic_commit["getFCProd"] = t1_6.times
dic_commit["getAggHamiltonian"] = t1_7.times
dic_commit["getAggHamiltonianInteraction"] = t1_8.times

## Simulations
HR = 0.01
shift = (2.0 * HR)
modes = [Mode(180.0, shift)]
mols = [
    Molecule([Mode(300.0, shift)], 3, [12500.0, 12710.0]),
    Molecule([Mode(300.0, shift)], 3, [12500.0, 12690.0]),
]

agg = Aggregate(mols)
for mol_i = 2:length(agg.molecules)
    agg.coupling[mol_i, mol_i+1] = 200
    agg.coupling[mol_i+1, mol_i] = 200
end

aggIndices = getIndices(agg)
vibindices = getVibIndices(agg, aggIndices)
aggIndLen = length(aggIndices)
base = GenericBasis([aggIndLen])
FCFact = getFranckCondonFactors(agg, aggIndices)
FCProd = getFCProd(agg, FCFact, aggIndices, vibindices)
Ham = getAggHamiltonian(agg, aggIndices, FCFact)
basis = GenericBasis([length(aggIndices)])

Ham_bath = getAggHamiltonianBath(agg)
Ham_sys = getAggHamiltonianSystem(agg; groundEnergy = false)
b = GenericBasis([aggIndLen])
b_sys = GenericBasis([size(Ham_sys, 1)])
b_bath = GenericBasis([size(Ham_bath, 1)])

Ham_int = getAggHamiltonianInteraction(agg, aggIndices, FCFact)
Ham_S = Ham - Ham_int
Ham_B = tensor(Ham_bath, OneDenseOperator(b_sys))
E0 = Ham_B.data[1, 1]
for I = 1:aggIndLen
    Ham_B.data[I, I] -= E0
end
Ham_B = DenseOperator(b, b, Ham_B.data)

H_lambda, H_S = eigen(Ham_S.data)
H_Sinv = inv(H_S)
H_lambda = diagm(H_lambda)

t_max = 0.01
t_count = 200
t0 = 0.0
t_step = (t_max - t0) / (t_count)
tspan = [t0:t_step:t_max;]
T = 300
mu_array = [[2, 1]]
W0_1 = thermal_state(T, mu_array, Ham, vibindices, aggIndices; diagonalize = false)
mu_array = [[1, 2]]
W0_2 = thermal_state(T, mu_array, Ham, vibindices, aggIndices; diagonalize = false)
W0 = 0.8 * W0_1 + 0.2 * W0_2
W0 = DenseOperator(W0.basis_l, W0.basis_r, complex(W0.data))
W0_bath = get_rho_bath(W0, agg, FCProd, aggIndices, vibindices)
W0_bath = DenseOperator(W0_bath.basis_l, W0_bath.basis_r, complex(W0_bath.data))
rho0 = trace_bath(W0, agg, FCProd, aggIndices, vibindices)
rho0 = DenseOperator(rho0.basis_l, rho0.basis_r, complex(rho0.data))

### trace_bath
t2_0 = @benchmark begin
    rho0 = trace_bath(W0, agg, FCProd, aggIndices, vibindices)
end

### evolutionExact
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 10.0
t2_1 = @benchmark begin
    op_array = evolutionExact(W0, tspan, Ham)
end

### evolutionApproximate
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 2.0
t2_2 = @benchmark begin
    op_array = evolutionApproximate(W0, tspan, Ham)
end

### schroedinger
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 1.0
psi0 = randstate(basis)
t2_3 = @benchmark begin
    Tspan, psi_t = schroedinger(
        psi0,
        tspan,
        Ham;
        reltol = 1e-4,
        abstol = 1e-4,
        alg = OrdinaryDiffEq.Tsit5(),
    )
end

### schroedinger
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 10.0
t2_4 = @benchmark begin
    Tspan, rho_t = liouvilleVonNeumann(
        W0,
        tspan,
        Ham;
        reltol = 1e-4,
        abstol = 1e-4,
        alg = OrdinaryDiffEq.Tsit5(),
    )
end

### evolutionOperatorIterator
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 10.0
elLen = length(agg.molecules)
rho_t_exact = zeros(ComplexF64, length(tspan), elLen + 1, elLen + 1)
t2_5 = @benchmark begin
    t_i = 0
    foreach(
        evolutionOperatorIterator(Ham, tspan; diagonalize = false, approximate = true),
    ) do U_op
        t_i = t_i + 1
        # println(t_i / 150. * 100)
        W = U_op * W0 * U_op'
        rho_traced = trace_bath(W.data, agg, FCProd, aggIndices, vibindices)
        rho_t_exact[t_i, :, :] = rho_traced
    end
end

### master_int
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 20.0
p = (
    Ham_S,
    Ham_int,
    H_lambda,
    H_S,
    H_Sinv,
    W0,
    W0_bath,
    agg,
    FCProd,
    aggIndices,
    vibindices,
    ComplexF64,
)
t2_6 = @benchmark begin
    Tspan, W_t_master_int = master_int(
        W0,
        tspan,
        Ham_S,
        Ham_int;
        reltol = 1.0e-4,
        abstol = 1.0e-4,
        alg = DelayDiffEq.MethodOfSteps(DelayDiffEq.Tsit5()),
    )

    elLen = length(agg.molecules)
    rho_t_master = zeros(ComplexF64, length(tspan), elLen + 1, elLen + 1)
    for t_i = 1:length(tspan)
        t = tspan[t_i]
        U_op_S = evolutionOperator(Ham_sys, t)
        rho = trace_bath(W_t_master_int[t_i], agg, FCProd, aggIndices, vibindices)
        rho = U_op_S * rho * U_op_S'
        rho_t_master[t_i, :, :] = rho.data
    end
end

### evolutionOperatorIterator
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 20.0
p = (
    Ham_S,
    Ham_int,
    H_lambda,
    H_S,
    H_Sinv,
    Ham_B,
    W0,
    W0_bath,
    agg,
    FCProd,
    aggIndices,
    vibindices,
    ComplexF64,
)
t2_7 = @benchmark begin
    Tspan, rho_t_ansatz_int = master_ansatz(
        rho0,
        tspan,
        p;
        reltol = 1.0e-4,
        abstol = 1.0e-4,
        alg = DelayDiffEq.MethodOfSteps(DelayDiffEq.Tsit5()),
    )

    elLen = length(agg.molecules)
    rho_t_ansatz = zeros(ComplexF64, length(tspan), elLen + 1, elLen + 1)
    for t_i = 1:length(tspan)
        t = tspan[t_i]
        U_op_S = evolutionOperator(Ham_sys, t)
        rho = U_op_S * rho_t_ansatz_int[t_i] * U_op_S'
        rho_t_ansatz[t_i, :, :] = rho.data
    end
end

dic_commit["trace_bath"] = t2_0.times
dic_commit["evolutionExact"] = t2_1.times
dic_commit["evolutionApproximate"] = t2_2.times
dic_commit["schroedinger"] = t2_3.times
dic_commit["liouvilleVonNeumann"] = t2_4.times
dic_commit["evolutionOperatorIterator"] = t2_5.times
dic_commit["master_int"] = t2_6.times
dic_commit["master_ansatz"] = t2_7.times

dic_commits[commit_hash] = dic_commit
json_string = JSON.json(dic_commits)

open("benchmark/benchmark_commits.json", "w") do f
    write(f, json_string)
end

max_stars_repo_name: detrin/OpenQuantumSystems

max_issues_repo_name: detrin/OpenQuantumSystems

max_forks_repo_name: detrin/OpenQuantumSystems


text:
using Manifolds, Manopt, BenchmarkTools, Test

include("../test/solvers/trust_region_model.jl")

n = size(A, 1)
p = 2
N = Grassmann(n, p)
M = PowerManifold(N, ArrayPowerRepresentation(), 2)
x = random_point(M)

x_opt = trust_regions(M, cost, rgrad, rhess, x; max_trust_region_radius=8.0)
@btime x_opt = trust_regions($M, $cost, $rgrad, $rhess, $x; max_trust_region_radius=8.0)

h = RHess(M, A, p)
g = RGrad(M, A)
x_opt2 = trust_regions(
    M, cost, g, h, x; max_trust_region_radius=8.0, evaluation=MutatingEvaluation()
)

@btime trust_regions(
    $M, $cost, $g, $h, x2; max_trust_region_radius=8.0, evaluation=$(MutatingEvaluation())
) setup = (x2 = deepcopy($x))

x3 = deepcopy(x)
trust_regions!(
    M, cost, g, h, x3; max_trust_region_radius=8.0, evaluation=MutatingEvaluation()
)

@btime trust_regions!(
    $M, $cost, $g, $h, x3; max_trust_region_radius=8.0, evaluation=$(MutatingEvaluation())
) setup = (x3 = deepcopy($x))

@test distance(M, x_opt, x_opt2) ≈ 0
@test distance(M, x_opt, x3) ≈ 0