text:
#-------------------------------------------------------------------------------
# Set of useful functions related to LinearModel struct.
#-------------------------------------------------------------------------------

export LinearModel, trsbox!

"""

Define the type for linear models.

    - 'n': model dimension.
    - 'm': number of interpolation points.
    - 'imin': index of the function interpolated.
    - 'kopt': index of the best point so far in 'Y'.
    - 'c': constant of the model.
    - 'g': gradient of the model.
    - 'xbase': origin of the sample set and center of the model.
    - 'xopt': best point so far (in terms of function values).
    - 'fval': set of the function values of the interpolation points.
    - 'dst': distances between 'xopt' and other interpolation points.
    - 'τY': auxiliary vector of QR factorization.
    - 'jpvtY': auxiliary vector with the permutations of QR factorization.
    - 'factorsY': auxiliary matrix of QR factorization.
    - 'Y': set of interpolation points, shifted from the center of the sample set 'xbase'.

"""
struct LinearModel <: AbstractModel

    n        :: Int64
    m        :: Int64
    imin     :: Ref{Int64}
    kopt     :: Ref{Int64}
    c        :: Ref{Float64}
    g        :: AbstractVector
    xbase    :: AbstractVector
    xopt     :: AbstractVector
    fval     :: AbstractVector
    dst      :: AbstractVector
    τY       :: AbstractVector
    jpvtY    :: AbstractVector
    factorsY :: AbstractMatrix
    Y        :: AbstractMatrix

end

"""

    LinearModel(n::Int64)

Create a linear interpolation model for a function (not specified yet) of R^{n}.

    - 'n': dimension of the model.

Returns a model of LinearModel type.

"""
LinearModel(n::Int64) = LinearModel(
    n, n + 1, Ref{Int64}(), Ref{Int64}(), Ref{Float64}(),
    zeros(Float64, n), zeros(Float64, n), zeros(Float64, n), zeros(Float64, n + 1),
    zeros(Float64, n), zeros(Float64, n), zeros(Int64, n), zeros(Float64, n, n), zeros(Float64, n, n))

( model::LinearModel )( x::AbstractVector ) = model.c[] + dot( model.g, x )

"""

    rebuild_model!(model::LinearModel)

Computes the QR factorization of the matrix 'Y', the gradient and the constant of the model.

The function modifies the argument:

    - 'model': model of LinearModel type.

"""
function rebuild_model!(
                        model::LinearModel
                        )

    # Computes the QR-factorization of the matrix Y and stores the information in 'model.factorsY', 'model.τY' and 'model.jpvtY'.
    @. model.factorsY = model.Y
    qrY = qr!(model.factorsY, Val(true))
    @. model.τY = qrY.τ
    @. model.jpvtY = qrY.jpvt

    # Solves the system to determine the gradient of the linear model 'model.g'.
    @. model.g = model.fval[2:end] - model.fval[1]
    ldiv!(qrY, model.g)

    # Determines the constant of the linear model 'model.c'
    model.c[] = model.fval[1] - dot(model.g, model.xbase)

end

"""

    construct_model!(func_list::Array{Function, 1}, imin_idx::Int64,
                        δ::Float64, fbase::Float64, xbase::Vector{Float64},
                        ao::Vector{Float64}, bo::Vector{Float64},
                        model::LinearModel)

Constructs the model based in the given information.

    - 'func_list': list containing the functions that determine the objective
    function fmin.

    - 'imin_idx': index belonging to the set I_{min}('xbase').

    - 'δ': radius of the sample set.

    - 'fbase': objective function value in 'xbase'.

    - 'xbase': n-dimensional vector (origin of the sample set).

    - 'ao': n-dimensional vector with the shifted lower bounds.

    - 'bo': n-dimensional vector with the shifted upper bounds.

The function modifies the argument:

    - 'model': model of LinearModel type.

"""
function construct_model!(
                            func_list::Array{Function, 1},
                            imin_idx::Int64,
                            δ::Float64,
                            fbase::Float64,
                            xbase::Vector{Float64},
                            ao::Vector{Float64},
                            bo::Vector{Float64},
                            model::LinearModel
                            )

    # Saves information about the index i ∈ I_{min}('xbase') and f_{i}('xbase').
    model.imin[] = imin_idx
    model.fval[1] = fbase

    kopt = 0

    for i = 1:model.n

        # Saves information about 'model.xbase' and 'model.xopt'.
        model.xbase[i] = xbase[i]
        model.xopt[i] = xbase[i]

        # Constructs the set interpolation points 'model.Y', shifted from 'xbase'.
        # Evaluates the function values and saves the information in 'model.fval'.
        if bo[i] == 0.0

            model.Y[i, i] = - δ
            xbase[i] -= δ
            model.fval[i + 1] = fi_eval(func_list, imin_idx, xbase)
            xbase[i] += δ

        else

            model.Y[i, i] = δ
            xbase[i] += δ
            model.fval[i + 1] = fi_eval(func_list, imin_idx, xbase)
            xbase[i] -= δ

        end

        # Saves the distance between the i-th point in 'model.Y' and 'model.xbase'.
        model.dst[i] = δ

        # Searches for the least function value to determine 'kopt'.
        if model.fval[i + 1] < fbase

            kopt = i
            fbase = model.fval[i + 1]

        end

    end

    # Saves the index of the best point in 'model.kopt'
    model.kopt[] = kopt

    # If the best point is not 'xbase', updates the vector 'xopt' and saves the information in 'model.xopt'.
    # Also updates the vector of 'model.dst', to store the distance between the interpolation points and 'xopt'.
    # Note that in the 'kopt' position of the 'model.dst' vector we keep the distance between 'xbase' and 'xopt' points.
    if kopt != 0

        model.xopt[kopt] += model.Y[kopt, kopt]
        @. model.dst *= sqrt(2.0)
        model.dst[kopt] = δ

    end

    # Computes the constant and gradient of the linear model and the QR-factorization of the system.
    rebuild_model!( model )

end

"""

    update_model!(t::Int64, fnew::Float64, xnew::Vector{Float64},
                    model::LinearModel; trsbox_step::Bool=false)

Updates the model based in the given information.

    - 't': index of the point that must leave the sample set.

    - 'fnew: funtion value of the nw point.

    - 'xnew': n-dimensional vector (new sample point).

The function modifies the argument:

    - 'model': model of LinearModel type.

The optional argument is:

    - 'trsbox_step': Indicates if the step is of trust-region type. Set to 'false' by default.

"""
function update_model!(
                        t::Int64,
                        fnew::Float64,
                        xnew::Vector{Float64},
                        model::LinearModel;
                        trsbox_step::Bool=false
                        )

    if trsbox_step

        # If the step is of the TRSBOX type, the center of the sample set 'xbase' is never changed. Note that 't' is different from 'model.kopt[]'.
        # Since such a step is a descent direction for the model, 'fnew' is always less than the value of the function in 'xopt'.
        # The information about 'xopt' is also useful for the new model, so we maintain 'xopt' in the interpolation set 'Y', overwrite the point 't' by 'xnew',
        # and update 'xopt' to 'xnew'.

        model.kopt[] = t
        model.fval[t + 1] = fnew

        for i=1:model.n

            model.Y[t, i] = xnew[i] - model.xbase[i]
            model.xopt[i] = xnew[i]

        end

        @views for i=1:model.n

            if i == t

                model.dst[i] = norm( model.Y[i, :] )

            else

                model.dst[i] = norm( model.Y[i, :] - model.Y[t, :] )

            end

        end

    else

        # On the other hand, if the step is of type ALTMOV, the center of the sample set 'xbase' can be dropped or not from 'Y'.
        # Also, the best point so far 'xopt' can be changed or not, although the information about 'xopt' is always used for the new model.
        # Note that in the case where the points 'xbase' and 'xopt' coincide, the index 't' will never be 0.

        if t == 0

            @assert model.kopt[] != 0 "In ALTMOV iterations, if the index 't' is set to 0, so the index 'kopt' cannot be 0."

            # The point 'xbase' is dropped.

            if fnew < model.fval[ model.kopt[] + 1 ]

                # The point 'xnew' is set to the new 'xbase' and 'xopt' is updated to 'xnew'.
                # The old information about the point 'xopt' is used in the new model.

                for i=1:model.n

                    @. model.Y[i, :] += model.xbase - xnew

                end

                for i=1:model.n

                    model.xbase[i] = xnew[i]
                    model.xopt[i] = xnew[i]
                    model.dst[i] = norm( model.Y[i, :] )

                end

                model.fval[1] = fnew
                model.kopt[] = 0

            else

                # The point 'xopt' is set to be the new 'xbase' and 'xnew' is incorporated as a point in the set 'Y'.

                for i=1:model.n

                    if i == model.kopt[]

                        @. model.Y[i, :] = xnew - model.xopt

                    else

                        @. model.Y[i, :] += model.xbase - model.xopt

                    end

                end

                for i=1:model.n

                    model.xbase[i] = model.xopt[i]
                    model.dst[i] = norm( model.Y[i, :] )

                end

                model.fval[1] = model.fval[ model.kopt[] + 1 ]
                model.fval[ model.kopt[] + 1 ] = fnew
                model.kopt[] = 0

            end

        else

            # The point 'xbase' is maintained.

            model.fval[ t + 1 ] = fnew
            @. model.Y[t, :] = xnew - model.xbase

            if fnew < model.fval[ model.kopt[] + 1 ]

                # The point 'xopt' is updated to 'xnew'.

                model.dst[t] = norm( model.Y[t, :] )
                model.kopt[] = t

                for i=1:model.n

                    model.xopt[i] = xnew[i]

                    if i != t

                        model.dst[i] = norm( model.Y[i, :] - model.Y[t, :] )

                    end

                end

            else

                model.dst[t] = norm( xnew - model.xopt )

            end

        end

    end

    # Computes the constant and gradient of the linear model and the QR-factorization of the system.
    rebuild_model!( model )

end

"""

    construct_new_model!(func_list::Array{Function, 1}, imin_idx::Int64, δ::Float64,
                            fbase::Float64, xbase::Vector{Float64}, a::Vector{Float64},
                            b::Vector{Float64}, ao::Vector{Float64}, bo::Vector{Float64},
                            model::LinearModel)

Constructs the model based in the given information.

    - 'func_list': list containing the functions that determine the objective
    function fmin.

    - 'imin_idx': index belonging to the set I_{min}('xbase').

    - 'δ': radius of the sample set.

    - 'fbase': objective function value in 'xbase'.

    - 'xbase': n-dimensional vector (origin of the sample set).

    - 'a': n-dimensional vector with the lower bounds.

    - 'b': n-dimensional vector with the upper bounds.

The function modifies the argument:

    - 'ao': n-dimensional vector with the shifted lower bounds.

    - 'bo': n-dimensional vector with the shifted upper bounds.

    - 'model': model of LinearModel type.

"""
function construct_new_model!(
                            func_list::Array{Function, 1},
                            imin_idx::Int64,
                            δ::Float64,
                            fbase::Float64,
                            xbase::Vector{Float64},
                            a::Vector{Float64},
                            b::Vector{Float64},
                            ao::Vector{Float64},
                            bo::Vector{Float64},
                            model::LinearModel
                            )

    # Sets Y to a empty matrix.
    @. model.Y = 0.0

    # Saves information about the index i ∈ I_{min}('xbase') and f_{i}('xbase').
    model.imin[] = imin_idx
    model.fval[1] = fbase

    kopt = 0

    for i = 1:model.n

        # Saves information about 'model.xbase' and 'model.xopt'.
        model.xbase[i] = xbase[i]
        model.xopt[i] = xbase[i]

        # Computes the new relative bounds 'ao' and 'bo'.
        ao[i] = a[i] - xbase[i]
        bo[i] = b[i] - xbase[i]

        # Computes the new interpolation points 'model.Y', shifted from 'xbase'.
        # Evaluates the function values and saves the information in 'model.fval'.
        if bo[i] > δ

            model.Y[i, i] += δ
            xbase[i] += δ
            model.fval[i + 1] = fi_eval(func_list, imin_idx, xbase)
            xbase[i] -= δ

        else

            model.Y[i, i] -= δ
            xbase[i] -= δ
            model.fval[i + 1] = fi_eval(func_list, imin_idx, xbase)
            xbase[i] += δ

        end

        # Saves the distance between the i-th point in 'model.Y' and 'model.xbase'.
        model.dst[i] = δ

        # Searches for the least function value to determine 'kopt'.
        if model.fval[i + 1] < fbase

            kopt = i
            fbase = model.fval[kopt]

        end

    end

    # Saves the index of the best point in 'model.kopt'
    model.kopt[] = kopt

    # If the best point is not 'xbase', updates the vector 'xopt' and saves the information in 'model.xopt'.
    # Also updates the vector of 'model.dst', to store the distance between the interpolation points and 'xopt'.
    # Note that in the 'kopt' position of the 'model.dst' vector we keep the distance between 'xbase' and 'xopt' points.
    if kopt != 0

        model.xopt[ kopt ] += model.Y[ kopt, kopt ]
        @. model.dst *= sqrt(2.0)
        model.dst[kopt] = δ

    end

    # Computes the constant and gradient of the linear model and the QR-factorization of the system.
    rebuild_model!( model )

end

"""

    stationarity_measure(model::LinearModel, a::Vector{Float64},
                            b::Vector{Float64})

Calculates the stationarity measure π_{k} = || P_{Ω}( x_{k} - g_{k} ) - x_{k} ||.

    - 'model': model of LinearModel type.

    - 'a': n-dimensional vector with the lower bounds.

    - 'b': n-dimensional vector with the upper bounds.

Returns the stationarity measure.

"""
function stationarity_measure(
                                model::LinearModel,
                                a::Vector{Float64},
                                b::Vector{Float64}
                                )

    aux = 0.0
    for i = 1:model.n

        aux += ( min( max( a[i], model.xopt[i] - model.g[i] ), b[i] ) - model.xopt[i] ) ^ 2.0

    end

    return sqrt(aux)

end

"""

    compute_active_set!(model::LinearModel, a::Vector{Float64},
                        b::Vector{Float64}, active_set::Vector{Bool})

Computes the set of active constraints.

    - 'model': model of LinearModel type.

    - 'a': n-dimensional vector with the lower bounds.

    - 'b': n-dimensional vector with the upper bounds.

The function modifies the argument:

    - 'active_set': boolean n-dimensional vector with the active constraints.

"""
function compute_active_set!(
                        model::LinearModel,
                        a::Vector{Float64},
                        b::Vector{Float64},
                        active_set::Vector{Bool}
                        )

    for i=1:model.n

        if ( model.xopt[i] == a[i] && model.g[i] ≥ 0.0 ) || ( model.xopt[i] == b[i] && model.g[i] ≤ 0.0 )

            active_set[i] = true

        else

            active_set[i] = false

        end

    end

end

"""

    trsbox!(model::LinearModel, Δ::Float64, a::Vector{Float64}, b::Vector{Float64},
                active_set::Vector{Bool}, x::Vector{Float64}, d::Vector{Float64},
                s::Vector{Float64})

Computes a descent direction for the model based on TRSBOX routine of BOBYQA.

    - 'model': model of LinearModel type.

    - 'Δ': trust-region radius.

    - 'a': n-dimensional vector with the lower bounds.

    - 'b': n-dimensional vector with the upper bounds.

The function modifies the argument:

    - 'active_set': boolean n-dimensional vector with the active constraints.

    - 'x': n-dimensional vector (candidate for new sampling point 'x': 'xopt' + 'd').

    - 'd': n-dimensional vector (descent direction).

    - 's': n-dimensional vector (auxiliary vector for workspace).

Returns a flag indicating the stopping criteria.

"""
function trsbox!(
                    model::LinearModel,
                    Δ::Float64,
                    a::Vector{Float64},
                    b::Vector{Float64},
                    active_set::Vector{Bool},
                    x::Vector{Float64},
                    d::Vector{Float64},
                    s::Vector{Float64}
                    )

    # Copies vector 'xopt' to 'x' and sets vectors 'd' and 's' to zero.
    for i = 1:model.n

        x[i] = model.xopt[i]
        d[i] = 0.0
        s[i] = 0.0

    end

    # Saves the best function value so far.
    fopt = model.fval[ model.kopt[] + 1 ]

    # Creates the active constraint vector 'active_set'.
    compute_active_set!(model, a, b, active_set)

    # If the set of active constraints is full, 'd' is returned as the null vector.
    if sum(active_set) == model.n

        return :full_active_set

    end

    # Computes the symmetric vector of the projection of the gradient of the model
    # by the set of active constraints and stores in 's'.
    projection_active_set!(model.g, active_set, s, sym = true)

    # If 's' is a null vector, we won't get any decreases in that direction.
    # Therefore, 'd' is returned as the null vector.
    if norm(s) ≈ 0.0

        return :null_search_direction

    end

    while true

        # Computes the step 'α' along the direction 's'.
        α, chose_αΔ, idx_B = compute_alpha_linear(Δ, a, b, x, d, s)

        # Updates the search direction 'd'.
        @. d += α * s

        if chose_αΔ

            # If the trust-region boundary has been reached (case 'α' = 'α_Δ'),
            # prepares for an auxiliary procedure that tries to increase the reduction
            # of the obtained model value.

            while true

                # Uses 'x' as work vector to store projection of direction 'd' by a set of active constraints 'active_set'.
                projection_active_set!(d, active_set, x)

                # Computes P_{I}(∇model(xopt+d)) = P_{I}('model.g') and stores in 's'.
                projection_active_set!(model.g, active_set, s)

                pdTpd = dot( x, x )
                pdTpg = dot( x, s )
                pgTpg = dot( s, s )

                # Tests the stopping criteria.
                if ( pdTpd * pgTpg - pdTpg ^ 2.0 ) ≤ 1.0e-4 * ( fopt - model.c[] - dot( model.g, model.xopt ) - dot( model.g, d ) )

                    step_projection!(model, a, b, x, d)

                    return :trust_region_boundary

                else

                    # Otherwise, computes a new search direction 's', spanned by P_{I}(d) and P_{I}(∇model(xopt+d)),
                    # that satisfies ||s|| = ||P_{I}(d)||, s'P_{I}(d) = 0, and s'P_{I}(∇model(xopt+d)) < 0.0.
                    # Note that 'x' stores P_{I}(d) and 's' stores P_{I}(∇model(xopt+d)).
                    # After the call of 'new_search_direction', the desired direction is allocated in 's'.
                    new_search_direction!(pdTpd, pdTpg, pgTpg, x, s)

                    # Verifies if the direction 's' is null.
                    if norm(s) ≈ 0.0

                        step_projection!(model, a, b, x, d)

                        return :null_new_search_direction

                    end

                    # Computes the parameter θ that satisfies a ≤ xopt + d(θ) ≤ b, and g'd(θ) < 0.0.
                    # Note that the vector 's' holds d(θ) - d.
                    chose_θQ = compute_theta_linear!(model, a, b, d, x, s)

                    if norm(s) ≈ 0.0

                        step_projection!(model, a, b, x, d)

                        return :null_angle_search_direction

                    end

                    # Updates the direction 'd' to be d(θ).
                    @. d += s

                    # Updates 'index_set' with indices for which xopt + d is restricted by one of the bounds on the variables.
                    update_active_set!( a, b, model.xopt, d, active_set )

                    if chose_θQ

                        # Computes P_{I}(∇model(xopt+d)) = P_{I}('model.g') and stores in 's'.
                        projection_active_set!(model.g, active_set, s)

                        # Verifies the stopping criteria.
                        if ( norm(s) * Δ ) ≤ 1.0e-2 * ( fopt - model.c[] - dot( model.g, model.xopt ) - dot( model.g, d ) )

                            step_projection!(model, a, b, x, d)

                            return :tired_trust_region

                        end

                    end

                end

            end

        else

            # Otherwise, the current line search is restricted by a bound constraint. If the stooping criteria is true, the computational effort
            # to perform a new iteration is not promising and the procedure is terminated.

            # Computes - P_{I}(∇model(xopt+d)) = - P_{I}('model.g') and stores in 's'.
            projection_active_set!(model.g, active_set, s, sym = true)

            # Verifies the stopping criteria.
            if ( norm(s) * Δ ) ≤ 1.0e-2 * ( fopt - model.c[] - dot( model.g, x ) - dot( model.g, d ) )

                step_projection!(model, a, b, x, d)

                return :tired_bound_constraint

            else

                # Adds the index of the achieved bound constraint to 'active_set',
                # and updates 's' to be the symmetric projection of 'model.g' in the new set.

                if idx_B != 0

                    active_set[ idx_B ] = true
                    s[ idx_B ] = 0.0

                end

            end

        end

        # If the set of active constraints is full, the procedure terminates.
        if sum(active_set) == model.n

            step_projection!(model, a, b, x, d)

            return :full_active_set

        end

        # If the direction 's' is null, the procedure also terminates.
        if norm(s) ≈ 0.0

            step_projection!(model, a, b, x, d)

            return :null_search_direction

        end

    end

end

"""

    choose_index_trsbox(model::LinearModel, xnew::Vector{Float64}, v::Vector{Float64})

Chooses the point that must leave the sampling set for trust-region steps.

    - 'model': model of LinearModel type.

    - 'xnew': n-dimensional vector (new sample point).

The function modifies the argument:

    - 's': n-dimensional vector (auxiliary vector for workspace).

Returns an index 't' which indicates the point that must leave the sample set.

"""
function choose_index_trsbox(
                                model::LinearModel,
                                xnew::Vector{Float64},
                                v::Vector{Float64}
                                )

    qrY = QRPivoted( model.factorsY, model.τY, model.jpvtY )
    α = - Inf
    t = - 1

    for i = 1:model.n

        @. v = 0.0
        v[i] = 1.0
        ldiv!(qrY, v)
        α_t = abs( 1.0 + dot( xnew, v ) - dot( model.xbase, v ) - dot( model.Y[i, :], v ) )

        if α_t > α

            α = α_t
            t = i

        end

    end

    return t

end

"""

    lagrange_pol_t(model::LinearModel, idx_t::Int64, ∇Λ_t::Vector{Float64}, y::Vector{Float64})

Computes the t-th Lagrange polynomial at 'y', i.e., computes Λ_{t}(y).

    - 'model': model of LinearModel type.

    - 'idx_t': index of the Lagrange polynomial.

    - '∇Λ_t': n-dimensional vector (gradient of the t-th Lagrange polynomial).

    - 'y': n-dimensional vector (point of interest).

Returns the function value Λ_{t}(y).

"""
function lagrange_pol_t(
                        model::LinearModel,
                        idx_t::Int64,
                        ∇Λ_t::Vector{Float64},
                        y::Vector{Float64}
                        )

    if idx_t == 0

        return 1.0 + dot( y, ∇Λ_t ) - dot( model.xbase, ∇Λ_t )

    else

        return dot( y, ∇Λ_t ) - dot( model.xbase, ∇Λ_t )

    end

end

"""

    altmov!(model::LinearModel, idx_t::Int64, Δ::Float64, a::Vector{Float64},
                b::Vector{Float64}, x::Vector{Float64}, d::Vector{Float64},
                v::Vector{Float64}, w::Vector{Float64}, altmov_set::Vector{Bool})

Calculates a direction to improve the poisedness of the sample set based on ALTMOV routine of BOBYQA.

    - 'model': model of LinearModel type.

    - 'idx_t': index of the point that must leave the sample set.

    - 'Δ': trust-region radius.

    - 'a': n-dimensional vector with the lower bounds.

    - 'b': n-dimensional vector with the upper bounds.

The function modifies the argument:

    - 'x': n-dimensional vector (candidate for new sampling point 'x': 'xopt' + 'd').

    - 'd': n-dimensional vector (descent direction).

    - 'v': n-dimensional vector (auxiliary vector for workspace).

    - 'w': n-dimensional vector (auxiliary vector for workspace).

    - 'altmov_set': boolean n-dimensional vector with the active constraints.

Returns a flag indicating the stopping criteria.

"""
function altmov!(
                    model::LinearModel,
                    idx_t::Int64,
                    Δ::Float64,
                    a::Vector{Float64},
                    b::Vector{Float64},
                    x::Vector{Float64},
                    d::Vector{Float64},
                    v::Vector{Float64},
                    w::Vector{Float64},
                    altmov_set::Vector{Bool}
                    )

    # Reconstruct the QR-factorization of the system.
    qrY = QRPivoted( model.factorsY, model.τY, model.jpvtY )

    # Initializes some usefull variables and constants
    best_abs_ϕ = - 1.0
    best_idx = 0
    best_α = Inf
    a_0 = - Δ ^ 2.0
    a_1 = 0.0

    # Computes the gradient of the "idx_t"-th Lagrange polynomial and stores the information in vector 'v'.
    if idx_t == 0.0

        @. x = - 1.0
        ldiv!( v, qrY, x )

    else

        @. x = 0.0
        x[idx_t] = 1.0
        ldiv!( v, qrY, x )

    end

    # Calculates the "Usual" Alternative Step
    for j = 1:model.n

        # Uses 'd' as a workspace to save the difference between the j-th interpolation point and 'xopt'.
        # If 'xopt' and 'xbase' coincide, such differences are already calculated and available in the matrix 'model.Y'.
        # Otherwise, if 'xopt' is not the center of the sample set, and 'j' = 'kopt', 'd' holds the difference between 'xbase' and 'xopt'.
        if model.kopt[] == 0

            @. d = model.Y[ j, : ]

        else

            if j == model.kopt[]

                @. d = - model.Y[ model.kopt[], : ]

            else

                @. d = model.Y[ j, : ] - model.Y[ model.kopt[], : ]

            end

        end

        # Gets the bounds for α_{j} relative to the trust-region.
        roots = []
        a_2 = dot( d, d )
        solve_quadratic!(a_2, a_1, a_0, roots)

        α_upper = maximum(roots)

        if ( α_upper == NaN )

            for i = 1:model.n

                x[i] = model.xopt[i]
                d[i] = 0.0

            end

            return :invalid_steplenght

        else

            α_lower = minimum(roots)

        end

        # Adjusts the bounds to the box constraints.
        for i = 1:model.n

            if d[i] > 0.0

                α_lower = max( α_lower, ( a[i] - model.xopt[i] ) / d[i] )
                α_upper = min( α_upper, ( b[i] - model.xopt[i] ) / d[i] )

            elseif d[i] < 0.0

                α_lower = max( α_lower, ( b[i] - model.xopt[i] ) / d[i] )
                α_upper = min( α_upper, ( a[i] - model.xopt[i] ) / d[i] )

            end

        end

        # Since the Lagrange polynomial is linear, the optimal α is either lower or upper value.
        # Uses 'x' as a workspace to hold 'xopt + α * (y_{j} - xopt)'
        @. x = model.xopt + α_lower * d
        ϕ_lower = lagrange_pol_t( model, idx_t, v, x )

        @. x = model.xopt + α_upper * d
        ϕ_upper = lagrange_pol_t( model, idx_t, v, x )

        if abs( ϕ_lower ) > abs( ϕ_upper )

            α_j = α_lower
            abs_ϕ_j = abs( ϕ_lower )

        else

            α_j = α_upper
            abs_ϕ_j = abs( ϕ_upper )

        end

        # Compares the best ϕ_j value so far
        if abs_ϕ_j > best_abs_ϕ

            best_abs_ϕ = abs_ϕ_j
            best_idx = j
            best_α = α_j

        end

    end

    # Calculates the "Cauchy" Alternative Step

    Λ_1 = altmov_cauchy!(model, idx_t, Δ, a, b, v, d, w, altmov_set)
    Λ_2 = altmov_cauchy!(model, idx_t, Δ, a, b, v, x, w, altmov_set, sym = true)

    # Chooses the best step
    ( val, idx ) = findmax( [ best_abs_ϕ, abs( Λ_1 ), abs( Λ_2 ) ] )

    if idx == 1

        if model.kopt[] == 0

            @. d = best_α * model.Y[ best_idx, : ]

        else

            if best_idx == model.kopt[]

                @. d = - best_α * model.Y[ best_idx, :]

            else

                @. d = best_α * ( model.Y[ best_idx, : ] - model.Y[ model.kopt[], : ] )

            end

        end

        step_projection!(model, a, b, x, d)

        return :usual_altmov

    elseif idx == 2

        step_projection!(model, a, b, x, d)

    else

        @. d = x
        step_projection!(model, a, b, x, d)

    end

    return :cauchy_altmov

end

"""

    altmov_cauchy!(model::LinearModel, idx_t::Int64, Δ::Float64, a::Vector{Float64},
                    b::Vector{Float64}, ∇Λ_t::Vector{Float64}, s::Vector{Float64},
                    z::Vector{Float64}, active_set::Vector{Bool}; sym::Bool=false)

Calculates a direction to improve the poisedness of the sample set based on the Cauchy step of ALTMOV routine from BOBYQA.

    - 'model': model of LinearModel type.

    - 'idx_t': index of the point that must leave the sample set.

    - 'Δ': trust-region radius.

    - 'a': n-dimensional vector with the lower bounds.

    - 'b': n-dimensional vector with the upper bounds.

    - '∇Λ_t': n-dimensional vector (gradient of the t-th Lagrange polynomial).

The function modifies the argument:

    - 's': n-dimensional vector (direction).

    - 'z': n-dimensional vector (candidate for new sampling point 'x': 'xopt' + 's').

    - 'active_set': boolean n-dimensional vector with the active constraints.

The optional argument is:

    - 'sym': indicates if it is calculated based on the symmetric direction of '∇Λ_t'. Set to 'false' by default.

Returns the function value of the t-th Lagrange polynomial in 'z', i.e., Λ_t(z).

"""
function altmov_cauchy!(
                        model::LinearModel,
                        idx_t::Int64,
                        Δ::Float64,
                        a::Vector{Float64},
                        b::Vector{Float64},
                        ∇Λ_t::Vector{Float64},
                        s::Vector{Float64},
                        z::Vector{Float64},
                        active_set::Vector{Bool};
                        sym::Bool=false
                        )

    if sym

        p = 1.0

    else

        p = 2.0

    end

    norm2_s = 0.0

    for i = 1:model.n

        if ( (- 1.0) ^ p ) * ∇Λ_t[i] > 0.0

            s[i] = a[i] - model.xopt[i]
            norm2_s += s[i] ^ 2.0

        elseif ( (- 1.0) ^ p ) * ∇Λ_t[i] < 0.0

            s[i] = b[i] - model.xopt[i]
            norm2_s += s[i] ^ 2.0

        else

            s[i] == 0.0

        end

    end

    if norm2_s ≤ ( Δ ^ 2.0 )

        @. z = model.xopt + s

        return lagrange_pol_t( model, idx_t, ∇Λ_t, z )

    end

    # Saves a copy of 's' in 'z'.
    copyto!( z, s )

    sum_free = 0.0
    sum_fixed = 0.0

    for i = 1:model.n

        if s[i] ≈ 0.0

            active_set[i] = true

        else

            sum_free += ∇Λ_t[i] ^ 2.0
            active_set[i] = false

        end

    end

    if sum_free ≈ 0.0

        s .= 0.0
        return 0.0

    end

    while true

        num_μ = Δ ^ 2.0 - sum_fixed

        if num_μ > 0.0

            μ = sqrt( num_μ / sum_free )

            stop = true
            sum_free = 0.0
            sum_fixed = 0.0

            for i = 1:model.n

                if active_set[i]

                    s[i] = z[i]
                    sum_fixed += z[i] ^ 2.0

                else

                    s[i] = ( (- 1.0) ^ ( p + 1 ) ) * μ * ∇Λ_t[i]

                    step_i = s[i] + model.xopt[i]

                    if ( step_i < a[i] ) || ( step_i > b[i] )

                        stop = false
                        active_set[i] = true
                        sum_fixed += z[i] ^ 2.0

                    else

                        sum_free += ∇Λ_t[i] ^ 2.0

                    end

                end

            end

            if stop || ( sum_free ≈ 0.0 )

                break

            end

        else

            break

        end

    end

    @. z = model.xopt + s

    return lagrange_pol_t( model, idx_t, ∇Λ_t, z )

end

"""
    compute_theta_linear!(model::LinearModel, a::Vector{Float64}, b::Vector{Float64},
                            d::Vector{Float64}, proj_d::Vector{Float64}, s::Vector{Float64})

Calculates the parameter θ that satisfies 'a' ≤ 'xopt' + 'd(θ)' ≤ 'b', and 'g'd(θ)' < 0.0, in 'trsbox!' routine.

    - 'model': model of LinearModel type.

    - 'a': n-dimensional vector with the lower bounds.

    - 'b': n-dimensional vector with the upper bounds.

    - 'd': n-dimensional vector (direction).

    - 'proj_d': n-dimensional vector (projection of the direction 'd' on the set of active constraints).

The function modifies the argument:

    - 's': n-dimensional vector (new direction).

Returns 'true' if the parameter 'θQ' is chosen, and 'false' otherwise.

"""
function compute_theta_linear!(
                                model::LinearModel,
                                a::Vector{Float64},
                                b::Vector{Float64},
                                d::Vector{Float64},
                                proj_d::Vector{Float64},
                                s::Vector{Float64}
                                )

    cond_B(θ) = cond_θB(θ, a, b, model.xopt, d, proj_d, s)

    θB = binary_search(0.0, 0.25 * π , cond_B, 1.0e-2)

    if θB == NaN

        @. s = 0.0

    else

        gTd = dot(model.g, model.g)
        gTpd = dot(model.g, proj_d)
        gTs = dot(model.g, s)

        cond_Q(θ) = cond_θQ_linear(θ, gTd, gTpd, gTs)

        θQ = binary_search(0.0, 0.25 * π , cond_Q, 1.0e-2)

        if θQ == NaN

            @. s = 0.0

        else

            θ = min(θB, θQ)

            @. s = - proj_d + cos(θ) * proj_d + sin(θ) * s

            if θ == θQ

                return true

            else

                return false

            end

        end

    end

end

max_stars_repo_name: aschwertner/LOWDER

max_issues_repo_name: aschwertner/LOWDER

max_forks_repo_name: aschwertner/LOWDER


text:
using MakieLayout
using AbstractPlotting
using Test

include("debugrect.jl")

# have BBoxes show as such because the default is very verbose
Base.show(io::IO, bb::BBox) = print(io, "BBox(l: $(left(bb)), r: $(right(bb)), b: $(bottom(bb)), t: $(top(bb)))")

@testset "GridLayout Zero Outside AlignMode" begin
    bbox = BBox(0, 1000, 0, 1000)
    layout = GridLayout(bbox = bbox, alignmode = Outside(0))
    dr = layout[1, 1] = DebugRect()

    @test dr.layoutnodes.computedbbox[] == bbox

    dr.topprot = 100
    @test dr.layoutnodes.computedbbox[] == BBox(0, 1000, 0, 900)
    dr.bottomprot = 100
    @test dr.layoutnodes.computedbbox[] == BBox(0, 1000, 100, 900)
    dr.leftprot = 100
    @test dr.layoutnodes.computedbbox[] == BBox(100, 1000, 100, 900)
    dr.rightprot = 100
    @test dr.layoutnodes.computedbbox[] == BBox(100, 900, 100, 900)

    dr2 = layout[1, 2] = DebugRect()
    @test layout.nrows == 1 && layout.ncols == 2
    colgap!(layout, 1, Fixed(0))

    @test dr.layoutnodes.computedbbox[].widths == dr2.layoutnodes.computedbbox[].widths == Float32[400.0, 800.0]
end

@testset "GridLayout Outside AlignMode" begin
    bbox = BBox(0, 1000, 0, 1000)
    layout = GridLayout(bbox = bbox, alignmode = Outside(100, 200, 50, 150))
    dr = layout[1, 1] = DebugRect()

    @test dr.layoutnodes.computedbbox[] == BBox(100, 800, 50, 850)

    dr.topprot = 100
    @test dr.layoutnodes.computedbbox[] == BBox(100, 800, 50, 750)
    dr.bottomprot = 100
    @test dr.layoutnodes.computedbbox[] == BBox(100, 800, 150, 750)
    dr.leftprot = 100
    @test dr.layoutnodes.computedbbox[] == BBox(200, 800, 150, 750)
    dr.rightprot = 100
    @test dr.layoutnodes.computedbbox[] == BBox(200, 700, 150, 750)
end

@testset "GridLayout Inside AlignMode" begin
    bbox = BBox(0, 1000, 0, 1000)
    layout = GridLayout(bbox = bbox, alignmode = Inside())
    dr = layout[1, 1] = DebugRect()

    @test dr.layoutnodes.computedbbox[] == BBox(0, 1000, 0, 1000)

    dr.topprot = 100
    @test dr.layoutnodes.computedbbox[] == BBox(0, 1000, 0, 1000)
    dr.bottomprot = 100
    @test dr.layoutnodes.computedbbox[] == BBox(0, 1000, 0, 1000)
    dr.leftprot = 100
    @test dr.layoutnodes.computedbbox[] == BBox(0, 1000, 0, 1000)
    dr.rightprot = 100
    @test dr.layoutnodes.computedbbox[] == BBox(0, 1000, 0, 1000)
end

@testset "GridLayout Mixed AlignMode" begin
    bbox = BBox(0, 1000, 0, 1000)
    layout = GridLayout(bbox = bbox, alignmode = Mixed(left = 0, top = 100))
    dr = layout[1, 1] = DebugRect()

    @test MakieLayout.protrusion(layout, Left()) == 0
    @test MakieLayout.protrusion(layout, Right()) == 0
    @test MakieLayout.protrusion(layout, Bottom()) == 0
    @test MakieLayout.protrusion(layout, Top()) == 0

    @test dr.layoutnodes.computedbbox[] == BBox(0, 1000, 0, 900)

    dr.topprot = 100
    @test MakieLayout.protrusion(layout, Top()) == 0
    @test dr.layoutnodes.computedbbox[] == BBox(0, 1000, 0, 800)

    dr.bottomprot = 100
    @test MakieLayout.protrusion(layout, Bottom()) == 100
    @test dr.layoutnodes.computedbbox[] == BBox(0, 1000, 0, 800)

    dr.leftprot = 100
    @test MakieLayout.protrusion(layout, Left()) == 0
    @test dr.layoutnodes.computedbbox[] == BBox(100, 1000, 0, 800)

    dr.rightprot = 100
    @test MakieLayout.protrusion(layout, Right()) == 100
    @test dr.layoutnodes.computedbbox[] == BBox(100, 1000, 0, 800)
end

@testset "Layoutables constructors" begin
    scene, layout = layoutscene()
    ax = layout[1, 1] = LAxis(scene)
    cb = layout[1, 2] = LColorbar(scene)
    gl2 = layout[2, :] = GridLayout()
    bu = gl2[1, 1] = LButton(scene)
    sl = gl2[1, 2] = LSlider(scene)

    scat = scatter!(ax, rand(10))
    le = gl2[1, 3] = LLegend(scene, [scat], ["scatter"])

    to = gl2[1, 4] = LToggle(scene)
    te = layout[0, :] = LText(scene, "A super title")
end

ext: jl


text:
module StructJuMP

export StructuredModel

using MathOptInterface
const MOI = MathOptInterface

using JuMP # To reexport, should be using (not import)
# Macro to exportall
macro exportall(pkg)
    Expr(:export, names(JuMP)...)
end
@exportall JuMP

# The following is largely inspired from JuMP/test/JuMPExtension.jl

mutable struct StructuredModel <: JuMP.AbstractModel
    # Structured data
    parent::Union{StructuredModel, Nothing}
    children::Dict{Int, StructuredModel}
    probability::Dict{Int, Float64}
    num_scen::Int

    # Model data
    nextvaridx::Int                                 # Next variable index is nextvaridx+1
    variables::Dict{Int, JuMP.AbstractVariable}     # Map varidx -> variable
    varnames::Dict{Int, String}                     # Map varidx -> name
    nextconidx::Int                                 # Next constraint index is nextconidx+1
    constraints::Dict{Int, JuMP.AbstractConstraint} # Map conidx -> variable
    connames::Dict{Int, String}                     # Map varidx -> name
    objective_sense::MOI.OptimizationSense
    objective_function::JuMP.AbstractJuMPScalar
    objdict::Dict{Symbol, Any}                      # Same that JuMP.Model's field `objdict`
    function StructuredModel(; parent=nothing, same_children_as=nothing, id=0,
                             num_scenarios::Int=0,
                             prob::Float64=default_probability(parent))
        if same_children_as !== nothing
            if !isa(same_children_as, StructuredModel)
                error("The JuMP model given for the argument `same_children_as' is not valid. Please create it using the `StructuredModel' function.")
            end
            probability = same_children_as.probability
            children = same_children_as.children
        else
            probability = Dict{Int, Float64}()
            children = Dict{Int, JuMP.Model}()
        end

        model = new(parent, children, probability, num_scenarios,                    # Structured
                    0, Dict{Int, JuMP.AbstractVariable}(),   Dict{Int, String}(),    # Model Variables
                    0, Dict{Int, JuMP.AbstractConstraint}(), Dict{Int, String}(),    # Model Constraints
                    MOI.FEASIBILITY_SENSE, zero(JuMP.GenericAffExpr{Float64, StructuredVariableRef}), # Model objective
                    Dict{Symbol, Any}())                                             # Model objects

        if parent === nothing
            id = 0
        else
            @assert id != 0
            parent.children[id] = model
            parent.probability[id] = prob
        end

        return model
    end
end

#### Structured ####

getparent(model::StructuredModel)      = model.parent
getchildren(model::StructuredModel)    = model.children
getprobability(model::StructuredModel) = model.probability
num_scenarios(model::StructuredModel)  = model.num_scen

default_probability(model::StructuredModel) = 1 / num_scenarios(model)
default_probability(::Nothing) = 1.0

#### Model ####

JuMP.object_dictionary(m::StructuredModel) = m.objdict

# Variables
struct StructuredVariableRef <: JuMP.AbstractVariableRef
    model::StructuredModel # `model` owning the variable
    idx::Int       # Index in `model.variables`
end
Base.broadcastable(v::StructuredVariableRef) = Ref(v)
Base.copy(v::StructuredVariableRef) = v
Base.:(==)(v::StructuredVariableRef, w::StructuredVariableRef) = v.model === w.model && v.idx == w.idx
JuMP.owner_model(v::StructuredVariableRef) = v.model
JuMP.isequal_canonical(v::StructuredVariableRef, w::StructuredVariableRef) = v == w
JuMP.variable_type(::StructuredModel) = StructuredVariableRef
function JuMP.add_variable(m::StructuredModel, v::JuMP.AbstractVariable, name::String="")
    m.nextvaridx += 1
    vref = StructuredVariableRef(m, m.nextvaridx)
    m.variables[vref.idx] = v
    JuMP.set_name(vref, name)
    vref
end
function MOI.delete!(m::StructuredModel, vref::StructuredVariableRef)
    delete!(m.variables, vref.idx)
    delete!(m.varnames, vref.idx)
end
MOI.is_valid(m::StructuredModel, vref::StructuredVariableRef) = vref.idx in keys(m.variables)
JuMP.num_variables(m::StructuredModel) = length(m.variables)

# Internal function
variable_info(vref::StructuredVariableRef) = vref.model.variables[vref.idx].info
function update_variable_info(vref::StructuredVariableRef, info::JuMP.VariableInfo)
    vref.model.variables[vref.idx] = JuMP.ScalarVariable(info)
end

JuMP.has_lower_bound(vref::StructuredVariableRef) = variable_info(vref).has_lb
function JuMP.lower_bound(vref::StructuredVariableRef)
    @assert !JuMP.is_fixed(vref)
    variable_info(vref).lower_bound
end
function JuMP.set_lower_bound(vref::StructuredVariableRef, lower)
    info = variable_info(vref)
    update_variable_info(vref,
                         JuMP.VariableInfo(true, lower,
                                           info.has_ub, info.upper_bound,
                                           info.has_fix, info.fixed_value,
                                           info.has_start, info.start,
                                           info.binary, info.integer))
end
function JuMP.delete_lower_bound(vref::StructuredVariableRef)
    info = variable_info(vref)
    update_variable_info(vref,
                         JuMP.VariableInfo(false, info.lower_bound,
                                           info.has_ub, info.upper_bound,
                                           info.has_fix, info.fixed_value,
                                           info.has_start, info.start,
                                           info.binary, info.integer))
end
JuMP.has_upper_bound(vref::StructuredVariableRef) = variable_info(vref).has_ub
function JuMP.upper_bound(vref::StructuredVariableRef)
    @assert !JuMP.is_fixed(vref)
    variable_info(vref).upper_bound
end
function JuMP.set_upper_bound(vref::StructuredVariableRef, upper)
    info = variable_info(vref)
    update_variable_info(vref,
                         JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                           true, upper,
                                           info.has_fix, info.fixed_value,
                                           info.has_start, info.start,
                                           info.binary, info.integer))
end
function JuMP.delete_upper_bound(vref::StructuredVariableRef)
    info = variable_info(vref)
    update_variable_info(vref,
                         JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                           false, info.upper_bound,
                                           info.has_fix, info.fixed_value,
                                           info.has_start, info.start,
                                           info.binary, info.integer))
end
JuMP.is_fixed(vref::StructuredVariableRef) = variable_info(vref).has_fix
JuMP.fix_value(vref::StructuredVariableRef) = variable_info(vref).fixed_value
function JuMP.fix(vref::StructuredVariableRef, value)
    info = variable_info(vref)
    update_variable_info(vref,
                         JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                           info.has_ub, info.upper_bound,
                                           true, value,
                                           info.has_start, info.start,
                                           info.binary, info.integer))
end
function JuMP.unfix(vref::StructuredVariableRef)
    info = variable_info(vref)
    update_variable_info(vref,
                         JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                           info.has_ub, info.upper_bound,
                                           false, info.fixed_value,
                                           info.has_start, info.start,
                                           info.binary, info.integer))
end
JuMP.start_value(vref::StructuredVariableRef) = variable_info(vref).start
function JuMP.set_start_value(vref::StructuredVariableRef, start)
    info = variable_info(vref)
    update_variable_info(vref,
                         JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                           info.has_ub, info.upper_bound,
                                           info.has_fix, info.fixed_value,
                                           true, start,
                                           info.binary, info.integer))
end
JuMP.is_binary(vref::StructuredVariableRef) = variable_info(vref).binary
function JuMP.set_binary(vref::StructuredVariableRef)
    @assert !JuMP.is_integer(vref)
    info = variable_info(vref)
    update_variable_info(vref,
                         JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                           info.has_ub, info.upper_bound,
                                           info.has_fix, info.fixed_value,
                                           info.has_start, info.start,
                                           true, info.integer))
end
function JuMP.unset_binary(vref::StructuredVariableRef)
    info = variable_info(vref)
    update_variable_info(vref,
                         JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                           info.has_ub, info.upper_bound,
                                           info.has_fix, info.fixed_value,
                                           info.has_start, info.start,
                                           false, info.integer))
end
JuMP.is_integer(vref::StructuredVariableRef) = variable_info(vref).integer
function JuMP.set_integer(vref::StructuredVariableRef)
    @assert !JuMP.is_binary(vref)
    info = variable_info(vref)
    update_variable_info(vref,
                         JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                           info.has_ub, info.upper_bound,
                                           info.has_fix, info.fixed_value,
                                           info.has_start, info.start,
                                           info.binary, true))
end
function JuMP.unset_integer(vref::StructuredVariableRef)
    info = variable_info(vref)
    update_variable_info(vref,
                         JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                           info.has_ub, info.upper_bound,
                                           info.has_fix, info.fixed_value,
                                           info.has_start, info.start,
                                           info.binary, false))
end

# Constraints
struct StructuredConstraintRef
    model::StructuredModel # `model` owning the constraint
    idx::Int       # Index in `model.constraints`
end
JuMP.constraint_type(::StructuredModel) = StructuredConstraintRef
function JuMP.add_constraint(m::StructuredModel, c::JuMP.AbstractConstraint, name::String="")
    m.nextconidx += 1
    cref = StructuredConstraintRef(m, m.nextconidx)
    m.constraints[cref.idx] = c
    JuMP.set_name(cref, name)
    cref
end
function MOI.delete!(m::StructuredModel, cref::StructuredConstraintRef)
    delete!(m.constraints, cref.idx)
    delete!(m.connames, cref.idx)
end
MOI.is_valid(m::StructuredModel, cref::StructuredConstraintRef) = cref.idx in keys(m.constraints)
function JuMP.constraint_object(cref::StructuredConstraintRef, F::Type, S::Type)
    c = cref.model.constraints[cref.idx]
    # `TypeError` should be thrown is `F` and `S` are not correct
    # This is needed for the tests in `constraints.jl`
    c.func::F
    c.set::S
    c
end

# Objective
function JuMP.set_objective(m::StructuredModel, sense::MOI.OptimizationSense,
                            f::JuMP.AbstractJuMPScalar)
    m.objective_sense = sense
    m.objective_function = f
end
JuMP.objective_sense(m::StructuredModel) = m.objective_sense
function JuMP.objective_function(m::StructuredModel, FT::Type)
    m.objective_function isa FT || error("The objective function is not of type $FT")
    m.objective_function
end

# Names
JuMP.name(vref::StructuredVariableRef) = vref.model.varnames[vref.idx]
function JuMP.set_name(vref::StructuredVariableRef, name::String)
    vref.model.varnames[vref.idx] = name
end
JuMP.name(cref::StructuredConstraintRef) = cref.model.connames[cref.idx]
function JuMP.set_name(cref::StructuredConstraintRef, name::String)
    cref.model.connames[cref.idx] = name
end

include("BendersBridge.jl")

end


text:
const ValidVarInput = Union{Pair{Symbol,<:ValType}, Vector{<:Pair{Symbol,<:ValType}},
    Dict{Symbol,<:ValType}}
const BlockOrVar = Union{AbstractBlock, Symbol}

struct SequenceSpaceModel <: AbstractGraph{Int}
    dag::SimpleDiGraph{Int}
    pool::Vector{BlockOrVar}
    invpool::Dict{BlockOrVar,Int}
    order::Vector{Int}
    srcs::Set{Int}
    sssrcs::Set{Int}
    dests::Set{Int}
    function SequenceSpaceModel(dag::SimpleDiGraph{Int}, pool::Vector{BlockOrVar},
            invpool::Dict{BlockOrVar,Int})
        order = topological_sort_by_dfs(dag)
        srcs = Set{Int}()
        sssrcs = Set{Int}()
        dests = Set{Int}()
        for v in vertices(dag)
            pool[v] isa AbstractBlock && continue
            ins = inneighbors(dag, v)
            if isempty(ins)
                push!(srcs, v)
                for o in outneighbors(dag, v)
                    if pool[v] in ssinputs(pool[o])
                        push!(sssrcs, v)
                        break
                    end
                end
            elseif length(ins) > 1
                var = pool[v]
                error("$var is the output of $(length(ins)) blocks")
            elseif isempty(outneighbors(dag, v))
                push!(dests, v)
            end
        end
        return new(dag, pool, invpool, order, srcs, sssrcs, dests)
    end
end

function model(blocks)
    blocks isa BlockOrVar && (blocks = (blocks,))
    pool = collect(BlockOrVar, blocks)
    n = length(pool)
    invpool = Dict{BlockOrVar,Int}(pool.=>1:n)
    edges = Edge{Int}[]
    for (ib, b) in enumerate(blocks)
        vins = inputs(b)
        vouts = outputs(b)
        nvin = length(vins)
        for (c, v) in enumerate((vins..., vouts...))
            iv = get(invpool, v, 0)
            if iv === 0
                n += 1
                push!(pool, v)
                invpool[v] = n
                iv = n
            end
            push!(edges, c<=nvin ? Edge(iv,ib) : Edge(ib,iv))
        end
    end
    dag = SimpleDiGraphFromIterator(edges)
    return SequenceSpaceModel(dag, pool, invpool)
end

SimpleDiGraph(m::SequenceSpaceModel) = m.dag

edgetype(m::SequenceSpaceModel) = edgetype(m.dag)
eltype(m::SequenceSpaceModel) = eltype(m.dag)
nv(m::SequenceSpaceModel) = nv(m.dag)
ne(m::SequenceSpaceModel) = ne(m.dag)
vertices(m::SequenceSpaceModel) = vertices(m.dag)
edges(m::SequenceSpaceModel) = edges(m.dag)
is_directed(::Type{SequenceSpaceModel}) = true
has_vertex(m::SequenceSpaceModel, v) = has_vertex(m.dag, v)
has_edge(m::SequenceSpaceModel, s, d) = has_edge(m.dag, s, d)
inneighbors(m::SequenceSpaceModel, v) = inneighbors(m.dag, v)
outneighbors(m::SequenceSpaceModel, v) = outneighbors(m.dag, v)
zero(::Type{SequenceSpaceModel}) =
    SequenceSpaceModel(zero(SimpleDiGraph), BlockOrVar[], Dict{BlockOrVar,Int}())

srcs(m::SequenceSpaceModel) = getfield(m, :srcs)
sssrcs(m::SequenceSpaceModel) = getfield(m, :sssrcs)
dests(m::SequenceSpaceModel) = getfield(m, :dests)
isblock(m::SequenceSpaceModel, v) = m.pool[v] isa AbstractBlock

show(io::IO, m::SequenceSpaceModel) = print(io, "{$(nv(m)), $(ne(m))} ", typeof(m).name.name)

function show(io::IO, ::MIME"text/plain", m::SequenceSpaceModel)
    isb = isa.(m.pool, AbstractBlock)
    nb = sum(isb)
    nvar = length(m.pool) - nb
    print(io, "{$(nv(m)), $(ne(m))} ", typeof(m).name.name)
    print(io, " with $nb block")
    nb > 1 && print(io, 's')
    print(io, " and $nvar variable")
    nvar > 1 && print(io, 's')
    if nb > 0
        print(io, ':')
        for b in view(m.pool, isb)
            print(io, "\n  ", b)
        end
    end
end

struct SteadyState{TF<:AbstractFloat}
    parent::SequenceSpaceModel
    blks::Vector{AbstractBlock}
    vars::Vector{Symbol}
    ins::Vector{Symbol}
    outs::Vector{Symbol}
    calis::Dict{Symbol,ValType{TF}}
    tars::Dict{Symbol,ValType{TF}}
    varvals::Dict{Symbol,ValType{TF}}
    inits::Vector{TF}
    resids::Vector{TF}
    targeted::BitVector
end

function _collapse!(m::SequenceSpaceModel, calis::Dict, varvals::Dict)
    unknowns = Set{Int}()
    ins = Symbol[]
    outs = Symbol[]
    blks = AbstractBlock[]
    visitedsrc = 0
    Nsrc = length(srcs(m))
    for v in m.order
        if isblock(m, v)
            b = m.pool[v]
            gins = inneighbors(m.dag, v)
            for gi in gins
                if gi in unknowns
                    push!(unknowns, v)
                    push!(blks, b)
                    break
                end
            end
            steadystate!(varvals, b)
        else
            var = m.pool[v]
            cali = get(calis, var, nothing)
            if visitedsrc < Nsrc && v in srcs(m)
                visitedsrc += 1
                if v in sssrcs(m) && cali === nothing
                    msg = "initial value for input variable $var is not assigned"
                    haskey(varvals, var) || @warn msg
                    push!(unknowns, v)
                    push!(ins, var)
                end
            else
                gin = inneighbors(m.dag, v)[1]
                if gin in unknowns
                    push!(unknowns, v)
                    v in dests(m) && push!(outs, var)
                    if cali !== nothing
                        @warn "calibrated value for $var is replaced by output from its parent block"
                        delete!(calis, var)
                    end
                end
            end
            # Only assign calibrated values after the checks
            cali === nothing || (varvals[var] = cali)
        end
    end
    return ins, outs, blks
end

function SteadyState(m::SequenceSpaceModel, calibrated::ValidVarInput,
        targets::ValidVarInput, initials::Union{ValidVarInput,Nothing}=nothing,
        TF::Type=Float64)
    calibrated isa Pair && (calibrated = (calibrated,))
    targets isa Pair && (targets = (targets,))
    initials !== nothing && initials isa Pair && (initials = (initials,))
    calis = Dict{Symbol,ValType{TF}}(calibrated...)
    vars = Symbol[v for v in m.pool if v isa Symbol]
    varvals = Dict{Symbol,ValType{TF}}()
    # Assign initials to varvals helps look up assigned values when filling inits
    # This is also the way to tell whether any unknown variable is a Vector
    if initials !== nothing
        for (k, v) in initials
            varvals[k] = v
        end
    end

    ins, outs, blks = _collapse!(m, calis, varvals)

    tars = Dict{Symbol,ValType{TF}}()
    for (k, v) in targets
        if k in outs
            tars[k] = v
        else
            @warn "target value for $k is ignored because $k is an input of a block; consider adding a block for the gap between the target"
        end
    end
    # Assign a scalar zero whenever the initial value is not provided
    z = zero(TF)
    inits = Vector{TF}(undef, length(ins))
    i0 = 0
    if initials !== nothing
        @inbounds for vi in ins
            val = get(varvals, vi, z)
            w = length(val)
            if w > 1
                resize!(inits, length(inits)+w-1)
                inits[i0+1:i0+w] .= val
            else
                inits[i0+1] = val
            end
            i0 += w
        end
    end
    # Some variables in outs might not be targeted
    resids = Vector{TF}(undef, length(tars))
    targeted = haskey.(Ref(tars), outs)
    return SteadyState(m, blks, vars, ins, outs, calis, tars, varvals, inits, resids, targeted)
end

function _inputs!(ss::SteadyState, inputs::AbstractVector)
    N = length(inputs)
    i0 = 0
    @inbounds for vi in ss.ins
        val = get(ss.varvals, vi, nothing)
        w = length(val)
        i0+w > N && error("length of inputs does not match varvals")
        if val isa Vector
            for i in 1:w
                val[i] = inputs[i0+i]
            end
        else
            ss.varvals[vi] = inputs[i0+1]
        end
        i0 += w
    end
end

function _resids!(resids::AbstractVector, ss::SteadyState)
    i0 = 0
    @inbounds for out in ss.outs
        tar = get(ss.tars, out, nothing)
        tar === nothing && continue
        val = ss.varvals[out]
        w = length(val)
        ilast = i0 + w
        ilast > length(resids) && resize!(resids, ilast)
        if val isa Vector
            # tar might be a scalar
            resids[i0+1:i0+w] .= val .- tar
        else
            resids[i0+1] = val - tar
        end
        i0 += w
    end
end

function residuals!(resids::AbstractVector, ss::SteadyState)
    for b in ss.blks
        steadystate!(ss.varvals, b)
    end
    _resids!(resids, ss)
    return resids
end

residuals!(ss::SteadyState) = residuals!(ss.resids, ss)

function residuals!(resids::AbstractVector, ss::SteadyState, inputs::AbstractVector)
    _inputs!(ss, inputs)
    return residuals!(resids, ss)
end

residuals!(ss::SteadyState, inputs::AbstractVector) =
    residuals!(ss.resids, ss, inputs)

function criterion!(ss::SteadyState; weight::Union{AbstractMatrix,UniformScaling}=I)
    resids = residuals!(ss)
    return resids'*weight*resids
end

function criterion!(ss::SteadyState, inputs::AbstractVector;
        weight::Union{AbstractMatrix,UniformScaling}=I)
    resids = residuals!(ss, inputs)
    return resids'*weight*resids
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -24.652335523355834 <= q <= 21.32472335033785)

                     add_NL_constraint(m, :(gelu(-0.6982139593401175 + -0.19504968807074974*gelu(0.6793369023557885 + 0.6523875537646115*gelu(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])) + -0.5142038874540145*gelu(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])) + -0.0456439541419269*gelu(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])) + 0.6096854328317232*gelu(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))) + -0.41803298572193937*gelu(-0.7048442383777926 + 0.47781322839334983*gelu(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])) + -0.713369260368808*gelu(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])) + -0.44964139924772173*gelu(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])) + -0.9876898123617992*gelu(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))) + 0.9542284451941447*gelu(0.9969387928417852 + -0.011395717836928387*gelu(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])) + 0.8103499615759371*gelu(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])) + -0.12667586987008894*gelu(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])) + 0.6252342222679741*gelu(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))) + 0.7916680056214345*gelu(-0.5156629799588863 + -0.22745979325095078*gelu(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])) + -0.8222624364306839*gelu(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])) + -0.4937419854366709*gelu(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])) + -0.1046356011041838*gelu(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3])))) + gelu(-0.4539405071565561 + 0.5320780591355936*gelu(0.6793369023557885 + 0.6523875537646115*gelu(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])) + -0.5142038874540145*gelu(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])) + -0.0456439541419269*gelu(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])) + 0.6096854328317232*gelu(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))) + -0.8489788623944037*gelu(-0.7048442383777926 + 0.47781322839334983*gelu(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])) + -0.713369260368808*gelu(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])) + -0.44964139924772173*gelu(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])) + -0.9876898123617992*gelu(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))) + -0.25620536200506905*gelu(0.9969387928417852 + -0.011395717836928387*gelu(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])) + 0.8103499615759371*gelu(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])) + -0.12667586987008894*gelu(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])) + 0.6252342222679741*gelu(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))) + 0.9264345634068367*gelu(-0.5156629799588863 + -0.22745979325095078*gelu(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])) + -0.8222624364306839*gelu(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])) + -0.4937419854366709*gelu(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])) + -0.1046356011041838*gelu(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3])))) + gelu(-0.4017849565754932 + 0.7092320156791541*gelu(0.6793369023557885 + 0.6523875537646115*gelu(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])) + -0.5142038874540145*gelu(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])) + -0.0456439541419269*gelu(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])) + 0.6096854328317232*gelu(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))) + -0.20657569919094643*gelu(-0.7048442383777926 + 0.47781322839334983*gelu(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])) + -0.713369260368808*gelu(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])) + -0.44964139924772173*gelu(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])) + -0.9876898123617992*gelu(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))) + 0.6260976974386958*gelu(0.9969387928417852 + -0.011395717836928387*gelu(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])) + 0.8103499615759371*gelu(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])) + -0.12667586987008894*gelu(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])) + 0.6252342222679741*gelu(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))) + 0.026773003462286393*gelu(-0.5156629799588863 + -0.22745979325095078*gelu(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])) + -0.8222624364306839*gelu(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])) + -0.4937419854366709*gelu(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])) + -0.1046356011041838*gelu(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3])))) + gelu(0.8541923459688499 + -0.7979821016654509*gelu(0.6793369023557885 + 0.6523875537646115*gelu(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])) + -0.5142038874540145*gelu(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])) + -0.0456439541419269*gelu(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])) + 0.6096854328317232*gelu(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))) + 0.2172669330681649*gelu(-0.7048442383777926 + 0.47781322839334983*gelu(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])) + -0.713369260368808*gelu(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])) + -0.44964139924772173*gelu(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])) + -0.9876898123617992*gelu(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))) + -0.6626207135044284*gelu(0.9969387928417852 + -0.011395717836928387*gelu(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])) + 0.8103499615759371*gelu(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])) + -0.12667586987008894*gelu(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])) + 0.6252342222679741*gelu(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))) + -0.5735396155482761*gelu(-0.5156629799588863 + -0.22745979325095078*gelu(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])) + -0.8222624364306839*gelu(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])) + -0.4937419854366709*gelu(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])) + -0.1046356011041838*gelu(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


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

using Markdown
using InteractiveUtils

# ╔═╡ 6387ea68-f280-40b7-be3f-7b8241a39624
using CairoMakie, Random, DataFrames, FileIO

# ╔═╡ 7a5f7917-e371-40a0-b684-f3feac6bec7c
md"## CairoMakie demo"

# ╔═╡ 9b4d3474-614b-4889-baa8-38471d4c76e0
let
    Random.seed!(123)
    fig = Figure(resolution = (400, 427))
    ax = Axis(fig,  aspect = 1, xlabel = "x", ylabel = "y")
    hmap = heatmap!(rand(20,20), colormap = :lajolla)
    cbar = Colorbar(fig, hmap, label = "values", height = 15, vertical = false,
     flipaxis = false, ticksize=15, tickalign = 1, width = Relative(3.55/4))
    fig[1, 1] = ax
    fig[2, 1] = cbar
    rowgap!(fig.layout, 7)
    fig
end

# ╔═╡ 5e7a5b8f-73fe-4ec0-8d28-79195ebaf8a9
md"My data to use on the same demo layout"

# ╔═╡ a86eee2c-7d8d-42fd-abaf-b9cecb56d065
begin
	L = [480, 1090]
# The logo dewetting data
	logo_data = load("E:\\JuliaStuff\\Data_PhD\\Tutorials\\Logo_dewetting.jld2") |> DataFrame
end

# ╔═╡ b764b630-889b-47de-a151-532203804cc8
"""
	corr_disp(arr)

Rotates and inverts an snapshot of the simulation to be correctly displayed in Makie.
"""
function corr_disp(arr; time=15000)
	new_arr = zeros(L[1], L[2])
	new_arr = reverse(rotr90(reshape(arr[!, Symbol("h_$(time)")],L[1],L[2])), dims=2)

	return new_arr
end

# ╔═╡ 8a108953-eb78-43fd-aac0-0afa0ba563a9
let
    fig = Figure(resolution = (860, 680))
    ax = Axis(fig,  aspect = 2, xlabel = "x", ylabel = "y")
    hmap = heatmap!(corr_disp(logo_data, time=20000), colormap = :viridis)
    cbar = Colorbar(fig, hmap, label = "thickness", height = 15, vertical = false,
     flipaxis = false, ticksize=15, tickalign = 1, width = Relative(3.55/4))
    fig[1, 1] = ax
    fig[2, 1] = cbar
    rowgap!(fig.layout, 7)
    fig
end

# ╔═╡ c0662743-134d-4904-8983-f3ee1bd26c59
let
    Random.seed!(123)
    n = 15
    x, y, color = rand(n), rand(n), rand(n)
    cmaps = [:cool, :viridis, :inferno, :thermal]

    function FigGridHeatSharedCbarH()
        fig = Figure(resolution = (600, 400))
        axes = []
        c = 1
		m = 3
        for i in 1:1, j in 1:m
            ax = Axis(fig[i, j],aspect=1,xgridvisible = false,ygridvisible = false)
            pnts = heatmap!(rand(10,10), colormap= cmaps[2], colorrange=(0, 1))
            #ax.xticks = [1,10]
            #ax.yticks = [1,10]
            cbar = Colorbar(fig, pnts, vertical = false, flipaxis = false,
            width = Relative(3/4), height = 15, tickwidth = 2, ticklabelsize = 14)
            cbar.ticks = [0,0.5,1]
            fig[2, 1:m] = cbar
            c+=1
            push!(axes, ax)
        end
		for k in 2:3
        	hideydecorations!(axes[k], ticks = true)
		end
        fig
    end
    fig = FigGridHeatSharedCbarH()
    fig
end

# ╔═╡ a2d0c8c6-2063-4ed1-ad91-3e7ff08c4743
let
	noto_sans = assetpath("fonts", "NotoSans-Regular.ttf")
	noto_sans_bold = assetpath("fonts", "NotoSans-Bold.ttf")
    function FigGridHeatSharedCbarHLogo()

		# Defining the figure window
        fig = Figure(resolution = (1300, 450), font=noto_sans)
        axes = []
        c = 1
		m = 3
		# Data for the three heatmap images
		imgs = [corr_disp(logo_data, time=5000), corr_disp(logo_data, time=20000), corr_disp(logo_data, time=35000)]
		# Labels of the heatmaps
		times = ["5000 Δt", "20000 Δt", "35000 Δt"]
		# Gauge of the colormap
		h_lim = maximum(imgs[3])
		# Loop over the three images
        for i in 1:1, j in 1:m
            ax = Axis(fig[i, j],
				      aspect = 2,
				      xgridvisible = false,
				      ygridvisible = false,
				      xlabel = "x Δx",
				      ylabel = "y Δx",
					  xticklabelsize=25,
					  yticklabelsize=25,
					  )
            pnts = heatmap!(imgs[j], colormap=:viridis, colorrange=(0, h_lim))
            ax.xticks = [1,512,1090]
            ax.yticks = [1,256,480]
			# ax.fontsize = 26
            # ax.yticks = [1,256,480]
			ax.title = times[j]
			ax.titlesize = 34
			ax.xlabelsize = 28
			ax.ylabelsize = 28
            cbar = Colorbar(fig, pnts, label="Film Thickness", vertical = false, flipaxis = false,
            width = Relative(4/5), height = 15, tickwidth = 2, ticklabelsize = 24, labelsize=26)
            cbar.ticks = [0,0.75,1.5,2.25]
            fig[2, 1:m] = cbar
            c+=1
            push!(axes, ax)
        end
		for k in 2:3
        	hideydecorations!(axes[k], ticks = true)
        	hidexdecorations!(axes[k], ticks = true)

		end
        fig
    end
    fig = FigGridHeatSharedCbarHLogo()
    fig
	# save("E:\\JuliaStuff\\Figures\\Hiern_logo.png", fig)
end

# ╔═╡ a198a1eb-c164-4ab9-b0ee-8296f080846e
md"## Let's try an animation

Would be cool to have an animation of changing heatmaps!"

# ╔═╡ 39c23799-6f5d-402a-b8dd-9404f78357f4
begin
	time = Node(0.0)

	xs = range(0, 7, length=40)

	ys_1 = @lift(sin.(xs .- $time))
	ys_2 = @lift(cos.(xs .- $time) .+ 3)

	fig = lines(xs, ys_1, color = :blue, linewidth = 4,
    	        axis = (title = @lift("t = $(round($time, digits = 1))"),))
	scatter!(xs, ys_2, color = :red, markersize = 15)

	framerate = 30
	timestamps = range(0, 2, step=1/framerate)

	record(fig, "time_animation.mp4", timestamps; framerate = framerate) do t
    	time[] = t
	end
end

# ╔═╡ cd57c594-a9ca-48cd-ac05-e4a547d93d05
begin
	time_n = Node(1000)

	time_st = collect(1000:500:50000)
	max_h = maximum(logo_data[!, :h_50000])
	fig_logo = Figure(resolution = (860, 680))
    ax = Axis(fig_logo,
		      aspect = 2,
		      xlabel = "x Δx",
		      ylabel = "y Δx",
		      title = @lift("t = $($time_n) Δt"),
			  xticklabelsize=25,
			  yticklabelsize=25,
			  xlabelsize = 28,
			  ylabelsize = 28,
		      titlesize = 34)
	data_logo = @lift(corr_disp(logo_data, time=Int($time_n)))
    hmap = heatmap!(data_logo, colormap = :viridis, colorrange=(0, max_h))
    cbar = Colorbar(fig_logo, hmap, label = "Film Thickness", height = 15, vertical = false,
     flipaxis = false, ticksize=15, tickalign = 1, width = Relative(3.55/4), ticklabelsize = 24, labelsize=26)
    fig_logo[1, 1] = ax
    fig_logo[2, 1] = cbar
    rowgap!(fig_logo.layout, 7)

    record(fig_logo, "logo_animation.gif", time_st; framerate = 30) do t
    	time_n[] = t
	end
end

# ╔═╡ Cell order:
# ╠═6387ea68-f280-40b7-be3f-7b8241a39624
# ╟─7a5f7917-e371-40a0-b684-f3feac6bec7c
# ╠═9b4d3474-614b-4889-baa8-38471d4c76e0
# ╟─5e7a5b8f-73fe-4ec0-8d28-79195ebaf8a9
# ╠═a86eee2c-7d8d-42fd-abaf-b9cecb56d065
# ╠═b764b630-889b-47de-a151-532203804cc8
# ╠═8a108953-eb78-43fd-aac0-0afa0ba563a9
# ╠═c0662743-134d-4904-8983-f3ee1bd26c59
# ╠═a2d0c8c6-2063-4ed1-ad91-3e7ff08c4743
# ╠═a198a1eb-c164-4ab9-b0ee-8296f080846e
# ╠═39c23799-6f5d-402a-b8dd-9404f78357f4
# ╠═cd57c594-a9ca-48cd-ac05-e4a547d93d05

max_stars_repo_name: Zitzeronion/NotebooksnStuff

max_issues_repo_name: Zitzeronion/NotebooksnStuff

max_forks_repo_name: Zitzeronion/NotebooksnStuff


text:
import PowerModelsACDC; const _PMACDC = PowerModelsACDC
import PowerModels; const _PM = PowerModels
import InfrastructureModels; const _IM = InfrastructureModels
import Ipopt
# using CPLEX
import SCS
import Juniper
import Mosek
import MosekTools
import JuMP
import Gurobi
import Cbc
import CPLEX

file = "./test/data/tnep/case4_original.m"
file_acdc = "./test/data/tnep/case4_acdc.m"
data = _PM.parse_file(file)
_PMACDC.process_additional_data!(data)

data_bf=data
scs = JuMP.optimizer_with_attributes(SCS.Optimizer, "max_iters" => 100000)
ipopt = JuMP.optimizer_with_attributes(Ipopt.Optimizer, "tol" => 1e-4, "print_level" => 0)
cplex = JuMP.optimizer_with_attributes(CPLEX.Optimizer)
cbc = JuMP.optimizer_with_attributes(Cbc.Optimizer, "tol" => 1e-4, "print_level" => 0)
gurobi = JuMP.optimizer_with_attributes(Gurobi.Optimizer)
mosek = JuMP.optimizer_with_attributes(Mosek.Optimizer)
juniper = JuMP.optimizer_with_attributes(Juniper.Optimizer, "nl_solver" => ipopt, "mip_solver" => cbc, "time_limit" => 7200)
#
#
s = Dict("output" => Dict("branch_flows" => true), "conv_losses_mp" => false, "process_data_internally" => false)

resultDC = _PMACDC.run_tnepopf(file, _PM.DCPPowerModel, gurobi, setting = s)
resultAC = _PMACDC.run_tnepopf(file, _PM.ACPPowerModel, juniper, setting = s)
resultSOCBF = _PMACDC.run_tnepopf_bf(file, _PM.SOCBFPowerModel, gurobi, setting = s)
resultSOCWR = _PMACDC.run_tnepopf(file, _PM.SOCWRPowerModel, gurobi, setting = s)
resultQC     =  _PMACDC.run_tnepopf(file, _PM.QCRMPowerModel, gurobi; setting = s)
# resultSDP     =  _PMACDC.run_tnepopf(file, _PM.SDPWRMPowerModel, mosek; setting = s)
resultLPAC     =  _PMACDC.run_tnepopf(file, _PM.LPACCPowerModel, gurobi; setting = s)
#
_PMACDC.display_results_tnep(resultDC)
_PMACDC.display_results_tnep(resultAC)
_PMACDC.display_results_tnep(resultSOCBF)
_PMACDC.display_results_tnep(resultSOCWR)
_PMACDC.display_results_tnep(resultLPAC)
_PMACDC.display_results_tnep(resultQC)
## TEST ACDC TNEP
s = Dict("output" => Dict("branch_flows" => true), "conv_losses_mp" => true,"process_data_internally" => false)
resultACDC_dcp = _PMACDC.run_acdctnepopf(file_acdc, _PM.DCPPowerModel, gurobi, setting = s)
resultACDC_acp = _PMACDC.run_acdctnepopf(file_acdc, _PM.ACPPowerModel, juniper, setting = s)
#resultACDC_socbf = _PMACDC.run_acdctnepopf_bf(file_acdc, _PM.SOCBFPowerModel, gurobi, setting = s)  # BF TNEP not implemented in PowerModels.jl
resultACDC_socwr = _PMACDC.run_acdctnepopf(file_acdc, _PM.SOCWRPowerModel, gurobi, setting = s)
resultACDC_qc = _PMACDC.run_acdctnepopf(file_acdc, _PM.QCRMPowerModel, gurobi, setting = s)
resultACDC_lpac = _PMACDC.run_acdctnepopf(file_acdc, _PM.LPACCPowerModel, gurobi, setting = s)

t = 1:2
function build_mn_data(file)
    mp_data = _PM.parse_file(file)
    return _IM.replicate(mp_data, length(t), Set{String}(["source_type", "name", "source_version", "per_unit"]))
end

data1 = build_mn_data(file)
_PMACDC.process_additional_data!(data1)
data_acdc = build_mn_data(file_acdc)
_PMACDC.process_additional_data!(data_acdc)

resultDC1 = _PMACDC.run_mp_tnepopf(data1, _PM.DCPPowerModel, gurobi, multinetwork=true; setting = s)
resultAC1 = _PMACDC.run_mp_tnepopf(data1, _PM.ACPPowerModel, juniper, multinetwork=true; setting = s)
resultBF1 = _PMACDC.run_mp_tnepopf_bf(data1, _PM.SOCBFConicPowerModel, mosek, multinetwork=true; setting = s)
resultLPAC1 = _PMACDC.run_mp_tnepopf(data1, _PM.LPACCPowerModel, juniper, multinetwork=true; setting = s)
resultSOCWR1 = _PMACDC.run_mp_tnepopf(data1, _PM.SOCWRPowerModel, gurobi, multinetwork=true; setting = s)
resultSOCBF1 = _PMACDC.run_mp_tnepopf_bf(data1, _PM.SOCBFConicPowerModel, mosek, multinetwork=true; setting = s)

mp_resultACDC_acp = _PMACDC.run_mp_acdctnepopf(data_acdc, _PM.ACPPowerModel, juniper, multinetwork=true; setting = s)
mp_resultACDC_dcp = _PMACDC.run_mp_acdctnepopf(data_acdc, _PM.DCPPowerModel, gurobi, multinetwork=true; setting = s)
mp_resultACDC_socwr = _PMACDC.run_mp_acdctnepopf(data_acdc, _PM.SOCWRPowerModel, gurobi, multinetwork=true; setting = s)
mp_resultACDC_qc = _PMACDC.run_mp_acdctnepopf(data_acdc, _PM.QCRMPowerModel, gurobi, multinetwork=true; setting = s)
# mp_resultACDC_socbf = _PMACDC.run_mp_acdctnepopf_bf(data_acdc, _PM.SOCBFPowerModel, gurobi, multinetwork=true; setting = s) # BF TNEP not implemented in PowerModels.jl
mp_resultACDC_lpac = _PMACDC.run_mp_acdctnepopf(data_acdc, _PM.LPACCPowerModel, gurobi, multinetwork=true; setting = s)

data1["nw"]["2"]["load"]["2"]["pd"] = 0.5
resultDC2 = _PMACDC.run_mp_tnepopf(data1, _PM.DCPPowerModel, gurobi, multinetwork=true; setting = s)
resultAC2 = _PMACDC.run_mp_tnepopf(data1, _PM.ACPPowerModel, juniper, multinetwork=true; setting = s)
resultBF2 = _PMACDC.run_mp_tnepopf_bf(data1, _PM.SOCBFConicPowerModel, mosek, multinetwork=true; setting = s)


text:
#= Tools for working with levels by name =#

"""
    levelname_map(indexer)::Vector{Pair{Symbol, Int}}

This should be overloaded by any indexer for a MultiResolutionIterator
To map from level names to numbers,

An indexer could be a type that represents a Corpus for example.
It could be the top level of the MultiResolutionIterator itself,
if that has a unique type.
But it doesn't have to be. It just needs to be a handle we can dispatch on

The returned vector maps Named levels to their numbers
Multiple names can be given to the name nmbered level (Many to One)
"""
function levelname_map end

# include the levels
lvls(indexer) = Int[]
lvls(indexer, nums::Vararg{<:Integer}) = nums
function lvls(indexer, names::Vararg{Symbol})
    name_set = first.(levelname_map(indexer))
    num_set = last.(levelname_map(indexer))
    map(names) do name
        ind = findfirst(isequal(name), name_set)
        num_set[ind]
    end
end

function lvls(indexer, named_level_ops::Dict)
    nums = lvls(indexer, keys(named_level_ops)...)
    Dict(zip(nums, values(named_level_ops)))
end

################
# exclude3 the levels using `(!lvls)(args...)`

not_lvls(indexer) = unique(last.(levelname_map(indexer)))

function not_lvls(indexer, nums::Vararg{<:Integer})
    all_levels = (!lvls)(indexer)
    setdiff(all_levels, nums)
    # all the levels, except the ones we are excluding
end

function not_lvls(indexer, names::Vararg{Symbol})
    (!lvls)(indexer, lvls(indexer, names...)...)
end

# make `(!lvls)`  the same as `not_lvls`
(Base.:!(::typeof(lvls))) = not_lvls

#


text:
"""
    backbone(name::AbstractString)

Returns the backbone entry for a taxon, matched by *exact* name.
"""
function backbone(name::AbstractString)
    return first(backbones(Pair("name", name)))
end


text:
## Activate environment and import packages
cd(@__DIR__)
using Pkg
Pkg.activate(".")
using DirectDetections, Distributions

###################################
## More sophisticated model: fit orbit and mass using both astrometry and GAIA-HIPPARCOS proper motion anomaly.
@named b = Planet(
    Priors(
        a = Uniform(1, 75),
        i = Sine(),
        e = Uniform(0,  1),
        Ω = Uniform(0,  π),
        ω = Uniform(0, 2π),
        τ = Uniform(0, 1),
        # Mass in Jupiter-masses. Log-uniform is a reasonable uninformative prior.
        mass = Uniform(0.1, 100)
    ),
    Astrometry(
        (epoch=57009.0, ra= 69.3356, dec=-448.917, σ_ra= 1.82477,  σ_dec= 1.8787 ),
        (epoch=57052.0, ra= 78.3783, dec=-444.96 , σ_ra= 2.05054,  σ_dec= 2.05971),
        (epoch=57053.0, ra= 77.8303, dec=-450.121, σ_ra= 2.39716,  σ_dec= 2.565  ),
        (epoch=57054.0, ra= 76.9638, dec=-455.037, σ_ra=24.1568 ,  σ_dec=23.9074 ),
        (epoch=57266.0, ra=100.052 , dec=-443.966, σ_ra= 2.07345,  σ_dec= 2.22231),
        (epoch=57332.0, ra=108.641 , dec=-439.656, σ_ra= 4.56234,  σ_dec= 5.36208),
        (epoch=57374.0, ra=112.918 , dec=-441.705, σ_ra= 4.65729,  σ_dec= 6.13954),
        (epoch=57376.0, ra=112.464 , dec=-440.892, σ_ra= 3.39255,  σ_dec= 3.0549 ),
        (epoch=57415.0, ra=110.406 , dec=-440.845, σ_ra= 4.18857,  σ_dec= 5.93133),
        (epoch=57649.0, ra=142.053 , dec=-432.057, σ_ra= 2.05962,  σ_dec= 2.02432),
        (epoch=57652.0, ra=141.521 , dec=-428.673, σ_ra= 2.46576,  σ_dec= 2.6485 ),
        (epoch=57739.0, ra=153.258 , dec=-422.449, σ_ra= 2.12148,  σ_dec= 2.14627),
        (epoch=58068.0, ra=187.509 , dec=-406.365, σ_ra= 3.04171,  σ_dec= 3.02463),
        (epoch=58442.0, ra=219.468 , dec=-374.674, σ_ra= 1.815  ,  σ_dec= 1.9453 ),
    )
)
@named cEri_pma = System(
    Priors(
        plx = gaia_plx(gaia_id=3205095125321700480),
        M  = TruncatedNormal(1.75, 0.005, 0, Inf),
    ),
    # We compute the proper motion anomaly by comparing the GAIA and HIPPARCOS proper motions against their
    # long-term trend in position. This is like radial velocity detection, but in the plane of the sky instead
    # of towards/away.
    ProperMotionAnomHGCA(gaia_id=3205095125321700480),
    # And here we list our planets:
    b
)

chains_pma = DirectDetections.hmc(
    cEri_pma, 0.95,
    adaptation=  1_000,
    iterations= 25_000,
    tree_depth=    13,
)

## Plot fitted model (will take ~1 minute the first time), coloured by ecentricity.
using Plots
plotmodel(chains_pma, color=:e, pma_scatter=:mass, xlims=(-750,750), ylims=(-500,750))

##
# using DirectOrbits
# function mplot(chain, planet_key, color, prop)
#     planet = chain.info.model.planets[planet_key]
#     if prop == :ra
#         ylabel = "RA"
#     elseif prop == :dec
#         ylabel = "DEC"
#     elseif prop == :rv
#         ylabel = "RV"
#     elseif prop == :pmra
#         ylabel = "RA/yr"
#     elseif prop == :pmdec
#         ylabel = "DEC/yr"
#     end
#     p1 = plot(;
#         ylabel,
#         legend=:none
#     )
#     N = 500
#     ii = rand(1:size(chain,1)*size(chain,3), N)
#     elements = DirectDetections.construct_elements(chain, planet_key, ii)
#     y = nothing
#     if prop == :ra
#         t = range((extrema(b.astrometry.epoch) .+ [-300, 300])..., length=100)
#         y = planet.astrometry.ra
#         yerr = planet.astrometry.σ_ra
#         fit = raoff.(elements, t')'
#         x = planet.astrometry.epoch
#     elseif prop == :dec
#         t = range((extrema(b.astrometry.epoch) .+ [-300, 300])..., length=100)
#         y = planet.astrometry.dec
#         yerr = planet.astrometry.σ_dec
#         fit = decoff.(elements, t')'
#         x = planet.astrometry.epoch
#     elseif prop == :pmra
#         t = range((extrema(chain.info.model.propermotionanom.ra_epoch) .+ [-300, 300])..., length=100)
#         y = chain.info.model.propermotionanom.pm_ra
#         yerr = chain.info.model.propermotionanom.σ_pm_ra
#         fit = getindex.(
#             propmotionanom.(elements, t', collect(chain["$planet_key[mass]"][ii]).*DirectDetections.mjup2msol),
#             1
#         )'
#         x = chain.info.model.propermotionanom.ra_epoch
#     elseif prop == :pmdec
#         t = range((extrema(chain.info.model.propermotionanom.dec_epoch) .+ [-300, 300])..., length=100)
#         y = chain.info.model.propermotionanom.pm_dec
#         yerr = chain.info.model.propermotionanom.σ_pm_dec
#         fit = getindex.(
#             propmotionanom.(elements, t', collect(chain["$planet_key[mass]"][ii]).*DirectDetections.mjup2msol),
#             2
#         )'
#         x = chain.info.model.propermotionanom.dec_epoch
#     elseif prop == :rv
#         # TODO
#         t = range((extrema(b.astrometry.epoch) .+ [-300, 300])..., length=100)
#         x = planet.astrometry.epoch
#         fit = radvel.(elements, t', collect(chain["$planet_key[mass]"][ii]))'
#     end
#     plot!(
#         t, fit,
#         line_z=repeat(
#             collect(chain["$planet_key[$color]"][ii]),
#             1, length(t)
#         )',
#         alpha=0.05
#     )
#     if !isnothing(y)
#         scatter!(
#             p1,
#             x,
#             y; yerr
#         )
#     end
#     p1
# end
using Plots
plot(
    DirectDetections.timeplot(chains_pma, :b, :mass, :ra),
    DirectDetections.timeplot(chains_pma, :b, :mass, :dec),
    DirectDetections.timeplot(chains_pma, :b, :mass, :rv),
    DirectDetections.timeplot(chains_pma, :b, :mass, :pmra),
    DirectDetections.timeplot(chains_pma, :b, :mass, :pmdec),
    layout = @layout([
        A
        B
        C
        D E
    ]),
    framestyle=:box,
    size=(500,700)
)

## Look at just the mass
histogram(chains_pma["b[mass]"], xlabel="mass (Mⱼᵤₚ)")

## Mass vs. semi-major axis
histogram2d(
    chains_pma["b[a]"   ],
    chains_pma["b[mass]"],
    xlabel="a (AU)",
    ylabel="mass (Mⱼᵤₚ)",
)

##
density(chains_pma["b[mass]"], xlabel="mass (Mⱼᵤₚ)")

##
using PairPlots
corner(chains_pma)

## investigate acceleration
elems = DirectDetections.construct_elements(chains_pma, :b, 1:size(chains_pma,1));

using ForwardDiff
using StaticArrays
using DirectOrbits
function accel(elems, t, mass)
    ∂²ra∂t² =  ForwardDiff.derivative(t->propmotionanom(elems, t, mass)[1], t)
    ∂²dec∂t² =  ForwardDiff.derivative(t->propmotionanom(elems, t, mass)[2], t)
    return SVector(∂²ra∂t²*DirectOrbits.year2days, ∂²dec∂t²*DirectOrbits.year2days)
end
@time accel(elems[1], 1000.0, chains_pma["b[mass]"][1].*DirectDetections.mjup2msol)

##

N = 250
ii = rand(1:size(chains_pma,1), N)

# ts = range(first(b.astrometry.epoch) - 10*365, last(b.astrometry.epoch) + 10*365, length=200)
ts = range(first(cEri_pma.propermotionanom.ra_epoch) - 365, last(cEri_pma.propermotionanom.ra_epoch) + 365, length=200)

##
line_z=repeat(
    collect(chains_pma["b[mass]"][ii]),
    1, length(ts)
)'

vel = propmotionanom.(elems[ii], ts', chains_pma["b[mass]"][ii].*DirectDetections.mjup2msol)';
dra = getindex.(vel, 1)
ddec = getindex.(vel, 2)
pvelra = plot(ts, dra, color=:black, alpha=0.5,title="vel ra", line_z=line_z,label="", legend=false)
scatter!(cEri_pma.propermotionanom.ra_epoch, cEri_pma.propermotionanom.pm_ra, yerr=cEri_pma.propermotionanom.σ_pm_ra)
pveldec = plot(ts, ddec, color=:black, alpha=0.5,title="vel dec", line_z=line_z,label="", colorbar_title="mass (mjup)")
scatter!(cEri_pma.propermotionanom.dec_epoch, cEri_pma.propermotionanom.pm_dec, yerr=cEri_pma.propermotionanom.σ_pm_dec, label="")

acc = accel.(elems[ii], ts', chains_pma["b[mass]"][ii].*DirectDetections.mjup2msol)';
ddra = getindex.(acc, 1)
dddec = getindex.(acc, 2)
paccra = plot(ts, ddra, color=:black, alpha=0.5, title="acc ra", line_z=line_z,label="", legend=false)
paccdec = plot(ts, dddec, color=:black, alpha=0.5, title="acc dec", line_z=line_z,label="", colorbar_title="mass (mjup)")
plts = pvelra, pveldec, paccra, paccdec
vline!.(plts, Ref([extrema(cEri_pma.propermotionanom.ra_epoch)...]), label="")
p=plot(plts..., layout=(2,2))

##
elem_earth = KeplerianElementsDeg(
    a = 1.0,
    μ = 1.0,
    i = 0,
    e = 0,
    ω = 0,
    Ω = 0,
    plx = 1000,
    τ = 0,
)
vel = propmotionanom.(elem_earth, ts, 1 .* DirectDetections.mjup2msol);
dra = getindex.(vel, 1)
acc = accel.(elem_earth, ts,  1 .* DirectDetections.mjup2msol);
ddra = getindex.(acc, 1)
plot(ts, dra, legend=false)
plot!(ts, ddra, legend=false)

max_stars_repo_name: sefffal/teps-problems

max_issues_repo_name: sefffal/teps-problems

max_forks_repo_name: sefffal/teps-problems


text:
@testset "Block" begin
    @testset "Header" begin
        header_raw = hex2bytes("020000208ec39428b17323fa0ddec8e887b4a7c53b8c0a0a220cfd0000000000000000005b0750fce0a889502d40508d39576821155e9c9e3f5c3157f961db38fd8b25be1e77a759e93c0118a4ffd71d")
        io = IOBuffer(header_raw)
        header = Header(io)
        @testset "Parse" begin
            @test header.version == 0x20000002
            want = hex2bytes("8ec39428b17323fa0ddec8e887b4a7c53b8c0a0a220cfd000000000000000000")
            @test header.prevhash == want
            want = hex2bytes("5b0750fce0a889502d40508d39576821155e9c9e3f5c3157f961db38fd8b25be")
            @test header.merkleroot == want
            @test header.time == 0x59a7771e
            @test header.bits == 0x18013ce9
            @test header.nonce == 0x1dd7ffa4
        end
        @testset "Serialize" begin
            @test serialize(header) == header_raw
        end
        @testset "Hash" begin
            @test hash256(header) == hex2bytes("0000000000000000007e9e4c586439b0cdbe13b1370bdd9435d76a644d047523")
        end
        @testset "BIP9" begin
            @test BitcoinPrimitives.bip9(header)
            header_raw = hex2bytes("0400000039fa821848781f027a2e6dfabbf6bda920d9ae61b63400030000000000000000ecae536a304042e3154be0e3e9a8220e5568c3433a9ab49ac4cbb74f8df8e8b0cc2acf569fb9061806652c27")
            io = IOBuffer(header_raw)
            header = Header(io)
            @test !BitcoinPrimitives.bip9(header)
        end
        @testset "BIP91" begin
            @test !BitcoinPrimitives.bip91(header)
            header_raw = hex2bytes("1200002028856ec5bca29cf76980d368b0a163a0bb81fc192951270100000000000000003288f32a2831833c31a25401c52093eb545d28157e200a64b21b3ae8f21c507401877b5935470118144dbfd1")
            io = IOBuffer(header_raw)
            header = Header(io)
            @test BitcoinPrimitives.bip91(header)
        end
        @testset "BIP141" begin
            @test BitcoinPrimitives.bip141(header)
            header_raw = hex2bytes("0000002066f09203c1cf5ef1531f24ed21b1915ae9abeb691f0d2e0100000000000000003de0976428ce56125351bae62c5b8b8c79d8297c702ea05d60feabb4ed188b59c36fa759e93c0118b74b2618")
            io = IOBuffer(header_raw)
            header = Header(io)
            @test !BitcoinPrimitives.bip141(header)
        end
        @testset "Target" begin
            @test target(header) == parse(BigInt, "13ce9000000000000000000000000000000000000000000", base=16)
            @test difficulty(header) == 888171856257
        end
        @testset "Check POW" begin
            header_raw = hex2bytes("04000000fbedbbf0cfdaf278c094f187f2eb987c86a199da22bbb20400000000000000007b7697b29129648fa08b4bcd13c9d5e60abb973a1efac9c8d573c71c807c56c3d6213557faa80518c3737ec1")
            io = IOBuffer(header_raw)
            header = Header(io)
            @test check_pow(header) == true
            header_raw = hex2bytes("04000000fbedbbf0cfdaf278c094f187f2eb987c86a199da22bbb20400000000000000007b7697b29129648fa08b4bcd13c9d5e60abb973a1efac9c8d573c71c807c56c3d6213557faa80518c3737ec0")
            io = IOBuffer(header_raw)
            header = Header(io)
            @test !check_pow(header)
        end
    end
    raw_block = hex2bytes("0100000000000000000000000000000000000000000000000000000000000000000000003ba3edfd7a7b12b27ac72c3e67768f617fc81bc3888a51323a9fb8aa4b1e5e4a29ab5f49ffff001d1dac2b7c0101000000010000000000000000000000000000000000000000000000000000000000000000ffffffff4d04ffff001d0104455468652054696d65732030332f4a616e2f32303039204368616e63656c6c6f72206f6e206272696e6b206f66207365636f6e64206261696c6f757420666f722062616e6b73ffffffff0100f2052a01000000434104678afdb0fe5548271967f1a67130b7105cd6a828e03909a67962e0ea1f61deb649f6bc3f4cef38c4f35504e51ec112de5c384df7ba0b8d578a4c702b6bf11d5fac00000000")
    block = Block(IOBuffer(raw_block))
    @test block.header.version == 0x00000001
    @test block.header.prevhash == fill(0x00, 32)
    @test block.header.bits == 0x1d00ffff
    @test length(block.transactions) == 1
    @test block.transactions[1].inputs[1].prevout.index == 0xffffffff
    @test block.transactions[1].outputs[1].value == 0x000000012a05f200
end


text:
function spgenspd( n, p )
  nzval = Float64[]
  rowval = Int64[]
  colptr = Array(Int64,n+1)
  k = 1
  for j = 1:n
    colptr[j] = k
    push!(rowval, j)
    push!(nzval, randn() + n)
    k = k+1
    for row in randsubseq(1:j-1,p)
      push!(nzval, randn())
      push!(rowval, row)
      k = k+1
    end
  end
  colptr[n+1] = k

  return SparseMatrixCSC(n,n,colptr,rowval,nzval)
end

function genspd( n )
  A = randn(n,n)
  A = A+A' + n*eye(n)
  return A
end


text:
"""
    Fourier

A functor struct representing the Fourier basis.
"""
struct Fourier end

"""
    fourier_basis(k,x)

Evaluate the `k`-th Fourier basis function at `x`.
"""
function fourier_basis(k,x)
    cos((div(k,2)*x - iseven(k)/2)*π)
end

fourier_basis(k,x::AbstractVector) = fourier_basis.(k,x)

(::Fourier)(k,x) = fourier_basis(k,x)