text:
struct OneHotTensor{N, T} <: Block
    classes::AbstractVector{T}
end

function checkblock(block::OneHotTensor{N}, a::AbstractArray{T, M}) where {M, N, T}
    return N + 1 == M && last(size(a)) == length(block.classes)
end

mockblock(block::OneHotTensor{0}) = encode(
    OneHot(), Validation(), Label(block.classes), rand(block.classes))

function mockblock(block::OneHotTensor{N}) where N
    maskblock = Mask{N}(block.classes)
    return encode(OneHot(), Validation(), maskblock, mockblock(maskblock))
end

struct OneHotTensorMulti{N, T} <: Block
    classes::AbstractVector{T}
end

function checkblock(block::OneHotTensorMulti{N}, a::AbstractArray{T, M}) where {M, N, T}
    return N + 1 == M && last(size(a)) == length(block.classes)
end

function mockblock(block::OneHotTensorMulti{0}) where N
    labelblock = LabelMulti(block.classes)
    return encode(OneHot(), Validation(), labelblock, mockblock(labelblock))
end

"""
    OneHot()
    OneHot(T, threshold)

`Encoding` that turns categorical labels into one-hot encoded arrays of type `T`.

Encodes
```
      `Mask{N, U}` -> `OneHotTensor{N, T}`
`LabelMulti{N, U}` -> `OneHotTensorMulti{N, T}`
        `Label{U}` -> `OneHotTensor{N, T}`
```
"""
struct OneHot{TT<:Type} <: Encoding
    T::TT
    threshold::Float32
end

OneHot() = OneHot(Float32, 0.5f0)

# ### `Label` implementation
encodedblock(::OneHot, block::Label{T}) where T = OneHotTensor{0, T}(block.classes)
decodedblock(::OneHot, block::OneHotTensor{0}) = Label(block.classes)

function encode(enc::OneHot, context, block::Label, data)
    idx = findfirst(isequal(data), block.classes)
    isnothing(idx) && error("$data could not be found in `block.classes`: $(block.classes).")
    return DataAugmentation.onehot(enc.T, idx, length(block.classes))
end

function decode(::OneHot, context, block::OneHotTensor{0}, data)
    return block.classes[argmax(data)]
end

# ### `LabelMulti` implementation

encodedblock(::OneHot, block::LabelMulti{T}) where T = OneHotTensorMulti{0, T}(block.classes)
decodedblock(::OneHot, block::OneHotTensorMulti{0}) = LabelMulti(block.classes)

function encode(enc::OneHot, context, block::LabelMulti, data)
    return collect(enc.T, (c in data for c in block.classes))
end

function decode(enc::OneHot, context, block::OneHotTensorMulti{0}, data)
    return block.classes[softmax(data) .> enc.threshold]
end

# ### `Mask` implementation

encodedblock(::OneHot, block::Mask{N, T}) where {N, T} = OneHotTensor{N, T}(block.classes)
decodedblock(::OneHot, block::OneHotTensor{N, T}) where {N, T} = Mask{N, T}(block.classes)

function encode(enc::OneHot, context, block::Mask, data)
    tfm = DataAugmentation.OneHot{enc.T}()
    return apply(tfm, DataAugmentation.MaskMulti(data, block.classes)) |> DataAugmentation.itemdata
end

function decode(enc::OneHot, context, block::OneHotTensor, data)
    Tidx = length(block.classes) >= 255 ? UInt16 : UInt8
    classidxs = reshape(
        map(I -> Tidx(I.I[end]), argmax(data; dims = ndims(data))),
        size(data)[1:end-1])
    return IndirectArray(classidxs, block.classes)
end


text:
using Pkg
Pkg.activate(pwd())

using titanic

traindata, trainlabels = loadtrain("../data/train.csv")
testdata = loadtest("../data/test.csv")

alpha = svm(traindata, trainlabels)

prediction = classify(testdata, traindata, trainlabels, alpha)

savepredictions("../results/predictions-svm.csv", prediction)

ext: jl


text:
np = length(p);

function make_prob(T0, P, phi, p; tfinal=1.0)
    X0 = zeros(ns);
    X0[species_index(gas, fuel)] = phi
    X0[species_index(gas, oxygen)] = fuel2air
    X0[species_index(gas, inert)] = fuel2air * 3.76
    X0 = X0 ./ sum(X0);
    Y0 = X2Y(gas, X0, dot(X0, gas.MW));
    u0 = vcat(Y0, T0);

    mean_MW = 1.0 / dot(Y0, 1 ./ gas.MW)
    ρ_mass = P / R / T0 * mean_MW
    global ρ_mass

    prob = ODEProblem(dudt!, u0, (0.0, tfinal), p);
    prob
end

@inbounds function dudt!(du, u, p, t)
    Y = @view(u[1:ns])
    T = u[end]
    mean_MW = 1.0 / dot(Y, 1 ./ gas.MW)
    # ρ_mass = P / R / T * mean_MW
    P = ρ_mass * R * T / mean_MW
    X = Y2X(gas, Y, mean_MW)
    C = Y2C(gas, Y, ρ_mass)
    cp_mole, cp_mass = get_cp(gas, T, X, mean_MW)
    cv_mole, cv_mass = get_cv(cp_mole, cp_mass, mean_MW)
    h_mole = get_H(gas, T, Y, X)
    u_mole = get_U(h_mole, Y)
    S0 = get_S(gas, T, P, X)
    _p = reshape(p, npr, 3)
    # kp = exp.(p)
    qdot = wdot_func(gas.reaction, T, C, S0, h_mole; get_qdot=true)
    @. qdot[1:npr] *= @views(exp(_p[:, 1] + _p[:, 2] * log(T) - _p[:, 3] * 1000.0 * 4184.0 / R / T))
    wdot = gas.reaction.vk * qdot
    Ydot = wdot / ρ_mass .* gas.MW
    Tdot = -dot(u_mole, wdot) / ρ_mass / cv_mass
    du .= vcat(Ydot, Tdot)
    return du
end

@inbounds function dudtp!(du, u, p, t)
    Y = @view(u[1:ns])
    T = u[end]
    mean_MW = 1.0 / dot(Y, 1 ./ gas.MW)
    # ρ_mass = P / R / T * mean_MW
    P = ρ_mass * R * T / mean_MW
    X = Y2X(gas, Y, mean_MW)
    C = Y2C(gas, Y, ρ_mass)
    cp_mole, cp_mass = get_cp(gas, T, X, mean_MW)
    cv_mole, cv_mass = get_cv(cp_mole, cp_mass, mean_MW)
    h_mole = get_H(gas, T, Y, X)
    u_mole = get_U(h_mole, Y)
    S0 = get_S(gas, T, P, X)
    _p = reshape(p, npr, 3)
    # kp = exp.(p)
    qdot::typeof(p) = wdot_func(gas.reaction, T, C, S0, h_mole; get_qdot=true)
    @. qdot[1:npr] *= @views(exp(_p[:, 1] + _p[:, 2] * log(T) - _p[:, 3] * 1000.0 * 4184.0 / R / T))
    wdot = gas.reaction.vk * qdot
    Ydot = wdot / ρ_mass .* gas.MW
    Tdot = -dot(u_mole, wdot) / ρ_mass / cv_mass
    du .= vcat(Ydot, Tdot)
    return du
end

function sensBVP_mthreadp(ts, pred, p)
    idt = ts[end]
    Tign = pred[end, end]
    ng = length(ts)
    Fy = BandedMatrix(Zeros(ng * nu, ng * nu), (nu, nu));
    Fp = zeros(ng * nu, np)
    i = 1
    i_F = 1 + (i - 1) * nu:i * nu - 1
    @view(Fy[i_F, i_F])[ind_diag] .= ones_nu
    Fy[i * nu, i * nu] = -1.0
    Fy[i * nu, (i + 1) * nu] = 1.0

    dts = @views(ts[2:end] .- ts[1:end - 1]) ./ idt

    @threads for i = 2:ng
        u = @view(pred[:, i])
        du = similar(u)
        i_F = 1 + (i - 1) * nu:i * nu - 1
        @view(Fp[i_F, :]) .= jacobian((du, x) -> dudtp!(du, u, x, 0.0),
                                    du, p)::Array{Float64,2} .* (-idt)
    end
    @threads for i = 2:ng
        u = @view(pred[:, i])
        du = similar(u)
        i_F = 1 + (i - 1) * nu:i * nu - 1
        @view(Fy[i_F, i_F]) .= jacobian((du, x) -> dudt!(du, x, p, 0.0),
                                    du, u)::Array{Float64,2} .* (-idt)
    end

    for i = 2:ng
        u = @view(pred[:, i])
        du = similar(u)
        i_F = 1 + (i - 1) * nu:i * nu - 1
        @view(Fy[i_F, i * nu]) .= - dudt!(du, u, p, 0.0)
        @view(Fy[i_F, i_F])[ind_diag] .+= ones_nu ./ (dts[i - 1])
        @view(Fy[i_F, i_F .- nu])[ind_diag] .+= ones_nu ./ (-dts[i - 1])
        if i < ng
            Fy[i * nu, i * nu] = -1.0
            Fy[i * nu, (i + 1) * nu] = 1.0
        else
            Fy[i * nu, i * nu - 1] = 1.0
        end
    end
    dydp = - Fy \ sparse(Fp)
    return @view(dydp[end, :]) ./ idt
end

max_stars_repo_name: DENG-MIT/DeepReduction

max_issues_repo_name: DENG-MIT/DeepReduction

max_forks_repo_name: DENG-MIT/DeepReduction


text:
"""
    proxNet1(v, lambda)

Proximal operator for the Elastic-net penalisation component updates in ADMM.

# Arguments

- v = float; value to update
- lambdaL1 = l1 penalty, a floating scalar
- lambdaL2 = l1 penalty, a floating scalar

# Value

2d array of floats

"""

function proxNet1(v::Float64, lambda::Float64)

    return max(0.0, abs(v)-lambda) * sign(v)
end

"""
    proxNet2(v, rho, u, l)

Proximal operator for the residual sum of squares component updates in ADMM.

# Arguments

- v = float; value to update
- rho = float; parameter that controls ADMM tuning.
- u = float; corresponding transformed Y
- l = float; corresponding eigenvalue

# Value

2d array of floats

"""

function proxNet2(v::Float64, rho::Float64,
               u::Float64, l::Float64)

    return (u + rho*v) / (l + rho)
end

"""
    update_admm!(B, B0, B2, resid, X, Y, Z, Qx, Qz, U, L, lambdaL1, lambdaL2,
                 regXidx, regZidx, rho, r, s, tau_incr, tau_decr, mu)

Updates coefficient estimates in place for each ADMM iteration.

# Arguments

- B = 2d array of floats consisting of coefficient estimates for L1 updates
- B0 = 2d array of floats consisting of coefficient estimates for L2 updates
- B2 = 2d array of floats consisting of coefficient estimates for dual updates
- resid = 2d array of floats consisting of the residuals
- X = 2d array of floats consisting of the row covariates, with all
  categorical variables coded in appropriate contrasts
- Y = 2d array of floats consisting of the multivariate response
- Z = 2d array of floats consisting of the column covariates, with all
  categorical variables coded in appropriate contrasts
- Qx = 2d array of floats consisting of the eigenvectors of X
- Qz = 2d array of floats consisting of the eigenvectors of Z
- U = 2d array of floats consisting of the transformed Y matrix
- L = 2d array of floats consisting of the kronecker product of the
  eigenvalues of X and Z
- lambdaL1 = l1 penalty, a floating scalar
- lambdaL2 = l2 penalty, a floating scalar
- regXidx = 1d array of indices corresponding to regularized X covariates
- regZidx = 1d array of indices corresponding to regularized Z covariates
- rho = float; parameter that controls ADMM tuning.
- r = 2d array of floats consisting of the primal residuals.
- s = 2d array of floats consisting of the dual residuals.
- tau_incr = float; parameter that controls the factor at which rho increases.
  Defaults to 2.0.
- tau_decr = float; parameter that controls the factor at which rho decreases.
  Defaults to 2.0.
- mu = float; parameter that controls the factor at which the primal and dual
  residuals should be within each other. Defaults to 10.0.

# Value

None; updates coefficients in place

# Some notes

Convergence is determined as when the log ratio of the current and previous
criteria is less than the threshold `thresh`.

`rho` controls ADMM tuning and can be specified by the user.

"""
function update_admm!(B::AbstractArray{Float64,2},
                      B0::AbstractArray{Float64,2},
                      B2::AbstractArray{Float64,2},
                      resid::AbstractArray{Float64,2},
                      X::AbstractArray{Float64,2},
                      Y::AbstractArray{Float64,2},
                      Z::AbstractArray{Float64,2},
                      Qx::AbstractArray{Float64,2},
                      Qz::AbstractArray{Float64,2},
                      U::AbstractArray{Float64,2},
                      L::AbstractArray{Float64,2},
                      lambdaL1::Float64, lambdaL2::Float64,
                      regXidx::AbstractArray{Int64,1},
                      regZidx::AbstractArray{Int64,1},
                      rho::AbstractArray{Float64,1},
                      r::AbstractArray{Float64,2},
                      s::AbstractArray{Float64,2},
                      tau_incr::Float64, tau_decr::Float64, mu::Float64)

    # Save a copy of the old values of B
    s .= copy(B)

    # L2 updates
    B0 .= B .- B2
    # Transform B0
    mul!(B0, transpose(Qx), B0*Qz)
    # Perform L2 updates and transform B0 back
    mul!(B0, Qx * proxNet2.(B0, rho[1], U, L), transpose(Qz))

    # L1 updates
    B .= B0 .+ B2
    B[regXidx,regZidx] .= proxNet1.(B[regXidx,regZidx], lambdaL1/rho[1])/(1+lambdaL2/rho[1])

    # Primal residuals
    r .= B0 .- B
    # Dual residuals
    s .-= B

    # Dual updates
    B2 .+= r

    # Update residuals
    calc_resid!(resid, X, Y, Z, B)

    # Update rho according to the relative size of the primal and dual residuals
    # Then re-scale B2 accordingly
    if sqrt(sum(r.^2)) > mu * rho[1] * sqrt(sum(s.^2))
        rho[:] .*= tau_incr
        B2 ./= tau_incr
    elseif rho[1] * sqrt(sum(s.^2)) > mu * sqrt(sum(r.^2))
        rho[:] ./= tau_decr
        B2 .*= tau_decr
    end
end

"""
    admm!(X, Y, Z, lambda, alpha, B, regXidx, regZidx, reg, norms, Qx, Qz, U, L;
          isVerbose, stepsize, rho, setRho, thresh, maxiter,
          tau_incr, tau_decr, mu)

Performs ADMM.

# Arguments

- X = 2d array of floats consisting of the row covariates, with all
  categorical variables coded in appropriate contrasts
- Y = 2d array of floats consisting of the multivariate response
- Z = 2d array of floats consisting of the column covariates, with all
  categorical variables coded in appropriate contrasts
- lambda = lambda penalty, a floating scalar
- alpha = parameter (ϵ[0, 1]) determining the mix of penalties between L1 and L2
- B = 2d array of floats consisting of starting coefficient estimates
- regXidx = 1d array of indices corresponding to regularized X covariates
- regZidx = 1d array of indices corresponding to regularized Z covariates
- reg = 2d array of bits, indicating whether or not to regularize each of the
  coefficients
- norms = 2d array of floats consisting of the norms corresponding to each
  coefficient or `nothing`
- Qx = 2d array of floats consisting of the eigenvectors of X
- Qz = 2d array of floats consisting of the eigenvectors of Z
- U = 2d array of floats consisting of the transformed Y matrix
- L = 2d array of floats consisting of the kronecker product of the
  eigenvalues of X and Z

# Keyword arguments

- isVerbose = boolean flag indicating whether or not to print messages.
  Defaults to `true`.
- stepsize = float; step size for updates (irrelevant for ADMM).
  Defaults to `0.01`.
- rho = float; parameter that controls ADMM tuning. Defaults to `1.0`.
- setRho = boolean flag indicating whether the ADMM tuning parameter `rho`
  should be calculated. Defaults to `true`.
- thresh = threshold at which the coefficients are considered to have
  converged, a floating scalar. Defaults to `10^(-7)`.
- maxiter = maximum number of update iterations. Defaults to `10^10`.
- tau_incr = float; parameter that controls the factor at which rho increases.
  Defaults to 2.0.
- tau_decr = float; parameter that controls the factor at which rho decreases.
  Defaults to 2.0.
- mu = float; parameter that controls the factor at which the primal and dual
  residuals should be within each other. Defaults to 10.0.

# Value

None; updates coefficients in place

# Some notes

Convergence is determined as when the log ratio of the current and previous
criteria is less than the threshold `thresh`.

`rho` controls ADMM tuning and can be specified by the user.

"""
function admm!(X::AbstractArray{Float64,2}, Y::AbstractArray{Float64,2},
               Z::AbstractArray{Float64,2}, lambda::Float64, alpha::Float64,
               B::AbstractArray{Float64,2},
               regXidx::AbstractArray{Int64,1},
               regZidx::AbstractArray{Int64,1}, reg::BitArray{2}, norms,
               Qx::AbstractArray{Float64,2}, Qz::AbstractArray{Float64,2},
               U::AbstractArray{Float64,2}, L::AbstractArray{Float64,2};
               isVerbose::Bool=true, stepsize::Float64=0.01,
               rho::Float64=1.0, setRho::Bool=true,
               thresh::Float64=10.0^(-7), maxiter::Int=10^10,
               tau_incr::Float64=2.0, tau_decr::Float64=2.0, mu::Float64=10.0)

    # Re-parametrize the Elastic-net tuning parameters
    lambdaL1 = lambda*alpha
    lambdaL2 = lambda*(1-alpha)

    # Set the ADMM tuning parameter, rho
    if setRho == true
        # Get smallest and largest eigenvalues
        mineig = minimum(L)
        maxeig = maximum(L)

        gamma = lambdaL1/sqrt(1+lambdaL2)
        # Set the value of rho
        if gamma < mineig
            rho = mineig
        elseif gamma > maxeig
            rho = gamma
        else
            rho = maxeig
        end
    end
    rho = [rho]

    # Calculate residuals
    resid = calc_resid(X, Y, Z, B)
    # Initialize values for L2 updates
    B0 = copy(B)
    # Initialize values for dual updates
    B2 = copy(B)

    # Initial primal residuals
    r = copy(B)
    # Initialize dual residuals
    s = copy(B)

    # Denominators of criterion
    crit_denom = [size(X,1)*size(Z,1), size(X,2)*size(Z,2)]
    # Placeholder to store the old criterion
    oldcrit = 1.0
    # Calculate the current criterion
    crit = criterion(B[regXidx, regZidx], resid, lambdaL1, lambdaL2, crit_denom)

    iter = 0
    # Iterate until coefficients converge or maximum iterations have been
    # reached.
    while (abs(log(crit/oldcrit)) > thresh) && (iter < maxiter)
        # Store the current criterion
        oldcrit = crit

        # Update the coefficient estimates
        update_admm!(B, B0, B2, resid, X, Y, Z, Qx, Qz, U, L, lambdaL1, lambdaL2,
                     regXidx, regZidx, rho, r, s, tau_incr, tau_incr, mu)

        # Calculate the criterion after updating
        crit = criterion(B[regXidx, regZidx], resid, lambdaL1, lambdaL2, crit_denom)

        iter += 1 # Increment the number of iterations
        # Print warning message if coefficient estimates do not converge.
        if iter==(maxiter)
            println("Estimates did not converge.")
        end
    end

    println_verbose(string("Criterion: ", crit), isVerbose)
    println_verbose(string("Number of iterations: ", iter), isVerbose)
end


text:
using Symbolics, SymbolicUtils
using ModelingToolkit, DifferentialEquations
using Plots

############ The Helper Functions ############################################

def_taylor(x, ps) = sum([a*x^i for (i,a) in enumerate(ps)])
def_taylor(x, ps, p₀) = p₀ + def_taylor(x, ps)

expand_sin(x, n) = sum([(isodd(k) ? -1 : 1)*(-x)^(2k-1)/factorial(2k-1) for k=1:n])

function collect_powers(eq, x, ns; max_power=100)
    eq = substitute(expand(eq), Dict(x^j => 0 for j=last(ns)+1:max_power))

    eqs = []
    for i in ns
        powers = Dict(x^j => (i==j ? 1 : 0) for j=1:last(ns))
        push!(eqs, substitute(eq, powers))
    end
    eqs
end

function solve_coef(eqs, ps)
    vals = Dict()

    for i = 1:length(ps)
        eq = substitute(eqs[i], vals)
        vals[ps[i]] = Symbolics.solve_for(eq ~ 0, ps[i])
    end
    vals
end

function solve_newton(f, x, x₀; abstol=1e-8, maxiter=50)
    xₙ = Float64(x₀)
    fₙ₊₁ = x - f / Symbolics.derivative(f, x)

    for i = 1:maxiter
        xₙ₊₁ = substitute(fₙ₊₁, Dict(x => xₙ))
        if abs(xₙ₊₁ - xₙ) < abstol
            return xₙ₊₁
        else
            xₙ = xₙ₊₁
        end
    end
    return xₙ₊₁
end

function calc_u0(sys, d; def_val=0.0)
    S = states(sys)
    u0 = Dict{typeof(S[1]), Float64}()

    for s in S
        name_withˍt = string(s)
        name_withoutˍt = replace(name_withˍt, r"\(.*\)" => "")

        if haskey(d, name_withˍt)
            u0[s] = d[name_withˍt]
        elseif haskey(d, name_withoutˍt)
            u0[s] = d[name_withoutˍt]
        else
            u0[s] = def_val
        end
    end

    return u0
end

#################### The examples ############################################

"""
    x^5 + x = 1
"""
function test_quintic(n=4)
    @variables ϵ z a[1:n]
    x = def_taylor(ϵ, a, 1)
    eq = x^5 + ϵ*x - 1
    eqs = collect_powers(eq, ϵ, 1:n)
    vals = solve_coef(eqs, a)

    x′ = substitute(x, vals)
    X = 𝜀 -> substitute(x′, Dict(ϵ => 𝜀))

    xₚ = X(1.0)
    xₙ = solve_newton(z^5 + z - 1, z, 1.0)

    return xₚ, xₙ
end

"""
    E - e * sin(E) = M
"""
function test_kepler(n=4)
    @variables ϵ z M a[1:n]
    x = def_taylor(ϵ, a, M)
    eq = x - ϵ * expand_sin(x, n) - M
    eqs = collect_powers(eq, ϵ, 1:n)
    vals = solve_coef(eqs, a)

    x′ = substitute(x, vals)
    X = (𝜀, 𝑀) -> substitute(x′, Dict(ϵ => 𝜀, M => 𝑀))

    𝑒 = 0.01671    # The Earth eccentricity
    M₀ = π/2

    xₚ = X(𝑒, M₀)
    xₙ = solve_newton(z - 𝑒*sin(z) - M₀, z, M₀)

    return xₚ, xₙ
end

function test_rocket(n=3)
    sys, y = calc_rocket_sys(n)
    u0 = calc_u0(sys, Dict("y₁ˍt" => 1.0))
    prob = ODEProblem(sys, u0, (0, 3.0))
    sol = solve(prob; dtmax=0.01)
    X = 𝜀 -> sum([𝜀^(i-1) * sol[y[i]] for i in eachindex(y)])
    plot(sol.t, hcat([X(ϵ) for ϵ = 0.0:0.1:0.5]...))
end

function calc_rocket_sys(n)
    @variables ϵ t y[1:n](t) ∂∂y[1:n]
    x = def_taylor(ϵ, y[2:end], y[1])
    ∂∂x = def_taylor(ϵ, ∂∂y[2:end], ∂∂y[1])
    eq = ∂∂x * (1 + ϵ*x)^2 + 1
    eqs = collect_powers(eq, ϵ, 0:n)
    vals = solve_coef(eqs, ∂∂y)

    D = Differential(t)
    subs = Dict(∂∂y[i] => D(D(y[i])) for i in eachindex(y))
    eqs = [substitute(first(v), subs) ~ substitute(last(v), subs) for v in vals]

    sys = ODESystem(eqs, t)
    sys = ode_order_lowering(sys)

    sys, y
end

function test_oscillator(n=3)
    sys, y = calc_oscillator_sys(n)
    u0 = calc_u0(sys, Dict("y₁ˍt" => 1.0))
    prob = ODEProblem(sys, u0, (0, 50.0))
    sol = solve(prob; dtmax=0.01)

    X = 𝜀 -> sum([𝜀^(i-1) * sol[y[i]] for i in eachindex(y)])
    T = sol.t
    Y = 𝜀 -> exp.(-𝜀*T) .* sin.(sqrt(1 - 𝜀^2)*T) / sqrt(1 - 𝜀^2)    # exact solution

    plot(sol.t, [Y(0.1), X(0.1)])
end

function calc_oscillator_sys(n)
    @variables ϵ t y[1:n](t) ∂y[1:n] ∂∂y[1:n]
    x = def_taylor(ϵ, y[2:end], y[1])
    ∂x = def_taylor(ϵ, ∂y[2:end], ∂y[1])
    ∂∂x = def_taylor(ϵ, ∂∂y[2:end], ∂∂y[1])

    eq = ∂∂x + 2*ϵ*∂x + x
    eqs = collect_powers(eq, ϵ, 0:n)
    vals = solve_coef(eqs, ∂∂y)

    D = Differential(t)
    subs1 = Dict(∂y[i] => D(y[i]) for i in eachindex(y))
    subs2 = Dict(∂∂y[i] => D(D(y[i])) for i in eachindex(y))
    subs = subs1 ∪ subs2
    eqs = [substitute(first(v), subs) ~ substitute(last(v), subs) for v in vals]

    sys = ODESystem(eqs, t)
    sys = ode_order_lowering(sys)

    sys, y
end


text:
using brainflow

# specify logging library to use
brainflow.enable_dev_brainflow_logger(brainflow.BOARD_CONTROLLER)

params = BrainFlowInputParams()
board_shim = brainflow.BoardShim(brainflow.SYNTHETIC_BOARD, params)

brainflow.prepare_session(board_shim)
brainflow.start_stream(board_shim)
sleep(5)
brainflow.stop_stream(board_shim)
data = brainflow.get_current_board_data(32, board_shim)
brainflow.release_session(board_shim)

eeg_channels = brainflow.get_eeg_channels(brainflow.SYNTHETIC_BOARD)
sampling_rate = brainflow.get_sampling_rate(brainflow.SYNTHETIC_BOARD)

data_first_channel = data[eeg_channels[1], :]
println("Original Data First Channel")
println(data_first_channel)
downsampled_data = brainflow.perform_downsampling(data_first_channel, 3, brainflow.EACH)
println("After Downsampling")
println(downsampled_data)

max_stars_repo_name: matthijscox/brainflow

max_issues_repo_name: matthijscox/brainflow

max_forks_repo_name: matthijscox/brainflow


text:
function get_area(obj::PropertiesComputer)
    return jcall(obj, "getArea", jdouble, ())
end

function get_barycenter(obj::PropertiesComputer)
    return jcall(obj, "getBarycenter", S2Point, ())
end

function get_convex_cells_inside_points(obj::PropertiesComputer)
    return jcall(obj, "getConvexCellsInsidePoints", List, ())
end

function visit_internal_node(obj::PropertiesComputer, arg0::BSPTree)
    return jcall(obj, "visitInternalNode", void, (BSPTree,), arg0)
end

function visit_leaf_node(obj::PropertiesComputer, arg0::BSPTree)
    return jcall(obj, "visitLeafNode", void, (BSPTree,), arg0)
end

function visit_order(obj::PropertiesComputer, arg0::BSPTree)
    return jcall(obj, "visitOrder", BSPTreeVisitor_Order, (BSPTree,), arg0)
end


text:
# this stuff pending PR into ProgressMeter.jl, for now copied here

using Serialization
using ProgressMeter: AbstractProgress
import ProgressMeter: next!, update!, finish!, cancel

"""
    RemoteProgressWrapper(pbar::AbstractProgress)

Wraps any `AbstractProgress` and returns another `AbstractProgress`
which is safe to ship to remote workers. Calls to `next!`, `update!`,
`cancel`, and `finish!` on the remote workers will be forwarded to the
master process and will be asynchronously applied to the original
`pbar` object.
"""
struct RemoteProgressWrapper <: AbstractProgress
    channel
    pbar
    function RemoteProgressWrapper(pbar::AbstractProgress)
        channel = RemoteChannel(()->Channel{Any}(), 1)
        @async while true
            msg = take!(channel)
            isnothing(msg) && break
            (func, args, kwargs) = msg
            func(pbar, args...; kwargs...)
        end
        new(channel, pbar)
    end
    RemoteProgressWrapper(channel::RemoteChannel) = new(channel)
end

next!(pbar::RemoteProgressWrapper, args...; kwargs...)   =  put!(pbar.channel, (next!,   args, kwargs))
update!(pbar::RemoteProgressWrapper, args...; kwargs...) =  put!(pbar.channel, (update!, args, kwargs))
finish!(pbar::RemoteProgressWrapper, args...; kwargs...) = (put!(pbar.channel, (finish!, args, kwargs)); close(pbar.channel))
cancel(pbar::RemoteProgressWrapper, args...; kwargs...)  = (put!(pbar.channel, (cancel,  args, kwargs)); close(pbar.channel))

function Serialization.serialize(s::AbstractSerializer, pbar::RemoteProgressWrapper)
    # only serializes the channel, remote workers never receive the
    # `pbar` field (which isn't needed / can't be serialized)
    Serialization.writetag(s.io, Serialization.OBJECT_TAG)
    Serialization.serialize(s, RemoteProgressWrapper)
    Serialization.serialize(s, pbar.channel)
end
function Serialization.deserialize(s::AbstractSerializer, ::Type{RemoteProgressWrapper})
    RemoteProgressWrapper(Serialization.deserialize(s))
end

RemoteProgress(args...; kwargs...) = RemoteProgressWrapper(Progress(args...; kwargs...))
RemoteProgressThresh(args...; kwargs...) = RemoteProgressWrapper(ProgressThresh(args...; kwargs...))
RemoteProgressUnknown(args...; kwargs...) = RemoteProgressWrapper(ProgressUnknown(args...; kwargs...))


text:
#= MIT License

Copyright (c) 2020, 2021, 2022 Uwe Fechner

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE. =#

#= Models of a kite-power system in implicit form: residual = f(y, yd)

This model implements a 3D mass-spring system with reel-out. It uses six tether segments (the number can be
configured in the file data/settings.yaml). Two kite models are provided, the one point and the four point
kite model. The spring constant and the damping decrease with the segment length. The aerodynamic kite forces are
calculated, depending on reel-out speed, depower and steering settings.

Scientific background: http://arxiv.org/abs/1406.6218 =#

module KiteModels

using Dierckx, StaticArrays, Rotations, LinearAlgebra, Parameters, NLsolve, DocStringExtensions, Sundials
using Reexport
@reexport using KitePodModels
@reexport using WinchModels
@reexport using AtmosphericModels
import Base.zero
import KiteUtils.calc_elevation
import KiteUtils.calc_azimuth
import KiteUtils.calc_heading
import KiteUtils.calc_course
import KiteUtils.SysState

export KPS3, KPS4, KVec3, SimFloat, ProfileLaw, EXP, LOG, EXPLOG                                  # constants and types
export calc_set_cl_cd!, copy_examples, copy_bin                                                   # helper functions
export clear!, find_steady_state!, residual!                                                      # low level worker functions
export init_sim!, next_step!                                                                      # hight level worker functions
export pos_kite, calc_height, calc_elevation, calc_azimuth, calc_heading, calc_course             # getters
export winch_force, lift_drag, lift_over_drag, unstretched_length, tether_length, v_wind_kite     # getters
export kite_ref_frame, orient_euler, spring_forces

set_zero_subnormals(true)           # required to avoid drastic slow down on Intel CPUs when numbers become very small
KiteUtils.set_data_path("")         # this statement is only executed during precompilation and ensures that the default settings.yaml
                                    # are used

# Constants
const G_EARTH = 9.81                # gravitational acceleration
const BRIDLE_DRAG = 1.1             # should probably be removed

# Type definitions
"""
    const SimFloat = Float64

This type is used for all real variables, used in the Simulation. Possible alternatives: Float32, Double64, Dual
Other types than Float64 or Float32 do require support of Julia types by the solver.
"""
const SimFloat = Float64

"""
   const KVec3    = MVector{3, SimFloat}

Basic 3-dimensional vector, stack allocated, mutable.
"""
const KVec3    = MVector{3, SimFloat}

"""
   const SVec3    = SVector{3, SimFloat}

Basic 3-dimensional vector, stack allocated, immutable.
"""
const SVec3    = SVector{3, SimFloat}

# the following two definitions speed up the function residual! from 940ns to 540ns
# disadvantage: changing the cl and cd curves requires a restart of the program
const calc_cl = Spline1D(se().alpha_cl, se().cl_list)
const calc_cd = Spline1D(se().alpha_cd, se().cd_list)

"""
    abstract type AbstractKiteModel

All kite models must inherit from this type. All methods that are defined on this type must work
with all kite models. All exported methods must work on this type.
"""
abstract type AbstractKiteModel end

"""
    const AKM = AbstractKiteModel

Short alias for the AbstractKiteModel.
"""
const AKM = AbstractKiteModel

include("KPS4.jl") # include code, specific for the four point kite model
include("KPS3.jl") # include code, specific for the one point kite model
include("init.jl") # functions to calculate the inital state vector, the inital masses and initial springs

# Calculate the lift and drag coefficient as a function of the angle of attack alpha.
function set_cl_cd!(s::AKM, alpha)
    angle =  alpha * 180.0 / π
    if angle > 180.0
        angle -= 360.0
    end
    if angle < -180.0
        angle += 360.0
    end
    s.param_cl = calc_cl(angle)
    s.param_cd = calc_cd(angle)
    nothing
end

# Calculate the angle of attack alpha from the apparend wind velocity vector
# v_app and the z unit vector of the kite reference frame.
function calc_alpha(v_app, vec_z)
    π/2.0 - acos(-(v_app ⋅ vec_z) / norm(v_app))
end

"""
    set_depower_steering!(s::AKM, depower, steering)

Setter for the depower and steering model inputs.

Parameters:
- depower:   Relative depower,  must be between 0 .. 1.0
- steering:  Relative steering, must be between -1.0 .. 1.0.

This function sets the variables s.depower, s.steering and s.alpha_depower.

It takes the depower offset c0 and the dependency of the steering sensitivity from
the depower settings into account.
"""
function set_depower_steering!(s::AKM, depower, steering)
    s.depower  = depower
    s.alpha_depower = calc_alpha_depower(s.kcu, depower)
    s.steering = (steering - s.set.c0) / (1.0 + s.set.k_ds * (s.alpha_depower / deg2rad(s.set.alpha_d_max)))
    nothing
end

"""
    set_v_reel_out!(s::AKM, v_reel_out, t_0, period_time = 1.0 / s.set.sample_freq)

Setter for the reel-out speed. Must be called on every timestep (before each simulation).
It also updates the tether length, therefore it must be called even if `v_reel_out` has
not changed.

- t_0 the start time of the next timestep relative to the start of the simulation [s]
"""
function set_v_reel_out!(s::KPS3, v_reel_out, t_0, period_time = 1.0 / s.set.sample_freq)
    s.l_tether += 0.5 * (v_reel_out + s.last_v_reel_out) * period_time
    s.last_v_reel_out = s.v_reel_out
    s.v_reel_out = v_reel_out
    s.t_0 = t_0
end

function set_v_reel_out!(s::KPS4, v_reel_out, t_0, period_time = 1.0 / s.set.sample_freq)
    s.sync_speed = v_reel_out
    s.last_v_reel_out = s.v_reel_out
    s.t_0 = t_0
end

"""
    unstretched_length(s::AKM)

Getter for the unstretched tether reel-out lenght (at zero force).
"""
function unstretched_length(s::AKM) s.l_tether end

"""
    lift_drag(s::AKM)

Return a tuple of the scalar lift and drag forces.

**Example:**

    lift, drag = lift_drag(s)
"""
function lift_drag(s::AKM) return (norm(s.lift_force), norm(s.drag_force)) end

"""
    lift_over_drag(s::AKM)

Return the lift-over-drag ratio.
"""
function lift_over_drag(s::AKM)
    lift, drag = lift_drag(s)
    return lift / drag
end

"""
    v_wind_kite(s::AKM)

Return the vector of the wind speed at the height of the kite.
"""
function v_wind_kite(s::AKM) s.v_wind end

"""
    set_v_wind_ground!(s::AKM, height, v_wind_gnd=s.set.v_wind, wind_dir=0.0)

Set the vector of the wind-velocity at the height of the kite. As parameter the height,
the ground wind speed [m/s] and the wind direction [radians] are needed.
Must be called every at each timestep.
"""
function set_v_wind_ground!(s::AKM, height, v_wind_gnd=s.set.v_wind, wind_dir=0.0)
    if height < 6.0
        height = 6.0
    end
    s.v_wind .= v_wind_gnd * calc_wind_factor(s.am, height) .* [cos(wind_dir), sin(wind_dir), 0]
    s.v_wind_gnd .= [v_wind_gnd * cos(wind_dir), v_wind_gnd * sin(wind_dir), 0.0]
    s.v_wind_tether .= v_wind_gnd * calc_wind_factor(s.am, height / 2.0) .* [cos(wind_dir), sin(wind_dir), 0]
    s.rho = calc_rho(s.am, height)
    nothing
end

"""
    tether_length(s::AKM)

Calculate and return the real, stretched tether lenght.
"""
function tether_length(s::AKM)
    length = 0.0
    for i in 1:s.set.segments
        length += norm(s.pos[i+1] - s.pos[i])
    end
    return length
end

"""
    orient_euler(s::AKM)

Calculate and return the orientation of the kite in euler angles (roll, pitch, yaw)
as SVector.
"""
function orient_euler(s::AKM)
    x, y, z = kite_ref_frame(s)
    roll = atan(y[3], z[3]) - π/2
    if roll < -π/2
       roll += 2π
    end
    pitch = asin(-x[3])
    yaw = -atan(x[2], x[1]) - π/2
    if yaw < -π/2
        yaw += 2π
    end
    SVector(roll, pitch, yaw)
end

"""
    calc_elevation(s::AKM)

Determine the elevation angle of the kite in radian.
"""
function calc_elevation(s::AKM)
    KiteUtils.calc_elevation(pos_kite(s))
end

"""
    calc_azimuth(s::AKM)

Determine the azimuth angle of the kite in radian.
"""
function calc_azimuth(s::AKM)
    KiteUtils.azimuth_east(pos_kite(s))
end

"""
    calc_heading(s::AKM)

Determine the heading angle of the kite in radian.
"""
function calc_heading(s::AKM)
    orientation = orient_euler(s)
    elevation = calc_elevation(s)
    azimuth = calc_azimuth(s)
    KiteUtils.calc_heading(orientation, elevation, azimuth)
end

"""
    calc_course(s::AKM)

Determine the course angle of the kite in radian.
Undefined if the velocity of the kite is near zero.
"""
function calc_course(s::AKM)
    elevation = calc_elevation(s)
    azimuth = calc_azimuth(s)
    KiteUtils.calc_course(s.vel_kite, elevation, azimuth)
end

"""
    SysState(s::AKM, zoom=1.0)

Constructor for creating a SysState object from a kite model (KPS3 or KPS4).
The SysState object can be used either for logging or for displaying the
system state in a viewer. Optionally the position arrays can be zoomed
according to the requirements of the viewer.
"""
function SysState(s::AKM, zoom=1.0)
    pos = s.pos
    P = length(pos)
    X = zeros(MVector{P, MyFloat})
    Y = zeros(MVector{P, MyFloat})
    Z = zeros(MVector{P, MyFloat})
    for i in 1:P
        X[i] = pos[i][1] * zoom
        Y[i] = pos[i][2] * zoom
        Z[i] = pos[i][3] * zoom
    end

    x, y, z = kite_ref_frame(s)
    pos_kite_ = pos_kite(s)
    pos_before = pos_kite_ + z

    rotation = rot(pos_kite_, pos_before, -x)
    q = QuatRotation(rotation)
    orient = MVector{4, Float32}(Rotations.params(q))

    elevation = calc_elevation(s)
    azimuth = calc_azimuth(s)
    force = winch_force(s)
    heading = calc_heading(s)
    course = calc_course(s)
    v_app_norm = norm(s.v_apparent)
    KiteUtils.SysState{P}(s.t_0, orient, elevation, azimuth, s.l_tether, s.v_reel_out, force, s.depower, s.steering, heading, course, v_app_norm, s.vel_kite, X, Y, Z)
end

function calc_pre_tension(s::AKM)
    forces = spring_forces(s)
    av_force = 0.0
    for i in 1:s.set.segments
        av_force += forces[i]
    end
    av_force /= s.set.segments
    res = av_force/s.set.c_spring
    if res < 0.0 res = 0.0 end
    if isnan(res) res = 0.0 end
    return res + 1.0
end

"""
    init_sim!(s; t_end=1.0, stiffness_factor=0.035, prn=false)

Initialises the integrator of the model.

Parameters:
- s:     an instance of an abstract kite model
- t_end: end time of the simulation; normally not needed
- stiffness_factor: factor applied to the tether stiffness during initialisation
- prn: if set to true, print the detailed solver results

Returns:
An instance of a DAE integrator.
"""
function init_sim!(s::AKM; t_end=1.0, stiffness_factor=0.035, prn=false)
    clear!(s)
    s.stiffness_factor = stiffness_factor
    y0, yd0 = KiteModels.find_steady_state!(s; stiffness_factor=stiffness_factor, prn=prn)

    differential_vars = ones(Bool, length(y0))
    solver  = IDA(linear_solver=Symbol(s.set.linear_solver), max_order = s.set.max_order)
    tspan   = (0.0, t_end)
    abstol  = s.set.abs_tol # max error in m/s and m
    prob    = DAEProblem(residual!, yd0, y0, tspan, s, differential_vars=differential_vars)
    integrator = Sundials.init(prob, solver, abstol=abstol, reltol=s.set.rel_tol)
end

"""
    next_step!(s::AKM, integrator; v_ro = 0.0, v_wind_gnd=s.set.v_wind, wind_dir=0.0, dt=1/s.set.sample_freq)

Calculates the next simulation step.

Parameters:
- s:            an instance of an abstract kite model
- integrator:   an integrator instance as returned by the function [`init_sim!`](@ref)
- v_ro:         set value of reel out speed in m/s
- `v_wind_gnd`: wind speed at reference height in m/s
- wind_dir:     wind direction in radians
- dt:           time step in seconds

Only the first two parameters are required.

Returns:
The end time of the time step in seconds.
"""
function next_step!(s::AKM, integrator; v_ro = 0.0, v_wind_gnd=s.set.v_wind, wind_dir=0.0, dt=1/s.set.sample_freq)
    KitePodModels.on_timer(s.kcu)
    KiteModels.set_depower_steering!(s, get_depower(s.kcu), get_steering(s.kcu))
    set_v_reel_out!(s, v_ro, integrator.t)
    set_v_wind_ground!(s, calc_height(s), v_wind_gnd, wind_dir)
    Sundials.step!(integrator, dt, true)
    if s.stiffness_factor < 1.0
        s.stiffness_factor+=0.01
        if s.stiffness_factor > 1.0
            s.stiffness_factor = 1.0
        end
    end
    integrator.t
end

"""
    copy_examples()

Copy the example scripts to the folder "examples"
(it will be created if it doesn't exist).
"""
function copy_examples()
    PATH = "examples"
    if ! isdir(PATH)
        mkdir(PATH)
    end
    src_path = joinpath(dirname(pathof(@__MODULE__)), "..", PATH)
    cp(joinpath(src_path, "compare_kps3_kps4.jl"), joinpath(PATH, "compare_kps3_kps4.jl"), force=true)
    cp(joinpath(src_path, "plot2d.jl"), joinpath(PATH, "plot2d.jl"), force=true)
    cp(joinpath(src_path, "simulate.jl"), joinpath(PATH, "simulate.jl"), force=true)
    chmod(joinpath(PATH, "compare_kps3_kps4.jl"), 0o664)
    chmod(joinpath(PATH, "plot2d.jl"), 0o664)
    chmod(joinpath(PATH, "simulate.jl"), 0o664)
end

"""
    copy_bin()

Copy the scripts create_sys_image and run_julia to the folder "bin"
(it will be created if it doesn't exist).
"""
function copy_bin()
    PATH = "bin"
    if ! isdir(PATH)
        mkdir(PATH)
    end
    src_path = joinpath(dirname(pathof(@__MODULE__)), "..", PATH)
    cp(joinpath(src_path, "create_sys_image2"), joinpath(PATH, "create_sys_image"), force=true)
    cp(joinpath(src_path, "run_julia2"), joinpath(PATH, "run_julia"), force=true)
    chmod(joinpath(PATH, "create_sys_image"), 0o774)
    chmod(joinpath(PATH, "run_julia"), 0o774)
    PATH = "test"
    if ! isdir(PATH)
        mkdir(PATH)
    end
    src_path = joinpath(dirname(pathof(@__MODULE__)), "..", PATH)
    cp(joinpath(src_path, "create_sys_image2.jl"), joinpath(PATH, "create_sys_image.jl"), force=true)
    cp(joinpath(src_path, "test_for_precompile.jl"), joinpath(PATH, "test_for_precompile.jl"), force=true)
    cp(joinpath(src_path, "update_packages.jl"), joinpath(PATH, "update_packages.jl"), force=true)
    chmod(joinpath(PATH, "create_sys_image.jl"), 0o664)
    chmod(joinpath(PATH, "test_for_precompile.jl"), 0o664)
    chmod(joinpath(PATH, "update_packages.jl"), 0o664)
end

precompile(find_steady_state!, (KPS3{SimFloat, KVec3, se().segments+1},))
const particles = se().segments+KITE_PARTICLES+1
const springs = se().segments+KITE_SPRINGS
precompile(find_steady_state!, (KPS4{Float64, MVector{3, Float64}, particles, springs, KiteModels.Spring{Int16, Float64}},))
precompile(init_sim!, (KPS4{Float64, MVector{3, Float64}, particles, springs, KiteModels.Spring{Int16, Float64}}, Float64,))

end


text:
# This is an example code used to precompile the SpinLab GUI
# Using a precompiled version of MRIsim should accelerate startup times for each library
# In the future, this will enable to create an Standalone-App (exceutable)

using MRIsim
SpinLab()

##
phantom = MRIsim.brain_phantom2D(;axis="coronal")
EPI,_,_,_ = MRIsim.EPI_base(40/100, 100, 4e-6, 60e-3)
TE = 25e-3
d = MRIsim.delay(TE-MRIsim.dur(EPI)/2)
DELAY = Sequence([d;d])
seq = DELAY + EPI
println("EPI successfully loaded! (TE = $(TE*1e3) ms)")
scanner = []
signal = 0
kdata = [0.0im 0.; 0. 0.]
MRIsim.print_gpus_info()

Δt = 4e-6 #<- simulate param
t = collect(Δt:Δt:MRIsim.dur(seq))
Nphant, Nt = prod(size(phantom)), length(t)
N_parts = floor(Int, Nphant*Nt/2.7e6*1/3) #heuristic
println("Dividing simulation in Nblocks=$N_parts")
S = @time MRIsim.run_sim2D_times_iter(phantom,seq,t;N_parts) #run_sim2D_times_iter run_sim2D_spin
signal = S[MRIsim.get_DAC_on(seq,t)]/prod(size(phantom)) #Acquired data
S = nothing
Nx = Ny = 100 #hardcoded by now
kdata = reshape(signal,(Nx,Ny)) #Turning into kspace image
kdata[:,2:2:Ny,:] = kdata[Nx:-1:1,2:2:Ny] #Flip in freq-dir for the EPI
kdata = convert(Array{Complex{Float64},2},kdata);
##
using CUDA
a = cu([1,2])
b = a .* a

@info "GPU related operations do not work yet, wait until Julia 1.6"

##
using PackageCompiler
PackageCompiler.create_sysimage(:MRIsim;
sysimage_path="C:/Users/ccp/Documents/MRIsim.jl/src/examples/sysimageSpinLab.dll",
precompile_execution_file="C:/Users/ccp/Documents/MRIsim.jl/src/examples/precompileSpinLab.jl")
# julia --sysimage C:/Users/ccp/Documents/MRIsim.jl/src/examples/sysimageSpinLab.dll --project=@.

##
# using PackageCompiler
# PackageCompiler.create_sysimage([:Plots, :JuMP, :Ipopt, :Blink, :PlotlyJS], replace_default=true)


text:
function getMinAndMaxAss(borrowingAllowed, minIncome, maxIncome, startA)

# A fuction that givs the minimum and maximum on the asset grid in each
# year. The minimum is the natural borrowing constraint. The maximum is how
# much one would have if one saved everything

# ------------------------------------------------------------------------
# Initialise the output matrices
   BC   = Array(Float64,T+1);
   maxA = Array(Float64,T+1);
# ------------------------------------------------------------------------
# Iteratively, calculate the borrowing constraints, and maximum on asset
# grid

 # Borrowing constraints
 BC[T + 1] = 0;
 for ixt = T:-1:1
     BC[ixt] = BC[ixt+1]/(1+r) - minIncome[ixt, 1] + minCons;
 end

# if borrowing is not allowed, replace negative points in the borrowing
# constraint with zero
if (borrowingAllowed == 0 )
	for i=1:length(BC)
    	if BC[i]<0;
    		BC[i] = 0
    	end
    end
end

# Maximum Assets
maxA[1] = startA;
for ixt = 2:1:T+1
    maxA[ixt] = (maxA[ixt - 1] + maxIncome[ixt-1] ) * (1+r);
end

#If starting assets are 0 we will have maxA<0 in period 1 and BC>0.
#Therefore BC>maxA (minA<maxA). Ensure that this is not the case by
#replacing any element of maxA that is less than the corresponding element
#of BC with that corresponding element plus 1
for ixt = 1:1:T+1
  if maxA[ixt] <= BC[ixt]
      maxA[ixt] = BC[ixt] + 1;
  end
end

return BC, maxA

end

max_stars_repo_name: floswald/ucl-econ-julia

max_issues_repo_name: floswald/ucl-econ-julia

max_forks_repo_name: floswald/ucl-econ-julia


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Tests related to the text backend.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

data = Any[1    false      1.0     0x01 ;
           2     true      2.0     0x02 ;
           3    false      3.0     0x03 ;
           4     true      4.0     0x04 ;
           5    false      5.0     0x05 ;
           6     true      6.0     0x06 ;]

# Default
# ==============================================================================

@testset "Default" begin

    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│      1 │  false │    1.0 │      1 │
│      2 │   true │    2.0 │      2 │
│      3 │  false │    3.0 │      3 │
│      4 │   true │    4.0 │      4 │
│      5 │  false │    5.0 │      5 │
│      6 │   true │    6.0 │      6 │
└────────┴────────┴────────┴────────┘
"""
    result = pretty_table(String, data)
    @test result == expected
end

# Alignments
# ==============================================================================

@testset "Alignments" begin
    # Left
    # ==========================================================================
    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│ 1      │ false  │ 1.0    │ 1      │
│ 2      │ true   │ 2.0    │ 2      │
│ 3      │ false  │ 3.0    │ 3      │
│ 4      │ true   │ 4.0    │ 4      │
│ 5      │ false  │ 5.0    │ 5      │
│ 6      │ true   │ 6.0    │ 6      │
└────────┴────────┴────────┴────────┘
"""
    result = pretty_table(String, data; alignment = :l)
    @test result == expected

    # Center
    # ==========================================================================
    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│   1    │ false  │  1.0   │   1    │
│   2    │  true  │  2.0   │   2    │
│   3    │ false  │  3.0   │   3    │
│   4    │  true  │  4.0   │   4    │
│   5    │ false  │  5.0   │   5    │
│   6    │  true  │  6.0   │   6    │
└────────┴────────┴────────┴────────┘
"""
    result = pretty_table(String, data; alignment = :c)
    @test result == expected

    # Per column configuration
    # ==========================================================================

    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│ 1      │  false │  1.0   │      1 │
│ 2      │   true │  2.0   │      2 │
│ 3      │  false │  3.0   │      3 │
│ 4      │   true │  4.0   │      4 │
│ 5      │  false │  5.0   │      5 │
│ 6      │   true │  6.0   │      6 │
└────────┴────────┴────────┴────────┘
"""
    result = pretty_table(String, data; alignment = [:l,:r,:c,:r])
    @test result == expected

    # Cell override
    # ==========================================================================

    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│ 1      │  false │  1.0   │ 1      │
│ 2      │   true │  2.0   │      2 │
│      3 │ false  │  3.0   │   3    │
│ 4      │   true │  4.0   │   4    │
│ 5      │  false │  5.0   │      5 │
│ 6      │   true │  6.0   │ 6      │
└────────┴────────┴────────┴────────┘
"""
    result = pretty_table(String, data;
                          alignment = [:l,:r,:c,:r],
                          cell_alignment = Dict( (3,1) => :r,
                                                 (3,2) => :l,
                                                 (1,4) => :l,
                                                 (3,4) => :c,
                                                 (4,4) => :c,
                                                 (6,4) => :l ))
    @test result == expected
end

# Filters
# ==============================================================================

@testset "Filters" begin
    expected = """
┌─────┬────────┬────────┐
│ Row │ Col. 1 │ Col. 3 │
├─────┼────────┼────────┤
│   2 │      2 │    2.0 │
│   4 │      4 │    4.0 │
│   6 │      6 │    6.0 │
└─────┴────────┴────────┘
"""

    result = pretty_table(String, data;
                          filters_row     = ( (data,i) -> i%2 == 0,),
                          filters_col     = ( (data,i) -> i%2 == 1,),
                          formatter       = ft_printf("%.3",3),
                          show_row_number = true)
    @test result == expected

    expected = """
┌─────┬────────┬────────┐
│ Row │ Col. 1 │ Col. 3 │
├─────┼────────┼────────┤
│   2 │   2    │ 2.0    │
│   4 │   4    │ 4.0    │
│   6 │   6    │ 6.0    │
└─────┴────────┴────────┘
"""

    result = pretty_table(String, data;
                          filters_row     = ( (data,i) -> i%2 == 0,),
                          filters_col     = ( (data,i) -> i%2 == 1,),
                          formatter       = ft_printf("%.3",3),
                          show_row_number = true,
                          alignment       = [:c,:l,:l,:c])
    @test result == expected
end

# Formatter
# ==============================================================================

@testset "Formatter" begin
    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│      1 │  false │      1 │      1 │
│      0 │   true │      0 │      0 │
│      3 │  false │      3 │      3 │
│      0 │   true │      0 │      0 │
│      5 │  false │      5 │      5 │
│      0 │   true │      0 │      0 │
└────────┴────────┴────────┴────────┘
"""
    formatter = Dict(0 => (v,i) -> isodd(i) ? i : 0,
                     2 => (v,i) -> v)
    result = pretty_table(String, data; formatter = formatter)
    @test result == expected
end

# Show row number
# ==============================================================================

@testset "Show row number" begin
    expected = """
┌─────┬────────┬────────┬────────┬────────┐
│ Row │ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├─────┼────────┼────────┼────────┼────────┤
│   1 │ 1      │  false │  1.0   │      1 │
│   2 │ 2      │   true │  2.0   │      2 │
│   3 │ 3      │  false │  3.0   │      3 │
│   4 │ 4      │   true │  4.0   │      4 │
│   5 │ 5      │  false │  5.0   │      5 │
│   6 │ 6      │   true │  6.0   │      6 │
└─────┴────────┴────────┴────────┴────────┘
"""
    result = pretty_table(String, data;
                          alignment       = [:l,:r,:c,:r],
                          show_row_number = true)
    @test result == expected
end

# Pre-defined formats
# ==============================================================================

@testset "Pre-defined formats" begin

    # ascii_dots
    # ==========================================================================
    expected = """
.....................................
: Col. 1 : Col. 2 : Col. 3 : Col. 4 :
:........:........:........:........:
:      1 :  false :    1.0 :      1 :
:      2 :   true :    2.0 :      2 :
:      3 :  false :    3.0 :      3 :
:      4 :   true :    4.0 :      4 :
:      5 :  false :    5.0 :      5 :
:      6 :   true :    6.0 :      6 :
:........:........:........:........:
"""
    result = pretty_table(String, data, tf = ascii_dots)
    @test result == expected

    # ascii_rounded
    # ==========================================================================
    expected = """
.--------.--------.--------.--------.
| Col. 1 | Col. 2 | Col. 3 | Col. 4 |
:--------+--------+--------+--------:
|      1 |  false |    1.0 |      1 |
|      2 |   true |    2.0 |      2 |
|      3 |  false |    3.0 |      3 |
|      4 |   true |    4.0 |      4 |
|      5 |  false |    5.0 |      5 |
|      6 |   true |    6.0 |      6 |
'--------'--------'--------'--------'
"""
    result = pretty_table(String, data, tf = ascii_rounded)
    @test result == expected

    # borderless
    # ==========================================================================
    expected = """
  Col. 1   Col. 2   Col. 3   Col. 4

       1    false      1.0        1
       2     true      2.0        2
       3    false      3.0        3
       4     true      4.0        4
       5    false      5.0        5
       6     true      6.0        6
"""
    result = pretty_table(String, data, tf = borderless)
    @test result == expected

    # compact
    # ==========================================================================
    expected = """
 -------- -------- -------- --------
  Col. 1   Col. 2   Col. 3   Col. 4
 -------- -------- -------- --------
       1    false      1.0        1
       2     true      2.0        2
       3    false      3.0        3
       4     true      4.0        4
       5    false      5.0        5
       6     true      6.0        6
 -------- -------- -------- --------
"""
    result = pretty_table(String, data, tf = PrettyTables.compact)
    @test result == expected

    # markdown
    # ==========================================================================
    expected = """
| Col. 1 | Col. 2 | Col. 3 | Col. 4 |
|--------|--------|--------|--------|
|      1 |  false |    1.0 |      1 |
|      2 |   true |    2.0 |      2 |
|      3 |  false |    3.0 |      3 |
|      4 |   true |    4.0 |      4 |
|      5 |  false |    5.0 |      5 |
|      6 |   true |    6.0 |      6 |
"""
    result = pretty_table(String, data, tf = markdown)
    @test result == expected

    # matrix
    # ==========================================================================

    expected = """
┌                     ┐
│ 1   false   1.0   1 │
│ 2    true   2.0   2 │
│ 3   false   3.0   3 │
│ 4    true   4.0   4 │
│ 5   false   5.0   5 │
│ 6    true   6.0   6 │
└                     ┘
"""
    result = pretty_table(String, data, tf = matrix, noheader = true)
    @test result == expected

    # mysql
    # ==========================================================================
    expected = """
+--------+--------+--------+--------+
| Col. 1 | Col. 2 | Col. 3 | Col. 4 |
+--------+--------+--------+--------+
|      1 |  false |    1.0 |      1 |
|      2 |   true |    2.0 |      2 |
|      3 |  false |    3.0 |      3 |
|      4 |   true |    4.0 |      4 |
|      5 |  false |    5.0 |      5 |
|      6 |   true |    6.0 |      6 |
+--------+--------+--------+--------+
"""
    result = pretty_table(String, data, tf = mysql)
    @test result == expected

    # simple
    # ==========================================================================
    expected = """
========= ======== ======== =========
  Col. 1   Col. 2   Col. 3   Col. 4
========= ======== ======== =========
       1    false      1.0        1
       2     true      2.0        2
       3    false      3.0        3
       4     true      4.0        4
       5    false      5.0        5
       6     true      6.0        6
========= ======== ======== =========
"""
    result = pretty_table(String, data, tf = simple)
    @test result == expected

    # unicode_rounded
    # ==========================================================================
    expected = """
╭────────┬────────┬────────┬────────╮
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│      1 │  false │    1.0 │      1 │
│      2 │   true │    2.0 │      2 │
│      3 │  false │    3.0 │      3 │
│      4 │   true │    4.0 │      4 │
│      5 │  false │    5.0 │      5 │
│      6 │   true │    6.0 │      6 │
╰────────┴────────┴────────┴────────╯
"""
    result = pretty_table(String, data, tf = unicode_rounded)
    @test result == expected

    # Custom formats
    # ==========================================================================

    expected = """
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│      1 │  false │    1.0 │      1 │
│      2 │   true │    2.0 │      2 │
│      3 │  false │    3.0 │      3 │
│      4 │   true │    4.0 │      4 │
│      5 │  false │    5.0 │      5 │
│      6 │   true │    6.0 │      6 │
"""

    tf = PrettyTables.TextFormat(unicode, top_line = false, bottom_line = false)
    result = pretty_table(String, data, tf = tf)
    @test result == expected

    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
│      1 │  false │    1.0 │      1 │
│      2 │   true │    2.0 │      2 │
│      3 │  false │    3.0 │      3 │
│      4 │   true │    4.0 │      4 │
│      5 │  false │    5.0 │      5 │
│      6 │   true │    6.0 │      6 │
└────────┴────────┴────────┴────────┘
"""

    tf = PrettyTables.TextFormat(unicode, header_line = false)
    result = pretty_table(String, data, tf = tf)
    @test result == expected
end

# Pre-defined formatters
# ==============================================================================

@testset "Pre-defined formatters" begin

    # ft_round
    # ==========================================================================
    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│    1.0 │    0.0 │    1.0 │    1.0 │
│    2.0 │    1.0 │    2.0 │    2.0 │
│    3.0 │    0.0 │    3.0 │    3.0 │
│    4.0 │    1.0 │    4.0 │    4.0 │
│    5.0 │    0.0 │    5.0 │    5.0 │
│    6.0 │    1.0 │    6.0 │    6.0 │
└────────┴────────┴────────┴────────┘
"""
    result = pretty_table(String, data; formatter = ft_round(1))
    @test result == expected

    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│    1.0 │  false │    1.0 │      1 │
│    2.0 │   true │    2.0 │      2 │
│    3.0 │  false │    3.0 │      3 │
│    4.0 │   true │    4.0 │      4 │
│    5.0 │  false │    5.0 │      5 │
│    6.0 │   true │    6.0 │      6 │
└────────┴────────┴────────┴────────┘
"""
    result = pretty_table(String, data; formatter = ft_round(1,[3,1]))
    @test result == expected

    # ft_printf
    # ==========================================================================

    expected = """
┌──────────┬──────────┬──────────┬──────────┐
│   Col. 1 │   Col. 2 │   Col. 3 │   Col. 4 │
├──────────┼──────────┼──────────┼──────────┤
│    1.000 │    0.000 │    1.000 │    1.000 │
│    2.000 │    1.000 │    2.000 │    2.000 │
│    3.000 │    0.000 │    3.000 │    3.000 │
│    4.000 │    1.000 │    4.000 │    4.000 │
│    5.000 │    0.000 │    5.000 │    5.000 │
│    6.000 │    1.000 │    6.000 │    6.000 │
└──────────┴──────────┴──────────┴──────────┘
"""
    result = pretty_table(String,data; formatter = ft_printf("%8.3f"))
    @test result == expected

    expected = """
┌──────────┬────────┬────────┬──────────┐
│   Col. 1 │ Col. 2 │ Col. 3 │   Col. 4 │
├──────────┼────────┼────────┼──────────┤
│    1.000 │  false │    1.0 │    1.000 │
│    2.000 │   true │    2.0 │    2.000 │
│    3.000 │  false │    3.0 │    3.000 │
│    4.000 │   true │    4.0 │    4.000 │
│    5.000 │  false │    5.0 │    5.000 │
│    6.000 │   true │    6.0 │    6.000 │
└──────────┴────────┴────────┴──────────┘
"""
    result = pretty_table(String,data; formatter = ft_printf("%8.3f",[1,4]))
    @test result == expected

    expected = """
┌──────────┬────────┬────────┬──────────┐
│   Col. 1 │ Col. 2 │ Col. 3 │   Col. 4 │
├──────────┼────────┼────────┼──────────┤
│     1.00 │  false │    1.0 │   1.0000 │
│     2.00 │   true │    2.0 │   2.0000 │
│     3.00 │  false │    3.0 │   3.0000 │
│     4.00 │   true │    4.0 │   4.0000 │
│     5.00 │  false │    5.0 │   5.0000 │
│     6.00 │   true │    6.0 │   6.0000 │
└──────────┴────────┴────────┴──────────┘
"""
    result = pretty_table(String, data;
                          formatter = ft_printf(["%8.2f","%8.4f"],[1,4]))
    @test result == expected
end

# Sub-headers
# ==============================================================================

@testset "Sub-headers" begin
    expected = """
┌───┬───────┬─────┬───┐
│ 1 │     2 │   3 │ 4 │
├───┼───────┼─────┼───┤
│ 1 │ false │ 1.0 │ 1 │
│ 2 │  true │ 2.0 │ 2 │
│ 3 │ false │ 3.0 │ 3 │
│ 4 │  true │ 4.0 │ 4 │
│ 5 │ false │ 5.0 │ 5 │
│ 6 │  true │ 6.0 │ 6 │
└───┴───────┴─────┴───┘
"""
    result = pretty_table(String, data, [1;2;3;4])
    @test result == expected

    result = pretty_table(String, data, [1 2 3 4])
    @test result == expected

    expected = """
┌───┬───────┬─────┬───┐
│ 1 │     2 │   3 │ 4 │
│ 5 │     6 │   7 │ 8 │
├───┼───────┼─────┼───┤
│ 1 │ false │ 1.0 │ 1 │
│ 2 │  true │ 2.0 │ 2 │
│ 3 │ false │ 3.0 │ 3 │
│ 4 │  true │ 4.0 │ 4 │
│ 5 │ false │ 5.0 │ 5 │
│ 6 │  true │ 6.0 │ 6 │
└───┴───────┴─────┴───┘
"""
    result = pretty_table(String, data, [1 2 3 4; 5 6 7 8])
    @test result == expected
end

# Horizontal lines
# ==============================================================================

@testset "Horizontal lines" begin
    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│      1 │  false │    1.0 │      1 │
│      2 │   true │    2.0 │      2 │
├────────┼────────┼────────┼────────┤
│      3 │  false │    3.0 │      3 │
│      4 │   true │    4.0 │      4 │
├────────┼────────┼────────┼────────┤
│      5 │  false │    5.0 │      5 │
│      6 │   true │    6.0 │      6 │
└────────┴────────┴────────┴────────┘
"""
    result = pretty_table(String, data; hlines = findall(x->x == true, data[:,2]))
    @test result == expected

    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│      1 │  false │    1.0 │      1 │
│      2 │   true │    2.0 │      2 │
├........+........+........+........┤
│      3 │  false │    3.0 │      3 │
│      4 │   true │    4.0 │      4 │
├........+........+........+........┤
│      5 │  false │    5.0 │      5 │
│      6 │   true │    6.0 │      6 │
└────────┴────────┴────────┴────────┘
"""
    result = pretty_table(String, data;
                          hlines = findall(x->x == true, data[:,2]),
                          hlines_format = ('├','+','┤','.'))
    @test result == expected
end

# Line breaks inside cells
# ==============================================================================

@testset "Line breaks inside cells" begin
    data = ["This line contains\nthe velocity [m/s]" 10.0;
            "This line contains\nthe acceleration [m/s^2]" 1.0;
            "This line contains\nthe time from the\nbeginning of the simulation" 10;]

    header = ["Information", "Value"]

    expected = """
┌─────────────────────────────┬───────┐
│                 Information │ Value │
├─────────────────────────────┼───────┤
│          This line contains │  10.0 │
│          the velocity [m/s] │       │
│          This line contains │   1.0 │
│    the acceleration [m/s^2] │       │
│          This line contains │    10 │
│           the time from the │       │
│ beginning of the simulation │       │
└─────────────────────────────┴───────┘
"""

    result = pretty_table(String, data, header; linebreaks = true)
    @test result == expected

    expected = """
┌─────────────────────────────┬───────┐
│         Information         │ Value │
├─────────────────────────────┼───────┤
│     This line contains      │ 10.0  │
│     the velocity [m/s]      │       │
│     This line contains      │  1.0  │
│  the acceleration [m/s^2]   │       │
│     This line contains      │  10   │
│      the time from the      │       │
│ beginning of the simulation │       │
└─────────────────────────────┴───────┘
"""

    result = pretty_table(String, data, header; alignment = :c, linebreaks = true)
    @test result == expected

    expected = """
┌─────────────────────────────┬───────┐
│ Information                 │ Value │
├─────────────────────────────┼───────┤
│ This line contains          │ 10.0  │
│ the velocity [m/s]          │       │
│ This line contains          │ 1.0   │
│ the acceleration [m/s^2]    │       │
│ This line contains          │ 10    │
│ the time from the           │       │
│ beginning of the simulation │       │
└─────────────────────────────┴───────┘
"""

    result = pretty_table(String, data, header; alignment = :l, linebreaks = true)
    @test result == expected

    expected = """
┌────────────────────────────────────────────────────────────────────┬───────┐
│                                                        Information │ Value │
├────────────────────────────────────────────────────────────────────┼───────┤
│                             This line contains\\nthe velocity [m/s] │  10.0 │
│                       This line contains\\nthe acceleration [m/s^2] │   1.0 │
│ This line contains\\nthe time from the\\nbeginning of the simulation │    10 │
└────────────────────────────────────────────────────────────────────┴───────┘
"""

    result = pretty_table(String, data, header)
    @test result == expected
end

# Hiding header and sub-header
# ==============================================================================

@testset "Hide header and sub-header" begin

    # Header
    # --------------------------------------------------------------------------

    expected = """
┌───┬───────┬─────┬───┐
│ 1 │ false │ 1.0 │ 1 │
│ 2 │  true │ 2.0 │ 2 │
│ 3 │ false │ 3.0 │ 3 │
│ 4 │  true │ 4.0 │ 4 │
│ 5 │ false │ 5.0 │ 5 │
│ 6 │  true │ 6.0 │ 6 │
└───┴───────┴─────┴───┘
"""

    result = pretty_table(String, data; noheader = true)
    @test result == expected

    result = pretty_table(String, data, [1 2]; noheader = true)
    @test result == expected

    expected = """
┌───┬───────┬─────┬───┐
│ 1 │ false │ 1.0 │ 1 │
│ 2 │ true  │ 2.0 │ 2 │
│ 3 │ false │ 3.0 │ 3 │
│ 4 │ true  │ 4.0 │ 4 │
│ 5 │ false │ 5.0 │ 5 │
│ 6 │ true  │ 6.0 │ 6 │
└───┴───────┴─────┴───┘
"""
    result = pretty_table(String, data, [1 2]; alignment = :l, noheader = true)
    @test result == expected

    # Sub-header
    # --------------------------------------------------------------------------

    expected = """
┌───┬───────┬─────┬───┐
│ 1 │     2 │   3 │ 4 │
├───┼───────┼─────┼───┤
│ 1 │ false │ 1.0 │ 1 │
│ 2 │  true │ 2.0 │ 2 │
│ 3 │ false │ 3.0 │ 3 │
│ 4 │  true │ 4.0 │ 4 │
│ 5 │ false │ 5.0 │ 5 │
│ 6 │  true │ 6.0 │ 6 │
└───┴───────┴─────┴───┘
"""

    header = [1 2 3 4; "this is" "a very very" "big" "sub-header"]
    result = pretty_table(String, data, header; nosubheader = true)
    @test result == expected
end

# Print vectors
# ==============================================================================

@testset "Print vectors" begin

    vec = 0:1:10

    expected = """
┌────────┐
│ Col. 1 │
├────────┤
│      0 │
│      1 │
│      2 │
│      3 │
│      4 │
│      5 │
│      6 │
│      7 │
│      8 │
│      9 │
│     10 │
└────────┘
"""

    result = pretty_table(String, vec)
    @test result == expected

    expected = """
┌─────┬────────┐
│ Row │ Col. 1 │
├─────┼────────┤
│   1 │   0    │
│   2 │   1    │
│   3 │   2    │
│   4 │   3    │
│   5 │   4    │
│   6 │   5    │
│   7 │   6    │
│   8 │   7    │
│   9 │   8    │
│  10 │   9    │
│  11 │   10   │
└─────┴────────┘
"""

    result = pretty_table(String, vec; alignment = :c, show_row_number = true)
    @test result == expected

    expected = """
┌────────────┐
│     Header │
│ Sub-header │
├────────────┤
│          0 │
│          1 │
│          2 │
│          3 │
│          4 │
│          5 │
│          6 │
│          7 │
│          8 │
│          9 │
│         10 │
└────────────┘
"""

    result = pretty_table(String, vec, ["Header", "Sub-header"])
    @test result == expected

    @test_throws Exception pretty_table(vec, ["1" "1"])
end

# Cropping
# ==============================================================================

@testset "Cropping" begin
    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│      1 │  false │    1.0 │      1 │
│      2 │   true │    2.0 │      2 │
│      3 │  false │    3.0 │      3 │
│      4 │   true │    4.0 │      4 │
│      5 │  false │    5.0 │      5 │
│      6 │   true │    6.0 │      6 │
└────────┴────────┴────────┴────────┘
"""

    result = pretty_table(String, data, screen_size = (10,10), crop = :none)
    @test result == expected

    result = pretty_table(String, data, screen_size = (-1,-1), crop = :both)
    @test result == expected

    expected = """
┌────────┬────────┬────────┬ ⋯
│ Col. 1 │ Col. 2 │ Col. 3 │ ⋯
├────────┼────────┼────────┼ ⋯
│      1 │  false │    1.0 │ ⋯
│      2 │   true │    2.0 │ ⋯
│   ⋮    │   ⋮    │   ⋮    │ ⋯
└────────┴────────┴────────┴ ⋯
"""

    result = pretty_table(String, data, screen_size = (10,30), crop = :both)
    @test result == expected

    result = pretty_table(String, data, screen_size = (10,30))
    @test result == expected

    expected = """
┌────────┬────────┬────────┬ ⋯
│ Col. 1 │ Col. 2 │ Col. 3 │ ⋯
├────────┼────────┼────────┼ ⋯
│      1 │  false │    1.0 │ ⋯
│      2 │   true │    2.0 │ ⋯
│      3 │  false │    3.0 │ ⋯
│      4 │   true │    4.0 │ ⋯
│      5 │  false │    5.0 │ ⋯
│      6 │   true │    6.0 │ ⋯
└────────┴────────┴────────┴ ⋯
"""

    result = pretty_table(String, data, screen_size = (10,30), crop = :horizontal)
    @test result == expected

    result = pretty_table(String, data, screen_size = (-1,30), crop = :both)
    @test result == expected

    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│      1 │  false │    1.0 │      1 │
│      2 │   true │    2.0 │      2 │
│   ⋮    │   ⋮    │   ⋮    │   ⋮    │
└────────┴────────┴────────┴────────┘
"""

    result = pretty_table(String, data, screen_size = (10,30), crop = :vertical)
    @test result == expected

    result = pretty_table(String, data, screen_size = (10,-1), crop = :both)
    @test result == expected

    model = [1.123456789 for i = 1:8, j = 1:10]

    expected = """
┌───┬─────────────┬────┬─────┬──────┬───────┬────────┬─────────┬──────────┬───────────┐
│ … │ Col. 2      │ C… │ Co… │ Col… │ Col.… │ Col. 7 │ Col. 8  │  Col. 9  │   Col. 10 │
├───┼─────────────┼────┼─────┼──────┼───────┼────────┼─────────┼──────────┼───────────┤
│ … │ 1.123456789 │ 1… │ 1.… │ 1.1… │ 1.12… │ 1.123… │ 1.1234… │ 1.12345… │ 1.123456… │
│ … │ 1.123456789 │ 1… │ 1.… │ 1.1… │ 1.12… │ 1.123… │ 1.1234… │ 1.12345… │ 1.123456… │
│ … │ 1.123456789 │ 1… │ 1.… │ 1.1… │ 1.12… │ 1.123… │ 1.1234… │ 1.12345… │ 1.123456… │
│ … │ 1.123456789 │ 1… │ 1.… │ 1.1… │ 1.12… │ 1.123… │ 1.1234… │ 1.12345… │ 1.123456… │
│ … │ 1.123456789 │ 1… │ 1.… │ 1.1… │ 1.12… │ 1.123… │ 1.1234… │ 1.12345… │ 1.123456… │
│ … │ 1.123456789 │ 1… │ 1.… │ 1.1… │ 1.12… │ 1.123… │ 1.1234… │ 1.12345… │ 1.123456… │
│ … │ 1.123456789 │ 1… │ 1.… │ 1.1… │ 1.12… │ 1.123… │ 1.1234… │ 1.12345… │ 1.123456… │
│ … │ 1.123456789 │ 1… │ 1.… │ 1.1… │ 1.12… │ 1.123… │ 1.1234… │ 1.12345… │ 1.123456… │
└───┴─────────────┴────┴─────┴──────┴───────┴────────┴─────────┴──────────┴───────────┘
"""

    result = pretty_table(String, model,
                          columns_width = [1, 0, 2, 3, 4, 5, 6, 7, 8, 9],
                          alignment = [:l, :l, :r, :c, :r, :c, :l, :l, :c, :r])

    expected = """
┌────────┬────────┬────────┬────────┬────────┬────────┬────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │ Col. 5 │ Col. 6 │ Col. 7 │ Col. 8 │ Col. 9 │ Col. … │
├────────┼────────┼────────┼────────┼────────┼────────┼────────┼────────┼────────┼────────┤
│ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │
│ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │
│ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │
│ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │
│ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │
│ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │
│ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │
│ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │ 1.123… │
└────────┴────────┴────────┴────────┴────────┴────────┴────────┴────────┴────────┴────────┘
"""

    result = pretty_table(String, model,
                          columns_width = 6,
                          alignment = [:l, :l, :r, :c, :r, :c, :l, :l, :c, :r])

    @test result == expected
end

# Auto wrapping
# ==============================================================================

@testset "Auto wrapping" begin
    table = [1 """Ouviram do Ipiranga as margens plácidas
                  De um povo heróico o brado retumbante,
                  E o sol da Liberdade, em raios fúlgidos,
                  Brilhou no céu da Pátria nesse instante.""";
             2 """Se o penhor dessa igualdade
                  Conseguimos conquistar com braço forte,
                  Em teu seio, ó Liberdade,
                  Desafia o nosso peito a própria morte!""";
             3 """Ó Pátria amada, Idolatrada, Salve! Salve!
                  Brasil, um sonho intenso, um raio vívido
                  De amor e de esperança à terra desce,
                  Se em teu formoso céu, risonho e límpido,"""]

    header = ["Verse number", "Verse"]

    expected = """
┌──────────────┬────────────────────────────────┐
│ Verse number │                          Verse │
├──────────────┼────────────────────────────────┤
│            1 │         Ouviram do Ipiranga as │
│              │               margens plácidas │
│              │     De um povo heróico o brado │
│              │                    retumbante, │
│              │       E o sol da Liberdade, em │
│              │                raios fúlgidos, │
│              │       Brilhou no céu da Pátria │
│              │                nesse instante. │
├──────────────┼────────────────────────────────┤
│            2 │    Se o penhor dessa igualdade │
│              │     Conseguimos conquistar com │
│              │                   braço forte, │
│              │      Em teu seio, ó Liberdade, │
│              │        Desafia o nosso peito a │
│              │                 própria morte! │
├──────────────┼────────────────────────────────┤
│            3 │    Ó Pátria amada, Idolatrada, │
│              │                  Salve! Salve! │
│              │   Brasil, um sonho intenso, um │
│              │                    raio vívido │
│              │       De amor e de esperança à │
│              │                   terra desce, │
│              │         Se em teu formoso céu, │
│              │             risonho e límpido, │
└──────────────┴────────────────────────────────┘
"""

    result = pretty_table(String, table, header,
                          autowrap      = true,
                          linebreaks    = true,
                          hlines        = 1:4,
                          columns_width = [-1,30])

    @test result == expected

    expected = """
┌──────────────┬────────────────────────────────┐
│ Verse number │             Verse              │
├──────────────┼────────────────────────────────┤
│      1       │     Ouviram do Ipiranga as     │
│              │        margens plácidas        │
│              │   De um povo heróico o brado   │
│              │          retumbante,           │
│              │    E o sol da Liberdade, em    │
│              │        raios fúlgidos,         │
│              │    Brilhou no céu da Pátria    │
│              │        nesse instante.         │
├──────────────┼────────────────────────────────┤
│      2       │  Se o penhor dessa igualdade   │
│              │   Conseguimos conquistar com   │
│              │          braço forte,          │
│              │   Em teu seio, ó Liberdade,    │
│              │    Desafia o nosso peito a     │
│              │         própria morte!         │
├──────────────┼────────────────────────────────┤
│      3       │  Ó Pátria amada, Idolatrada,   │
│              │         Salve! Salve!          │
│              │  Brasil, um sonho intenso, um  │
│              │          raio vívido           │
│              │    De amor e de esperança à    │
│              │          terra desce,          │
│              │     Se em teu formoso céu,     │
│              │       risonho e límpido,       │
└──────────────┴────────────────────────────────┘
"""

    result = pretty_table(String, table, header,
                          alignment     = :c,
                          autowrap      = true,
                          linebreaks    = true,
                          hlines        = 1:4,
                          columns_width = [-1,30])

    @test result == expected

    expected = """
┌──────────────┬────────────────────────────────┐
│ Verse number │ Verse                          │
├──────────────┼────────────────────────────────┤
│ 1            │ Ouviram do Ipiranga as         │
│              │ margens plácidas               │
│              │ De um povo heróico o brado     │
│              │ retumbante,                    │
│              │ E o sol da Liberdade, em       │
│              │ raios fúlgidos,                │
│              │ Brilhou no céu da Pátria       │
│              │ nesse instante.                │
├──────────────┼────────────────────────────────┤
│ 2            │ Se o penhor dessa igualdade    │
│              │ Conseguimos conquistar com     │
│              │ braço forte,                   │
│              │ Em teu seio, ó Liberdade,      │
│              │ Desafia o nosso peito a        │
│              │ própria morte!                 │
├──────────────┼────────────────────────────────┤
│ 3            │ Ó Pátria amada, Idolatrada,    │
│              │ Salve! Salve!                  │
│              │ Brasil, um sonho intenso, um   │
│              │ raio vívido                    │
│              │ De amor e de esperança à       │
│              │ terra desce,                   │
│              │ Se em teu formoso céu,         │
│              │ risonho e límpido,             │
└──────────────┴────────────────────────────────┘
"""

    result = pretty_table(String, table, header,
                          alignment     = :l,
                          autowrap      = true,
                          linebreaks    = true,
                          hlines        = 1:4,
                          columns_width = [-1,30])

    @test result == expected
end

# Dictionaries
# ==============================================================================

@testset "Dictionaries" begin
    dict = Dict{Int64,String}(1 => "Jan", 2 => "Feb", 3 => "Mar", 4 => "Apr",
                              5 => "May", 6 => "Jun")

    expected = """
┌───────┬────────┐
│  Keys │ Values │
│ Int64 │ String │
├───────┼────────┤
│     1 │    Jan │
│     2 │    Feb │
│     3 │    Mar │
│     4 │    Apr │
│     5 │    May │
│     6 │    Jun │
└───────┴────────┘
"""

    result = pretty_table(String, dict, sortkeys = true)
    @test result == expected
end

# Test if we can print `missing` and `nothing`
# ==============================================================================

@testset "Print missing and nothing" begin
    matrix = [missing missing; nothing nothing; missing nothing]

    expected = """
┌─────────┬─────────┐
│  Col. 1 │  Col. 2 │
├─────────┼─────────┤
│ missing │ missing │
│ nothing │ nothing │
│ missing │ nothing │
└─────────┴─────────┘
"""

    result = pretty_table(String, matrix)
    @test result == expected
end

# Helpers
# ==============================================================================

@testset "@pt" begin
    # To get the output of the macro @pt, we must redirect the stdout.
    old_stdout = stdout
    in, out    = redirect_stdout()

    # The configurations must not interfere with the printings of any type
    # except for the `Dict`.
    @ptconf sortkeys = true

    # Test 1
    # --------------------------------------------------------------------------

    expected = """
.......................
: 1 :     2 :   3 : 4 :
:...:.......:.....:...:
: 1 : false : 1.0 : 1 :
: 2 :  true : 2.0 : 2 :
: 3 : false : 3.0 : 3 :
: 4 :  true : 4.0 : 4 :
: 5 : false : 5.0 : 5 :
: 6 :  true : 6.0 : 6 :
:...:.......:.....:...:
"""

    @ptconf tf = ascii_dots
    @pt :header = ["1","2","3","4"] data

    result = String(readavailable(in))
    @test result == expected

    # Test 2
    # --------------------------------------------------------------------------

    expected = """
╭────────┬────────┬────────┬────────╮
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│   1    │ false  │  1.0   │   1    │
│   2    │  true  │  2.0   │   2    │
│   3    │ false  │  3.0   │   3    │
│   4    │  true  │  4.0   │   4    │
│   5    │ false  │  5.0   │   5    │
│   6    │  true  │  6.0   │   6    │
╰────────┴────────┴────────┴────────╯
"""

    @ptconf tf = unicode_rounded
    @ptconf alignment = :c
    @pt data

    result = String(readavailable(in))
    @test result == expected

    # Test 3
    # --------------------------------------------------------------------------

    expected = """
.----------.----------.
| Column 1 | Column 2 |
|   Sub. 1 |   Sub. 2 |
:----------+----------:
|        1 |        2 |
|        3 |        4 |
|        5 |        6 |
'----------'----------'
"""

    header = ["Column 1" "Column 2"; "Sub. 1" "Sub. 2"]

    @ptconfclean
    @ptconf tf = ascii_rounded
    @pt :header = header [1 2; 3 4; 5 6]

    result = String(readavailable(in))
    @test result == expected

    # Test 4
    # --------------------------------------------------------------------------

    expected = """
.-------.---------------------.
| Keys  | Values              |
| Int64 | String              |
:-------+---------------------:
| 1     | São José dos Campos |
| 2     | SP                  |
| 3     | Brasil              |
'-------'---------------------'
"""

    @ptconf sortkeys = true
    @ptconf alignment = :l
    @pt d = Dict(Int64(1) => "São José dos Campos",
                 Int64(2) => "SP",
                 Int64(3) => "Brasil")

    result = String(readavailable(in))
    @test result == expected

    # Restore the original stdout.
    close(in)
    close(out)
    redirect_stdout(old_stdout)

    # The expression after `@pt` must be evaluated. Hence, `d` must be a dict
    # with the defined elements.
    @test d isa Dict
    @test d[1] == "São José dos Campos"
    @test d[2] == "SP"
    @test d[3] == "Brasil"
end

# Issue #4
# ==============================================================================

@testset "Issue #4 - Header with escape sequences" begin

    expected = """
┌───┬───────┬─────┬───┐
│ 1 │   2\\n │   3 │ 4 │
├───┼───────┼─────┼───┤
│ 1 │ false │ 1.0 │ 1 │
│ 2 │  true │ 2.0 │ 2 │
│ 3 │ false │ 3.0 │ 3 │
│ 4 │  true │ 4.0 │ 4 │
│ 5 │ false │ 5.0 │ 5 │
│ 6 │  true │ 6.0 │ 6 │
└───┴───────┴─────┴───┘
"""

    result = pretty_table(String, data, ["1" "2\n" "3" "4"])
    @test result == expected
end

# Issue #9
# ==============================================================================

@testset "Issue #9 - Printing table with row number column and sub-headers" begin
    expected = """
┌─────┬───┬───────┬─────┬────┐
│ Row │ 1 │     2 │   3 │  4 │
│     │ 5 │     6 │   7 │  8 │
│     │ 9 │    10 │  11 │ 12 │
├─────┼───┼───────┼─────┼────┤
│   1 │ 1 │ false │ 1.0 │  1 │
│   2 │ 2 │  true │ 2.0 │  2 │
│   3 │ 3 │ false │ 3.0 │  3 │
│   4 │ 4 │  true │ 4.0 │  4 │
│   5 │ 5 │ false │ 5.0 │  5 │
│   6 │ 6 │  true │ 6.0 │  6 │
└─────┴───┴───────┴─────┴────┘
"""
    result = pretty_table(String, data, [1  2  3  4;
                                         5  6  7  8;
                                         9 10 11 12];
                          show_row_number = true)
    @test result == expected
end

# Issue #10
# ==============================================================================

@testset "Issue #10 - Size of the row number column without header" begin
    v = 0:1:6

    expected = """
┌───┬───┐
│ 1 │ 0 │
│ 2 │ 1 │
│ 3 │ 2 │
│ 4 │ 3 │
│ 5 │ 4 │
│ 6 │ 5 │
│ 7 │ 6 │
└───┴───┘
"""

    result = pretty_table(String, v; noheader = true, show_row_number = true)
    @test result == expected
end

# Issue #16
# ==============================================================================

@testset "Issue #16 - Printing DataFrames that contains strings" begin
    df = DataFrame(:a => Int64[1, 2], :b => ["A", "B"]);

    expected = """
┌───────┬────────┐
│     a │      b │
│ Int64 │ String │
├───────┼────────┤
│     1 │      A │
│     2 │      B │
└───────┴────────┘
"""

    result = pretty_table(String, df)
    @test result == expected
end

# Issue #19
# ==============================================================================

@testset "Issue #19 - ft_printf with cells that are not numbers" begin
    matrix = [1 1; 2 2; 3 3; "teste" "teste"; 4 4; 5 5; true true; :s :s]

    expected = """
┌────────────┬────────┐
│     Col. 1 │ Col. 2 │
├────────────┼────────┤
│       1.00 │      1 │
│       2.00 │      2 │
│       3.00 │      3 │
│      teste │  teste │
│       4.00 │      4 │
│       5.00 │      5 │
│       1.00 │   true │
│          s │      s │
└────────────┴────────┘
"""

    result = pretty_table(String, matrix; formatter = ft_printf("%10.2f",[1]))
    @test result == expected

end

# Issue #22
# ==============================================================================

@testset "Issue #22 - Strings with quotes" begin
    # Without linebreaks.
    matrix = [1 "teste\"teste"
              2 "teste\"\"teste"]

    expected = """
┌────────┬──────────────┐
│ Col. 1 │       Col. 2 │
├────────┼──────────────┤
│      1 │  teste"teste │
│      2 │ teste""teste │
└────────┴──────────────┘
"""

    result = pretty_table(String, matrix)
    @test result == expected

    # With linebreaks.
    matrix = [ 1 """
    function str(str = "one string")
        return str
    end"""
    2 """
    function str(str = "")
        if isempty(str)
            return "one string"
        else
            return str
        end
    end"""]

    expected = """
┌────────┬──────────────────────────────────┐
│ Col. 1 │ Col. 2                           │
├────────┼──────────────────────────────────┤
│ 1      │ function str(str = "one string") │
│        │     return str                   │
│        │ end                              │
│ 2      │ function str(str = "")           │
│        │     if isempty(str)              │
│        │         return "one string"      │
│        │     else                         │
│        │         return str               │
│        │     end                          │
│        │ end                              │
└────────┴──────────────────────────────────┘
"""

    result = pretty_table(String, matrix; alignment = :l, linebreaks = true)

    @test result == expected
end

# Issue #24
# ==============================================================================

@testset "Issue #24 - Tables compatibility" begin
    # Named tuple of vectors (satisfies Tables interface).
    ctable = (x = Int64[1, 2, 3, 4], y = ["a", "b", "c", "d"])

    cresult = sprint(pretty_table, ctable)

    expected = """
┌───────┬────────┐
│     x │      y │
│ Int64 │ String │
├───────┼────────┤
│     1 │      a │
│     2 │      b │
│     3 │      c │
│     4 │      d │
└───────┴────────┘
"""

  @test cresult == expected
end

# Issue #28
# ==============================================================================

@testset "Issue #28 - Tables.api must have priority when printing" begin
    # A DataFrame is compliant with Tables.jl API.
    df = DataFrame(x = Int64(1):Int64(3),
                   y = 'a':'c',
                   z = ["String 1";"String 2";"String 3"]);

    # Thus, the following 3 calls must provide the same results.
    result_1 = pretty_table(String, df)
    result_2 = pretty_table(String, Tables.rowtable(df))
    result_3 = pretty_table(String, Tables.columns(df))

    expected = """
┌───────┬──────┬──────────┐
│     x │    y │        z │
│ Int64 │ Char │   String │
├───────┼──────┼──────────┤
│     1 │    a │ String 1 │
│     2 │    b │ String 2 │
│     3 │    c │ String 3 │
└───────┴──────┴──────────┘
"""

    @test result_1 == result_2 == result_3 == expected
end


text:
module BracketingTests

import ..Algopt # ensure we are using the correct Algopt

using Test
using Algopt.Bracketing:
search, find_bracket, search_bracket,
Bracket,
GoldenSection

using Algopt.TestFunctions: Rosenbrock

quad2 = x->(x - 2)^2

ros = Rosenbrock(a = 1, b = 5)

@testset "Algopt.Bracketing.search(::GoldenSection)" begin
    @test 2 ≈ search(GoldenSection(), quad2)
    @test 2 ≈ search(GoldenSection(), quad2, -1e6)
    @test !(2 ≈ search(GoldenSection(max_iter = 2), quad2))

    @test abs(0 - search(GoldenSection(), x->x'x)) <= eps()

    # Line search of α on 2D Rosenbrock function
    ros_line = α->ros.f([10, 10] + α * [-1, -1])
    @test 9 ≈ search(GoldenSection(), ros_line, 0)
end

@testset "Algopt.Bracketing.find_bracket" begin
    bracket = find_bracket(quad2, 0)
    @test bracket.left < 2
    @test bracket.right > 2

    bracket = find_bracket(quad2, 10)
    @test bracket.left < 2
    @test bracket.right > 2

    bracket = find_bracket(quad2, -1e6)
    @test bracket.left < 2
    @test bracket.right > 2
end

@testset "Algopt.Bracketing.search_bracket" begin
    bracket = search_bracket(GoldenSection(), quad2, Bracket(0, 10))
    @test bracket.left ≈ 2
    @test bracket.right ≈ 2

    bracket = search_bracket(GoldenSection(), quad2, Bracket(-1e6, 1e6))
    @test bracket.left ≈ 2
    @test bracket.right ≈ 2

    bracket = search_bracket(GoldenSection(max_iter = 2), quad2, Bracket(-1e6, 1e6))
    @test !(bracket.left ≈ 2)
    @test !(bracket.right ≈ 2)
end

end # module


text:
@testset "ReificationDifferentFrom2Indicator: $(fct_type), type $(T)" for fct_type in ["vector of variables", "vector affine function"], T in [Int, Float64]
    base_model = if T == Int
        IntDifferentFromIndicatorMILPModel{Int}()
    elseif T == Float64
        FloatDifferentFromIndicatorMILPModel{Float64}()
    else
        @assert false
    end
    mock = MOIU.MockOptimizer(base_model)
    model = COIB.ReificationDifferentFrom2Indicator{T}(mock)

    if T == Int
        @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    end
    @test MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{T},
        MOI.EqualTo{T},
    )
    @test MOIB.supports_bridging_constraint(
        model,
        MOI.VectorAffineFunction{T},
        CP.Reification{CP.DifferentFrom{T}},
    )
    @test MOIB.supports_bridging_constraint(
        model,
        MOI.VectorOfVariables,
        CP.Reification{CP.DifferentFrom{T}},
    )

    x, _ = MOI.add_constrained_variable(model, MOI.ZeroOne())
    if T == Int
        y, _ = MOI.add_constrained_variable(model, MOI.Integer())
    elseif T == Float64
        y = MOI.add_variable(model)
    end

    fct = if fct_type == "vector of variables"
        MOI.VectorOfVariables([x, y])
    elseif fct_type == "vector affine function"
        MOIU.vectorize([x, y])
    else
        @assert false
    end
    c = MOI.add_constraint(model, fct, CP.Reification(CP.DifferentFrom(zero(T))))

    @test MOI.is_valid(model, x)
    @test MOI.is_valid(model, y)
    @test MOI.is_valid(model, c)

    bridge = MOIBC.bridges(model)[MOI.ConstraintIndex{MOI.VectorOfVariables, CP.Reification{CP.DifferentFrom{T}}}(-1)]

    @testset "Bridge properties" begin
        @test MOIBC.concrete_bridge_type(typeof(bridge), MOI.VectorOfVariables, CP.Reification{CP.DifferentFrom{T}}) == typeof(bridge)
        @test MOIB.added_constrained_variable_types(typeof(bridge)) == Tuple{Type}[]
        @test MOIB.added_constraint_types(typeof(bridge)) == [
            (MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.EqualTo{T}}),
            (MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ZERO, CP.DifferentFrom{T}}),
        ]

        @test MOI.get(bridge, MOI.NumberOfVariables()) == 0
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.EqualTo{T}}}()) == 1
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ZERO, CP.DifferentFrom{T}}}()) == 1

        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.EqualTo{T}}}()) == [bridge.indic_true]
        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ZERO, CP.DifferentFrom{T}}}()) == [bridge.indic_false]
    end

    @testset "Constraint: indicator if true" begin
        @test MOI.is_valid(model, bridge.indic_true)
        f = MOI.get(model, MOI.ConstraintFunction(), bridge.indic_true)
        @test length(f.terms) == 2
        @test MOI.get(model, MOI.ConstraintSet(), bridge.indic_true) == MOI.Indicator{MOI.ACTIVATE_ON_ONE}(MOI.EqualTo(zero(T)))
        @test f.constants == [one(T), zero(T)]

        t1 = f.terms[1]
        @test t1.output_index == 1
        @test t1.scalar_term.coefficient === -one(T)
        @test t1.scalar_term.variable == x

        t2 = f.terms[2]
        @test t2.output_index == 2
        @test t2.scalar_term.coefficient === one(T)
        @test t2.scalar_term.variable == y
    end

    @testset "Constraint: indicator if false" begin
        @test MOI.is_valid(model, bridge.indic_false)
        f = MOI.get(model, MOI.ConstraintFunction(), bridge.indic_false)
        @test length(f.terms) == 2
        @test MOI.get(model, MOI.ConstraintSet(), bridge.indic_false) == MOI.Indicator{MOI.ACTIVATE_ON_ZERO}(CP.DifferentFrom(zero(T)))
        @test f.constants == [one(T), zero(T)]

        t1 = f.terms[1]
        @test t1.output_index == 1
        @test t1.scalar_term.coefficient === -one(T)
        @test t1.scalar_term.variable == x

        t2 = f.terms[2]
        @test t2.output_index == 2
        @test t2.scalar_term.coefficient === one(T)
        @test t2.scalar_term.variable == y
    end
end


text:
module FactorLoadingMatrices
using LinearAlgebra

export nnz_loading, loading_matrix, varimax

"""
Gives the number of nonzero entries (i.e. the number of entries in the lower
triangle) of the loading matrix for the specified data dimension `nx` and number
of factors `nfactor`.
"""
function nnz_loading(nx, nfactor)
    nx >= nfactor || throw(ArgumentError("`nfactor` ($nfactor) must be <= `nx` ($nx)"))
    return ((2nx + 1) * nfactor - nfactor^2) ÷ 2
end

"""
Constructs a matrix of loadings to map `nfactor` latent variables to `nx` observed
variables.  The upper triangle is all zeros to enfactororce linear independence among the loading
vectors.
# Arguments
- `values::AbstractVector`: vector of values to put in the nonzero lower triangle. They are
filled in order running down the columns from left to right.
- `nx::Integer`: Dimension of the data, i.e. the number of rows in the loading matrix
- `nfactor::Integer`: Number of factors, i.e. the number of columns in the loading matrix
"""
function loading_matrix(values::AbstractVector{T}, nx::Integer, nfactor::Integer) where T
    if length(values) != nnz_loading(nx, nfactor)
        throw(ArgumentError("Wrong number of data values for a loading matrix of the specified size."))
    end
    L = zeros(T, nx, nfactor)
    k = 1
	for i in 1:nx, j in 1:min(i, nfactor)
		L[i, j] = values[k]
		k += 1
	end
    return L
end

"""
	varimax(A; gamma = 1.0, minit = 20, maxit = 1000, reltol = 1e-12)
VARIMAX perform varimax (or quartimax, equamax, parsimax) rotation to the column vectors of the input matrix.
# Input Arguments
- `A::Matrix{TA}`: input matrix, whose column vectors are to be rotated. d, m = size(A).
- `gamma`: default is 1. gamma = 0, 1, m/2, and d(m - 1)/(d + m - 2), corresponding to quartimax, varimax, equamax, and parsimax.
- `minit::Int`: default is 20. Minimum number of iterations, in case of the stopping criteria fails initially.
- `maxit::Int`: default is 1000. Maximum number of iterations.
- `reltol::Float64`: default is 1e-12. Relative tolerance for stopping criteria.
# Output Argument
- `B::Matrix{Float64}`: output matrix, whose columns are already been rotated.
Implemented by Haotian Li, Aug. 20, 2019
"""
function varimax(A::AbstractMatrix{TA}; gamma = 1.0, minit = 20, maxit = 1000,
		reltol = 1e-12) where TA
	d, m = size(A)
	m == 1 && return A
	# Warm up step: start with a good initial orthogonal matrix T by SVD and QR
	T = Matrix{TA}(I, m, m)
	B = A * T
	L,_,M = svd(A' * (d*B.^3 - gamma*B * Diagonal(sum(B.^2, dims = 1)[:])))
	T = L * M'
	if norm(T - Matrix{TA}(I, m, m)) < reltol
		T,_ = qr(randn(m,m)).QT
		B = A * T
	end

	# Iteration step: get better T to maximize the objective (as described in Factor Analysis book)
	D = 0
	for k in 1:maxit
		Dold = D
		L,s,M = svd(A' * (d*B.^3 - gamma*B * Diagonal(sum(B.^2, dims = 1)[:])))
		T = L * M'
		D = sum(s)
		B = A * T
		if (abs(D - Dold)/D < reltol) && k >= minit
			break
		end
	end

	# Adjust the sign of each rotated vector such that the maximum absolute value is positive.
	for i in 1:m
		if abs(maximum(B[:,i])) < abs(minimum(B[:,i]))
			B[:,i] .= - B[:,i]
		end
	end

	return B
end

end # module