text:
"""
```
estimate(m, data; verbose=:low, proposal_covariance=Matrix())
```

Estimate the DSGE parameter posterior distribution.

### Arguments:
- `m::AbstractDSGEModel` or `m::AbstractVARModel`: model object

### Estimation Settings
Please see the section on 'Estimation Settings' on the 'Advanced Usage'
page of the online documentation
or src/defaults.jl for a full description of
all the estimation settings for both
the Metropolis-Hastings (MH) and Sequential Monte Carlo (SMC) algorithms.
For the latter, also see `?DSGE.smc2`. Most of the
optional and keyword arguments described below are
not directly related to the behavior of the sampling algorithms
(e.g. tuning, number of samples).

### Optional Arguments:
- `data`: well-formed data as `Matrix` or `DataFrame`. If this is not provided, the `load_data` routine will be executed.

### Keyword Arguments:
- `verbose::Symbol`: The desired frequency of function progress messages printed to standard out.
   - `:none`: No status updates will be reported.
   - `:low` (default): Status updates will be provided in csminwel and at each block in
     Metropolis-Hastings. For SMC's verbosity settings, see `?smc2`.
   - `:high`: Status updates provided at each iteration in Metropolis-Hastings.
- `proposal_covariance::Matrix = []`: Used to test the metropolis_hastings algorithm with a precomputed
  covariance matrix for the proposal distribution. When the Hessian is singular,
  eigenvectors corresponding to zero eigenvectors are not well defined, so eigenvalue
  decomposition can cause problems. Passing a precomputed matrix allows us to ensure that
  the rest of the routine has not broken.
- `method::Symbol`: The method to use when sampling from the posterior distribution. Can
    be either `:MH` for standard Metropolis Hastings Markov Chain Monte Carlo, or `:SMC`
    for Sequential Monte Carlo. This should be specified by the setting `sampling_method` in `m`.
- `mle = false`: Set to true if parameters should be estimated by maximum likelihood directly.
    If this is set to true, this function will return after estimating parameters.
- `sampling = true`: Set to false to disable sampling from the posterior.
- `old_data::Matrix{Float64}` = []: A matrix containing the time series of observables of previous data
    (with `data` being the new data) for the purposes of a time tempered estimation
    (that is, using the posterior draws from a previous estimation as the initial set
    of draws for an estimation with new data). Running a bridge estimation
    requires both `old_data` and `old_cloud`.
- `old_cloud::Union{DSGE.ParticleCloud, DSGE.Cloud, SMC.Cloud} = DSGE.Cloud(m, 0)`: old Cloud object
    used to describe particles from a previous estimation with old data. Running a bridge estimation
    requires both `old_data` and `old_cloud`. If running a bridge estimation and
    no `old_cloud` is provided, then it will be loaded using
    the filepaths in `m`. If no `old_cloud` exists, then the bridge estimation will not run.
- `old_model::Union{AbstractDSGEModel, AbstractVARModel} = m`: model object from which we can build
    the old log-likelihood function for a time tempered SMC estimation. It should be possible
    to evaluate the old log-likelihood given `old_data` and the current draw of parameters.
    This may be nontrivial if, for example, *new* parameters have been added to `m` since the old
    estimation. In this case, `old_model` should include the *new* parameters but still return
    the old log-likelihood as the original estimation if given the same *old* parameters.
    By default, we assume the log-likelihood function has not changed
    and therefore coincides with the current one.
- `filestring_addl::Vector{String} = []`: Additional strings to add to the file name
    of estimation output as a way to distinguish output from each other.
- `continue_intermediate::Bool = false`: set to true if the estimation is starting
    from an intermediate stage that has been previously saved.
- `intermediate_stage_start::Int = 0`: number of the stage from which the user wants
    to continue the estimation (see `continue_intermediate`)
- `save_intermediate::Bool = true`: set to true to save intermediate stages when using SMC
- `intermediate_stage_increment::Int = 10`: number of stages that must pass before saving
    another intermediate stage.
-  `run_csminwel::Bool = true`: by default, csminwel is run after a SMC estimation finishes
    to recover the true mode of the posterior. Set to false to avoid this step
    (csminwel can take hours for medium-scale DSGE models).
-  `toggle::Bool = true`: when regime-switching, several functions assume regimes
    are toggled to regime 1. If the likelihood function is not
    written to toggle to regime 1 when done, then regime-switching estimation
    will not work properly. Set to `false` to reduce computation time if the
    user is certain that the likelihood is written properly.
- `log_prob_old_data::Float64 = 0.0`:Log p(\tilde y) which is the log marginal data density
    of the bridge estimation.
"""
function estimate(m::Union{AbstractDSGEModel,AbstractVARModel}, df::DataFrame;
                  verbose::Symbol = :low,
                  proposal_covariance::Matrix = Matrix(undef, 0, 0),
                  mle::Bool = false,
                  sampling::Bool = true,
                  filestring_addl::Vector{String} = Vector{String}(),
                  old_data::Matrix{Float64} = Matrix{Float64}(undef, size(data, 1), 0),
                  old_cloud::Union{DSGE.ParticleCloud, DSGE.Cloud,
                                   SMC.Cloud} = DSGE.ParticleCloud(m, 0),
                  old_model::Union{AbstractDSGEModel, AbstractVARModel} = m,
                  continue_intermediate::Bool = false,
                  intermediate_stage_start::Int = 0,
                  intermediate_stage_increment::Int = 10,
                  save_intermediate::Bool = false,
                  run_csminwel::Bool = true,
                  toggle::Bool = true,
                  log_prob_old_data::Float64 = 0.0)
    data = df_to_matrix(m, df)
    estimate(m, data; verbose = verbose, proposal_covariance = proposal_covariance,
             mle = mle, sampling = sampling,
             old_data = old_data, old_cloud = old_cloud,
             old_model = old_model,
             continue_intermediate = continue_intermediate,
             intermediate_stage_increment = intermediate_stage_increment,
             save_intermediate = save_intermediate,
             run_csminwel = run_csminwel, toggle = toggle, log_prob_old_data = log_prob_old_data)
end

function estimate(m::Union{AbstractDSGEModel,AbstractVARModel};
                  verbose::Symbol = :low,
                  proposal_covariance::Matrix = Matrix(undef, 0, 0),
                  mle::Bool = false,
                  sampling::Bool = true,
                  filestring_addl::Vector{String} = Vector{String}(),
                  old_data::Matrix{Float64} = Matrix{Float64}(undef, size(data, 1), 0),
                  old_cloud::Union{DSGE.ParticleCloud, DSGE.Cloud,
                                   SMC.Cloud} = DSGE.ParticleCloud(m, 0),
                  old_model::Union{AbstractDSGEModel, AbstractVARModel} = m,
                  continue_intermediate::Bool = false,
                  intermediate_stage_start::Int = 0,
                  intermediate_stage_increment::Int = 10,
		          save_intermediate::Bool = false,
                  run_csminwel::Bool = true,
                  toggle::Bool = true, log_prob_old_data::Float64 = 0.0)
    # Load data
    df = load_data(m; verbose = verbose)
    estimate(m, df; verbose = verbose, proposal_covariance = proposal_covariance,
             mle = mle, sampling = sampling,
             old_data = old_data, old_cloud = old_cloud,
             old_model = old_model,
             continue_intermediate = continue_intermediate,
             intermediate_stage_increment = intermediate_stage_increment,
	         save_intermediate = save_intermediate,
             run_csminwel = run_csminwel, toggle = toggle, log_prob_old_data = log_prob_old_data)
end

function estimate(m::Union{AbstractDSGEModel,AbstractVARModel}, data::AbstractArray;
                  verbose::Symbol = :low,
                  proposal_covariance::Matrix = Matrix(undef, 0,0),
                  mle::Bool = false,
                  sampling::Bool = true,
                  filestring_addl::Vector{String} = Vector{String}(),
                  old_data::Matrix{Float64} = Matrix{Float64}(undef, size(data, 1), 0),
                  old_cloud::Union{DSGE.ParticleCloud, DSGE.Cloud,
                                   SMC.Cloud} = DSGE.ParticleCloud(m, 0),
                  old_model::Union{AbstractDSGEModel, AbstractVARModel} = m,
                  continue_intermediate::Bool = false,
                  intermediate_stage_start::Int = 0,
                  intermediate_stage_increment::Int = 10,
		          save_intermediate::Bool = false,
                  run_csminwel::Bool = true,
                  toggle::Bool = true, log_prob_old_data::Float64 = 0.0)

    if !(get_setting(m, :sampling_method) in [:SMC, :MH])
        error("method must be :SMC or :MH")
    else
        method = get_setting(m, :sampling_method)
    end

    regime_switching = haskey(get_settings(m), :regime_switching) &&
        get_setting(m, :regime_switching)

    ########################################################################################
    ### Step 1: Find posterior/likelihood mode (if reoptimizing, run optimization routine)
    ########################################################################################

    # Specify starting mode

    vint = get_setting(m, :data_vintage)
    if reoptimize(m) && method == :MH
        println("Reoptimizing...")

        # Inputs to optimization algorithm
        n_iterations       = get_setting(m, :optimization_iterations)
        ftol               = get_setting(m, :optimization_ftol)
        xtol               = get_setting(m, :optimization_xtol)
        gtol               = get_setting(m, :optimization_gtol)
        step_size          = get_setting(m, :optimization_step_size)
        converged          = false

        # If the algorithm stops only because we have exceeded the maximum number of
        # iterations, continue improving guess of modal parameters
        total_iterations = 0
        optimization_time = 0
        max_attempts = get_setting(m, :optimization_attempts)
        attempts = 1

        while !converged
            begin_time = time_ns()
            out, H = optimize!(m, data;
                               method = get_setting(m, :optimization_method),
                               ftol = ftol, grtol = gtol, xtol = xtol,
                               iterations = n_iterations, show_trace = true, step_size = step_size,
                               mle = mle, toggle = toggle, verbose = verbose)

            attempts += 1
            total_iterations += out.iterations
            converged = out.converged || attempts > max_attempts

            end_time = (time_ns() - begin_time)/1e9
            println(verbose, :low, @sprintf "Total iterations completed: %d\n" total_iterations)
            println(verbose, :low, @sprintf "Optimization time elapsed: %5.2f\n" optimization_time += end_time)

            # Write params to file after every `n_iterations` iterations
            params = ModelConstructors.get_values(get_parameters(m); regime_switching = regime_switching)
            h5open(rawpath(m, "estimate", "paramsmode.h5"),"w") do file
                file["params"] = params
            end
        end

        # write parameters to file one last time so we have the final mode
        h5open(rawpath(m, "estimate", "paramsmode.h5"),"w") do file
            file["params"] = params
        end
    end

    params = ModelConstructors.get_values(get_parameters(m); regime_switching = regime_switching)

    # Sampling does not make sense if mle=true
    if mle || !sampling
        return nothing
    end

    if get_setting(m,:sampling_method) == :MH
        ########################################################################################
        ### Step 2: Compute proposal distribution for Markov Chain Monte Carlo (MCMC)
        ###
        ### In Metropolis-Hastings, we draw sample parameter vectors from
        ### the proposal distribution, which is a degenerate multivariate
        ### normal centered at the mode. Its variance is the inverse of
        ### the hessian. We find the inverse via eigenvalue decomposition.
        ########################################################################################

        ## Calculate the Hessian at the posterior mode
        hessian = if calculate_hessian(m)
            println(verbose, :low, "Recalculating Hessian...")

            hessian, _ = hessian!(m, params, data; toggle = toggle, verbose = verbose)

            h5open(rawpath(m, "estimate","hessian.h5"),"w") do file
                file["hessian"] = hessian
            end

            hessian

        ## Read in a pre-calculated Hessian
        else
            fn = hessian_path(m)
            println(verbose, :low, "Using pre-calculated Hessian from $fn")

            hessian = h5open(fn,"r") do file
                read(file, "hessian")
            end

            hessian
	    end

        # Compute inverse hessian and create proposal distribution, or
        # just create it with the given cov matrix if we have it
        propdist = if isempty(proposal_covariance)
            # Make sure the mode and hessian have the same number of parameters
            n = length(params)
            @assert (n, n) == size(hessian)

            # Compute the inverse of the Hessian via eigenvalue decomposition
            #S_diag, U = eigen(hessian)
            F = svd(hessian)

            big_eig_vals = findall(x -> x > 1e-6, F.S)
            hessian_rank = length(big_eig_vals)

            S_inv = zeros(n, n)
            #for i = (n-hessian_rank+1):n
            for i = 1:hessian_rank
                S_inv[i, i] = 1/F.S[i]
            end

            #hessian_inv = U*sqrt.(S_inv) # this is the inverse of the hessian
            hessian_inv = F.V * S_inv * F.U'#sqrt.(S_inv) * F.U'
            DegenerateMvNormal(params, hessian_inv, hessian, diag(S_inv))
        else
            DegenerateMvNormal(params, proposal_covariance, pinv(proposal_covariance),
                               eigen(proposal_covariance).values)
        end

        if rank(propdist) != n_parameters_free(m)
            println("problem –    shutting down dimensions")
        end

        ########################################################################################
        ### Step 3: Sample from posterior using Metropolis-Hastings algorithm
        ########################################################################################

        # Set the jump size for sampling
        cc0 = get_setting(m, :mh_cc0)
        cc  = get_setting(m, :mh_cc)

        metropolis_hastings(propdist, m, data, cc0, cc; regime_switching = regime_switching,
                            toggle = toggle, verbose = verbose, filestring_addl = filestring_addl);

    elseif get_setting(m, :sampling_method) == :SMC
        ########################################################################################
        ### Step 3: Run Sequential Monte Carlo (SMC)
        ###
        ### In Sequential Monte Carlo, a large number of Markov Chains
        ### are simulated iteratively to create a particle approximation
        ### of the posterior. Portions of this method are executed in
        ### parallel.
        ########################################################################################
        smc2(m, data; verbose = verbose, filestring_addl = filestring_addl,
             old_data = old_data, old_cloud = old_cloud,
             old_model = old_model,
             continue_intermediate = continue_intermediate,
             intermediate_stage_start = intermediate_stage_start,
             save_intermediate = save_intermediate,
             intermediate_stage_increment = intermediate_stage_increment,
             run_csminwel = run_csminwel,
             regime_switching = regime_switching, log_prob_old_data = log_prob_old_data)
    end

    ########################################################################################
    ### Step 4: Calculate and save parameter covariance matrix
    ########################################################################################

    compute_parameter_covariance(m, filestring_addl = filestring_addl)

    return nothing
end

"""
```
compute_parameter_covariance(m::Union{AbstractDSGEModel,AbstractVARModel})
```

Calculates the parameter covariance matrix from saved parameter draws, and writes it to the
parameter_covariance.h5 file in the `workpath(m, "estimate")` directory.

### Arguments
* `m::Union{AbstractDSGEModel,AbstractVARModel}`: the model object
- `filestring_addl::Vector{String} = []`: Additional strings to add to the file name
    of estimation output as a way to distinguish output from each other.
"""
function compute_parameter_covariance(m::Union{AbstractDSGEModel,AbstractVARModel};
                                      filestring_addl::Vector{String} = Vector{String}(undef, 0))

    sampling_method = get_setting(m, :sampling_method)
    if sampling_method ∉ [:MH, :SMC]
        throw("Invalid sampling method specified in setting :sampling_method")
    end

    prefix           = sampling_method == :MH ? "mh" : "smc"
    param_draws_path = rawpath(m, "estimate", prefix * "save.h5", filestring_addl)
    savepath         = workpath(m, "estimate", "parameter_covariance.h5", filestring_addl)

    return compute_parameter_covariance(param_draws_path, sampling_method; savepath = savepath)
end

"""
```
compute_parameter_covariance(param_draws_path::String, sampling_method::Symbol;
                             savepath = "parameter_covariance.h5")
```

Generic function calculates the parameter covariance matrix from saved parameter draws,
writes it to the specified savepath.

### Arguments
* `param_draws_path::String`: Path to file where parameter draws are stored.
* `sampling_method::Symbol`: Sampling method used for estimation.
```
   - `:MH`: Metropolis-Hastings
   - `:SMC`: Sequential Monte Carlo
```
### Optional Arguments
* `savepath::String`: Where parameter covariance matrix is to be saved. Will default
    to "parameter_covariance.h5" if unspecified.
"""
function compute_parameter_covariance(param_draws_path::String, sampling_method::Symbol;
                                      savepath::String = "parameter_covariance.h5")
    if sampling_method ∉ [:MH, :SMC]
        throw("Invalid sampling method specified in setting :sampling_method")
    end
    prefix = sampling_method == :MH ? "mh" : "smc"

    if !isfile(param_draws_path)
        @printf stderr "Saved parameter draws not found.\n"
        return
    end

    param_draws = h5open(param_draws_path, "r") do f
        read(f, prefix * "params")
    end

    # Calculate covariance matrix
    param_covariance = cov(param_draws)

    # Write to file
    h5open(savepath, "w") do f
        f[prefix * "cov"] = param_covariance
    end
end

"""
```
get_estimation_output_files(m)
```

Returns a `Dict{Symbol, String}` with all files created by `estimate(m)`.
"""
function get_estimation_output_files(m::Union{AbstractDSGEModel,AbstractVARModel})
    output_files = Dict{Symbol, String}()

    for file in [:paramsmode, :hessian, :mhsave]
        output_files[file] = rawpath(m, "estimate", "$file.h5")
    end

    for file in [:paramsmean, :parameter_covariance]
        output_files[file] = workpath(m, "estimate", "$file.h5")
    end

    for file in [:priors, :prior_posterior_means, :moments]
        output_files[file] = tablespath(m, "estimate", "$file.tex")
    end

    return output_files
end


text:
using RandomMatrices
using Test

@testset "TracyWidom" begin

# CDF lies in correct range
@test all(i->(0<i<1), cdf(TracyWidom, randn(5)))

#Test far outside support
@test cdf(TracyWidom, -10.1) ≈ 0 atol=1e-14
@test cdf(TracyWidom, 10.1) ≈ 1 atol=1e-14

# Test exact values
# See https://arxiv.org/abs/0904.1581
@test cdf(TracyWidom,0,beta=1) ≈ 0.83190806620295 atol=1e-14
@test cdf(TracyWidom,0,beta=2) ≈ 0.96937282835526 atol=1e-14

@test isfinite(rand(TracyWidom, 10))
@test isfinite(rand(TracyWidom, 100))

end # testset

ext: jl


text:
#####################################################
#
#         STRUCTURES FOR THE OPTIONS
#
#####################################################
mutable struct Options
    x_tol::Float64
    f_tol::Float64
    g_tol::Float64
    h_tol::Float64
    f_calls_limit::Float64
    g_calls_limit::Float64
    h_calls_limit::Float64
    time_limit::Float64

    iterations::Int
    store_convergence::Bool
    show_results::Bool
    debug::Bool
    search_type::Symbol
    seed::UInt
end

"""
    Options(;
        x_tol::Real = 0.0,
        f_tol::Real = 0.0,
        g_tol::Real = 0.0,
        h_tol::Real = 0.0,
        f_calls_limit::Real = 0,
        g_calls_limit::Real = 0,
        h_calls_limit::Real = 0,
        time_limit::Real = Inf,
        iterations::Int = 1,
        store_convergence::Bool = false,
        debug::Bool = false,
        seed = rand(UInt)
    )

`Options` stores common settings for metaheuristics such as the maximum number of iterations
debug options, maximum number of function evaluations, etc.

Main properties:

- `x_tol` tolerance to the true minimizer if specified in `Information`.
- `f_tol` tolerance to the true minimum if specified in `Information`.
- `f_calls_limit` is the maximum number of function evaluations limit.
- `time_limit` is the maximum time that `optimize` can spend in seconds.
- `iterations` is the maximum number iterationn permited.
- `store_convergence` if `true`, then push the current `State` in `State.convergence` at each generation/iteration
- `debug` if `true`, then `optimize` function reports the current `State` (and interest information) for each iterations.
- `seed` non-negative integer for the random generator seed.

# Example

```jldoctest
julia> options = Options(f_calls_limit = 1000, debug=true, seed=1)
Options(0.0, 0.0, 0.0, 0.0, 1000.0, 0.0, 0.0, 0, false, true, true, :minimize, 0x0000000000000001)

julia> f(x) = sum(x.^2)
f (generic function with 1 method)

julia> bounds = [  -10.0 -10 -10; # lower bounds
                    10.0  10 10 ] # upper bounds
2×3 Array{Float64,2}:
 -10.0  -10.0  -10.0
  10.0   10.0   10.0

julia> state = optimize(f, bounds, ECA(options=options))
[ Info: Initializing population...
[ Info: Starting main loop...
+=========== RESULT ==========+
| Iter.: 1
| f(x) = 6.97287
| solution.x = [-2.3628796262231875, -0.6781207370770752, -0.9642728360479853]
| f calls: 42
| Total time: 0.0004 s
+============================+

...

[ Info: Stopped since call_limit was met.
+=========== RESULT ==========+
| Iter.: 47
| f(x) = 1.56768e-08
| solution.x = [-2.2626761322304715e-5, -9.838697194048792e-5, 7.405966506272336e-5]
| f calls: 1000
| Total time: 0.0313 s
+============================+
```

"""
function Options(;
    x_tol = 0.0,
    f_tol = 0.0,
    g_tol = 0.0,
    h_tol = 0.0,
    f_calls_limit = 0.0,
    g_calls_limit = 0.0,
    h_calls_limit = 0.0,
    iterations::Int = 0,
    time_limit::Float64 = Inf,
    store_convergence = false,
    show_results = true,
    debug = false,
    search_type = :minimize,
    seed = rand(UInt)
)

    Options(
        promote(x_tol, f_tol, g_tol, h_tol)...,
        promote(f_calls_limit, g_calls_limit, h_calls_limit, time_limit)...,
        promote(iterations)...,

        # Results options
        promote(store_convergence, show_results, debug)...,
        Symbol(search_type),
        UInt(seed)
    )

end


text:
@testset "macro for generating fast functions" begin
	ε = 100*eps(Float64)
	f(x1,x2,x3) = 3*x3^2 + 2*x2*x1
	@polyvar x y z
	polynomial = 3*z^2 + 2*y*x
	g = @fastfunc polynomial

	t1,t2,t3 = randn(3)
	@test isapprox(f(t1,t2,t3),g(t1,t2,t3),atol=ɛ)
end


text:
import Base: +
"""
### parsetypes(data::AbstractString) -> T(data)
------------------
This method will turn strings of data read to their approapriate types. Notably
used by functions like getargs() in order to parse argument data into types.
Returns the data parsed into its implicit type. For reading in arguments as
specific types, see getargs(::HTTP.Stream, ::Symbol, ::Type)

"""
function parsetypes(data::AbstractString)
    x = nothing
    try
        x = parse(Int64, data)
    catch
        try
        x = parse(Bool, data)
        catch
            try
            x = parse(Float64, data)
            catch
            try
                x = parse(Array, data)
            catch
                try
                    x = parse(Dict, data)
                catch
                    x = data
                end
            end
        end
    end
end
    return(x)
end

"""
### getargs(::HTTP.Stream) -> ::Dict
------------------
The getargs method returns arguments from the HTTP header (GET requests.)
Returns a full dictionary of these values.

"""
function getargs(http::HTTP.Stream)
    target = split(http.message.target, '?')[2]
    args = split(target, '&')
    arg_dict = Dict()
    for arg in args
        keyarg = split(arg, '=')
        x = tryparse(keyarg[2])
        push!(arg_dict, Symbol(keyarg[1]) => x)
    end
    return(arg_dict)
end

"""
### getargs(::HTTP.Stream, ::Symbol) -> ::Vector
------------------
Returns the requested arguments from the target.

"""
function getarg(http::Any, s::Symbol)
    getargs(http)[s]
end

"""
### getargs(::HTTP.Stream, ::Symbol, ::Type) -> ::Vector
------------------
This method is the same as getargs(::HTTP.Stream, ::Symbol), however types are
parsed as type T(). Note that "Cannot convert..." errors are possible with this
method.

"""
function getarg(http::HTTP.Stream, s::Symbol, T::Type)
    parse(getargs(http)[s], T)
end

"""
### getpost(http::HTTP.Stream) -> _
------------------
Returns the post argument data of an HTTP stream.

"""
function getpost(http::HTTP.Stream)
    http.message.body
end

"""
### write_file(URI::String, http::HTTP.Stream) -> _
------------------
Writes a file to an HTTP.Stream.

"""
function write_file(URI::String, http::HTTP.Stream)
    open(URI, "r") do i
        write(http, i)
    end
end

"""
### lists(::Pair{String, String} ...) -> ::Vector{List}
------------------
Creates multiple Lists much more quickly. Takes pairs. keys should be
labels for the list elements, values should be the href of that button. For
    more information on List, see List.

"""
function lists(dct::Pair{String, String} ...)
    lists::Vector{List} = []
    for (key, value) in dct
        push!(lists, List(label = key, href = value))
    end
    lists
end

"""
### +(::Function, ::Function) -> ::Page
------------------
Creates a page from two function servables.

"""
+(f::Function, f2::Function) = Page([f, f2])

"""
### +(::Page, ::Function) -> ::Page
------------------
An operator binding for p.add()

"""
+(p::Page, f::Function) = p.add(f)

"""
### +(::FormComponent, ::FormComponent) -> ::Page
------------------
Easy compound forms with the + operator. See Form for more information.

"""
+(fc::FormComponent, fc2::FormComponent) = Form(fc, fc2)


text:
using Distributions, ForwardDiff
using Random: seed!
seed!(1)
Y = rand(Normal(4.0, 1.0), 100)
epochs = 200
cost(Y, μ) = loglikelihood(Normal(μ, 1.0), Y)

# adagrad
θ_adagrad = rand(Normal(), 1)
opt = Adagrad(η=1.0)

for i in 1:epochs
    g = ForwardDiff.gradient(μ -> cost(Y, μ[1]), θ_adagrad)

    δ = update(opt, g)
    θ_adagrad .+= δ
end

# rmsprop
θ_rmsprop = rand(Normal(), 1)
opt = RMSprop(η=0.1)

for i in 1:(epochs*2)
    g = ForwardDiff.gradient(μ -> cost(Y, μ[1]), θ_rmsprop)

    δ = update(opt, g)
    θ_rmsprop .+= δ
end

# adam
θ_adam = rand(Normal(), 1)
opt = Adam(α=1.0)

for i in 1:epochs
    g = ForwardDiff.gradient(μ -> cost(Y, μ[1]), θ_adam)

    δ = update(opt, g)
    θ_adam .+= δ
end

# adamax
θ_adamax = rand(Normal(), 1)
opt = Adamax(α=1.0)

for i in 1:epochs
    g = ForwardDiff.gradient(μ -> cost(Y, μ[1]), θ_adamax)

    δ = update(opt, g)
    θ_adamax .+= δ
end

## nadam
θ_nadam = rand(Normal(), 1)
opt = Nadam(η=1.0)

for i in 1:epochs
    g = ForwardDiff.gradient(μ -> cost(Y, μ[1]), θ_nadam)

    δ = update(opt, g)
    θ_nadam .+= δ
end

## vanilla
θ_vanilla = rand(Normal(), 1)
opt = VanillaGradDescent(η=0.001)

for i in 1:epochs
    g = ForwardDiff.gradient(μ -> cost(Y, μ[1]), θ_vanilla)

    δ = update(opt, g)
    θ_vanilla .+= δ
end

## inversedecay
θ_invdec = rand(Normal(), 1)
opt = Inversedecay(κ=0.9)

for i in 1:epochs
    g = ForwardDiff.gradient(μ -> cost(Y, μ[1]), θ_invdec)

    δ = update(opt, g)
    θ_invdec .+= δ
end

@testset "Normal MLE" begin
    @test mean(Y) ≈ θ_adagrad[1] atol=1e-3
    @test mean(Y) ≈ θ_rmsprop[1] atol=1e-1
    @test mean(Y) ≈ θ_adam[1] atol=1e-3
    @test mean(Y) ≈ θ_adamax[1] atol=1e-3
    @test mean(Y) ≈ θ_nadam[1] atol=1e-3
    @test mean(Y) ≈ θ_vanilla[1] atol=1e-3
    @test mean(Y) ≈ θ_invdec[1] atol=1e-3
end


text:
using QuasiCopula, LinearAlgebra, GLM
using Random, Distributions, DataFrames, ToeplitzMatrices
using Test, BenchmarkTools
BLAS.set_num_threads(1)
Threads.nthreads()

p = 3    # number of fixed effects, including intercept

# true parameter values
Random.seed!(12345)
# try next
βtrue = rand(Uniform(-2, 2), p)
σ2true = [0.5]
ρtrue = [0.5]

function get_V(ρ, n)
  vec = zeros(n)
  vec[1] = 1.0
  for i in 2:n
      vec[i] = vec[i - 1] * ρ
  end
  V = ToeplitzMatrices.SymmetricToeplitz(vec)
  V
end

trueparams = [βtrue; ρtrue; σ2true] #hold true parameters

#simulation parameters
samplesize = 10000
ni = 25

st = time()
currentind = 1
d = Bernoulli()
link = LogitLink()
D = typeof(d)
Link = typeof(link)
T = Float64

gcs = Vector{GLMCopulaARObs{T, D, Link}}(undef, samplesize)

V = get_V(ρtrue[1], ni)

# true Gamma
Γ = σ2true[1] * V

# for reproducibility I will simulate all the design matrices here
Random.seed!(12345)
X_samplesize = [randn(ni, p - 1) for i in 1:samplesize]

for i in 1:samplesize
  X = [ones(ni) X_samplesize[1]]
  η = X * βtrue
  μ = exp.(η) ./ (1 .+ exp.(η))
  vecd = Vector{DiscreteUnivariateDistribution}(undef, ni)
  for i in 1:ni
      vecd[i] = Bernoulli(μ[i])
  end
  nonmixed_multivariate_dist = NonMixedMultivariateDistribution(vecd, Γ)
  # simuate single vector y
  y = Vector{Float64}(undef, ni)
  res = Vector{Float64}(undef, ni)
  rand(nonmixed_multivariate_dist, y, res)
  V = [ones(ni, ni)]
  gcs[i] = GLMCopulaARObs(y, X, d, link)
end

# form model
gcm = GLMCopulaARModel(gcs)
# precompile
println("precompiling Bernoulli AR fit")
gcm2 = deepcopy(gcm);
QuasiCopula.fit!(gcm2, IpoptSolver(print_level = 0, max_iter = 20));

fittime = @elapsed QuasiCopula.fit!(gcm)
@show fittime
@show gcm.β
@show gcm.σ2
@show gcm.ρ
@show gcm.∇β
@show gcm.∇σ2
@show gcm.∇ρ

@test logl(gcm) == loglikelihood!(gcm, false, false)
@show get_CI(gcm)

# mse and time under our model
mseβ, mseρ, mseσ2 = MSE(gcm, βtrue, ρtrue, σ2true)
@show mseβ
@show mseσ2
@show mseρ

@test mseβ < 0.01
@test mseσ2 < 1
@test mseρ < 0.01

println("checking memory allocation for Bernoulli AR")
logl_gradient_memory = @benchmark loglikelihood!($gcm.data[1], $gcm.β, $gcm.ρ[1], $gcm.σ2[1], true, false)
@test logl_gradient_memory.memory == 0.0

# note for multi-threading we will allocate some at the model level loglikelihood!(gcm, true, false)


text:
#!/usr/bin/env julia
using Plots, Random, Printf

function initCondGen() #get random initial conditions for mass/radius, position, and velocity
    function getMass(nBodies) #generate random masses that better reflect actual stellar populations
        mList=zeros(nBodies)
        N=(0.5^(-1.3)-150^(-1.3))/1.3 #crude approximation of IMF integral assuming alpha = 2.3, stellar mass range of 0.5:150 solar masses
        rescale=1e6
        max=floor(Int,N*rescale)
        for i=1:nBodies
            intTarget=rand(0:max,1)[1]/rescale
            m=(0.5^(-1.3)-intTarget*1.3)^(-1/1.3) #just algebra from above
            mList[i]=round(m,digits=2)
        end
        return mList
    end
    m=rand(1:1500,3)./10 #3 random masses between 0.1 and 150 solar masses, uniform distribution
    #m=getMass(3) #get mass from IMF -- this way is kind of boring...so not using it, but left here in case I change my mind?
    rad=m.^0.8 #3 radii based on masses in solar units
    m=m.*2e30 #convert to SI kg
    rad=rad.*7e8 #convert to SI m
    pos1=rand(-10:10,3) #random initial coordinates x & y for first body, AU
    function genPos2(pos1)
        accept2=false
        while accept2==false
            pos2=rand(-10:10,3) #random initial coordinates for second body, AU
            dist21=sqrt((pos1[1]-pos2[1])^2+(pos1[2]-pos2[2])^2+(pos1[3]-pos2[3])^2)
            if (dist21*1.5e11)>(rad[1]+rad[2]) #they aren't touching
                accept2=true
                return pos2
            end
        end
    end
    pos2=genPos2(pos1)
    function genPos3(pos1,pos2)
        accept3=false
        while accept3==false
            pos3=rand(-10:10,3) #random initial coordinates for third body, AU
            dist31=sqrt((pos1[1]-pos3[1])^2+(pos1[2]-pos3[2])^2+(pos1[3]-pos3[3])^2)
            dist32=sqrt((pos2[1]-pos3[1])^2+(pos2[2]-pos3[2])^2+(pos2[3]-pos3[3])^2)
            if (dist31*1.5e11)>(rad[1]+rad[3]) && (dist32*1.5e11)>(rad[2]+rad[3]) #3rd isn't touching either
                accept3=true
                return pos3
            end
        end
    end
    pos3=genPos3(pos1,pos2)
    pos=[pos1[1],pos1[2],pos1[3],pos2[1],pos2[2],pos2[3],pos3[1],pos3[2],pos3[3]].*1.5e11 #convert accepted positions to SI, m
    v=rand(-7e3:7e3,9) #random xyz velocities with mag between -10 & 10 km/s, totally arbitrary...
    #r=[x1,y1,x2,y2,x3,y3,v1x,v1y,v2x,v2y,v3x,v3y]
    r=[pos[1],pos[2],pos[3],pos[4],pos[5],pos[6],pos[7],pos[8],pos[9],v[1],v[2],v[3],v[4],v[5],v[6],v[7],v[8],v[9]]
    open("initCond.txt","w") do f #save initial conditions to file in folder where script is run
        write(f,"m1=$(@sprintf("%.1f",(m[1]/2e30))) m2=$(@sprintf("%.1f",(m[2]/2e30))) m3=$(@sprintf("%.1f",(m[3]/2e30))) (solar masses)\nv1x=$(v[1]/1e3) v1y=$(v[2]/1e3) v1z=$(v[3]/1e3) v2x=$(v[4]/1e3) v2y=$(v[5]/1e3) v2z=$(v[6]/1e3) v3x=$(v[7]/1e3) v3y=$(v[8]/1e3) v3z=$(v[9]/1e3) (km/s)\nx1=$(pos1[1]) y1=$(pos1[2]) z1=$(pos1[3]) x2=$(pos2[1]) y2=$(pos2[2]) z2=$(pos2[3]) x3=$(pos3[1]) y3=$(pos3[2]) z3=$(pos3[3]) (AU from center)")
    end
    return r, rad, m
end

function dR(r,m) #function we will use RK4 on to approximate solution
    G=6.67408313131313e-11# Nm^2/kg^2
    M1,M2,M3=m[1],m[2],m[3] #kg
    x1,x2,x3=r[1],r[4],r[7] #m
    y1,y2,y3=r[2],r[5],r[8] #m
    z1,z2,z3=r[3],r[6],r[9]
    c1,c2,c3=G*M1,G*M2,G*M3 #Nm^2/kg
    r1_2=sqrt((x1-x2)^2+(y1-y2)^2+(z1-z2)^2) #distance from 1->2, m
    r1_3=sqrt((x1-x3)^2+(y1-y3)^2+(z1-z3)^2) #distance from 1->3, m
    r2_3=sqrt((x2-x3)^2+(y2-y3)^2+(z2-z3)^2) #distance from 2->3, m

    v1X,v2X,v3X=r[10],r[13],r[16] #these are our change in position after dt (dr/dt*dt=dr)
    v1Y,v2Y,v3Y=r[11],r[14],r[17] #m after * dt
    v1Z,v2Z,v3Z=r[12],r[15],r[18]

    #get change in velocity from accelerations (d^2r/dt^2*dt=dv/dt*dt=dv)
    dx1=-(c2*(x1-x2)/(r1_2^3))-(c3*(x1-x3)/(r1_3^3)) #d^2x/dt^2 for 1, m/s after * dt
    dx2=-(c1*(x2-x1)/(r1_2^3))-(c3*(x2-x3)/(r2_3^3)) #d^2x/dt^2 for 2, m/s
    dx3=-(c1*(x3-x1)/(r1_3^3))-(c2*(x3-x2)/(r2_3^3)) #d^2x/dt^2 for 3, m/s
    dy1=-(c2*(y1-y2)/(r1_2^3))-(c3*(y1-y3)/(r1_3^3)) #d^2y/dt^2 for 1, m/s
    dy2=-(c1*(y2-y1)/(r1_2^3))-(c3*(y2-y3)/(r2_3^3)) #d^2y/dt^2 for 2, m/s
    dy3=-(c1*(y3-y1)/(r1_3^3))-(c2*(y3-y2)/(r2_3^3)) #d^2y/dt^2 for 3, m/s
    dz1=-(c2*(z1-z2)/(r1_2^3))-(c3*(z1-z3)/(r1_3^3)) #d^2y/dt^2 for 1, m/s
    dz2=-(c1*(z2-z1)/(r1_2^3))-(c3*(z2-z3)/(r2_3^3)) #d^2y/dt^2 for 2, m/s
    dz3=-(c1*(z3-z1)/(r1_3^3))-(c2*(z3-z2)/(r2_3^3)) #d^2y/dt^2 for 3, m/s

    return [v1X,v1Y,v1Z,v2X,v2Y,v2Z,v3X,v3Y,v3Z,dx1,dy1,dz1,dx2,dy2,dz2,dx3,dy3,dz3]
end

function gen3Body(stopCond=[10,100],numSteps=10000) #default stop conditions of 10 yrs and 100 AU sep
    tStop=stopCond[1]*365*24*3600 #convert to SI s
    sepStop=stopCond[2]*1.5e11 #convert to SI m
    stop=false
    currentT=0
    t=range(0,stop=tStop,length=(numSteps+1)) #+1 because I don't want 0 to count
    stepSize=tStop/numSteps
    x1=zeros(length(t))
    y1=zeros(length(t))
    z1=zeros(length(t))
    x2=zeros(length(t))
    y2=zeros(length(t))
    z2=zeros(length(t))
    x3=zeros(length(t))
    y3=zeros(length(t))
    z3=zeros(length(t))
    r,rad,m=initCondGen()
    min12=rad[1]+rad[2]
    min13=rad[1]+rad[3]
    min23=rad[2]+rad[3]
    i=1
    stopT=maximum(t)

    #implement RK4 to model solutions to differential equations
    while stop==false
        if currentT==stopT || currentT>stopT #in case of rounding error or something
            stop=true
        elseif i>numSteps+1 #inf loop failsafe
            stop=true
            println("error: shouldn't have gotten here")
        else
            x1[i]=r[1]
            y1[i]=r[2]
            z1[i]=r[3]
            x2[i]=r[4]
            y2[i]=r[5]
            z2[i]=r[6]
            x3[i]=r[7]
            y3[i]=r[8]
            z3[i]=r[9]

            k1=stepSize*dR(r,m)
            k2=stepSize*dR(r.+0.5.*k1,m)
            k3=stepSize*dR(r.+0.5.*k2,m)
            k4=stepSize*dR(r.+k3,m)
            r+=(k1.+2.0*k2.+2.0.*k3.+k4)./6

            #check separation after each dt step
            sep12=sqrt((x1[i]-x2[i])^2+(y1[i]-y2[i])^2+(z1[i]-z2[i])^2)
            sep13=sqrt((x1[i]-x3[i])^2+(y1[i]-y3[i])^2+(z1[i]-z3[i])^2)
            sep23=sqrt((x3[i]-x2[i])^2+(y3[i]-y2[i])^2+(z3[i]-z2[i])^2)

            if sep12<min12 || sep13<min13 || sep23<min23 || sep12>sepStop || sep13>sepStop || sep23>sepStop
                stop=true #stop if collision happens or body is ejected
                t=range(0,stop=currentT,length=i) #t should match pos vectors
                x1=x1[1:i] #don't want trailing zeros
                y1=y1[1:i]
                z1=z1[1:i]
                x2=x2[1:i]
                y2=y2[1:i]
                z2=z2[1:i]
                x3=x3[1:i]
                y3=y3[1:i]
                z3=z3[1:i]
            end
            i+=1
            currentT+=stepSize #next step
        end
    end
    return [x1,y1,z1,x2,y2,z2,x3,y3,z3], t, m, rad
end

function getInteresting3Body(minTime=0) #in years, defaults to 0
    #sometimes random conditions result in a really short animation where things
    #just crash into each other/fly away, so this function throws away those
    yearSec=365*24*3600
    interesting=false
    i=1
    while interesting==false
        plotData,t,m,rad=gen3Body([60,150],600000)
        if (maximum(t)/yearSec)>minTime #only return if simulation runs for longer than minTime
            println(maximum(t)/yearSec) #tell me how many years we are simulating
            open("cron_log.txt","a") do f #for cron logging, a flag = append
                write(f,"$(maximum(t)/yearSec)\n")
            end
            return plotData,t,m,rad
            interesting=true
        elseif i>5000 #computationally expensive so don't want to go forever
            interesting=true #render it anyways I guess because sometimes it's fun?
            println("did not find interesting solution in number of tries allotted, running anyways")
            println(maximum(t)/yearSec) #how many years simulation runs for
            open("cron_log.txt","a") do f #for cron logging
                write(f,"$(maximum(t)/yearSec)\n")
            end
            return plotData,t,m,rad
        end
        i+=1
    end
end

function getLims(pos,padding) #determines plot limits at each frame, padding in units of pos
    x=[pos[1],pos[4],pos[7]]
    xMin=minimum(x)
    xMax=maximum(x)
    dx=xMax-xMin
    y=[pos[2],pos[5],pos[8]]
    yMin=minimum(y)
    yMax=maximum(y)
    dy=yMax-yMin
    z=[pos[3],pos[6],pos[9]]
    zMin=minimum(z)
    zMax=maximum(z)
    dz=zMax-zMin
    dList=[dx,dy,dz]
    if maximum(dList)==dx
        #use x for square
        xlims=[xMin-padding,xMax+padding]
        ylims=[yMin-padding,yMin+dx+padding]
        zlims=[zMin-padding,zMin+dx+padding]
    elseif maximum(dList)==dy
        #use y for square
        xlims=[xMin-padding,xMin+dy+padding]
        ylims=[yMin-padding,yMax+padding]
        zlims=[zMin-padding,zMin+dy+padding]
    else
        #use z for cube
        xlims=[xMin-padding,xMin+dz+padding]
        ylims=[yMin-padding,yMin+dz+padding]
        zlims=[zMin-padding,zMax+padding]
    end
    return xlims,ylims,zlims
end

function getColors(m,c) #places colors of objects according to mass/size
    #c=[:biggest,:medium,:smallest] (order of input colors)
    maxM=maximum(m)
    minM=minimum(m)
    colors=[:blue,:blue,:blue] #testing
    if m[1]==maxM
        colors[1]=c[1]
        if m[2]==minM
            colors[2]=c[3]
            colors[3]=c[2]
        else
            colors[3]=c[3]
            colors[2]=c[2]
        end
    elseif m[2]==maxM
        colors[2]=c[1]
        if m[1]==minM
            colors[1]=c[3]
            colors[3]=c[2]
        else
            colors[3]=c[3]
            colors[1]=c[2]
        end
    else
        colors[3]=c[1]
        if m[1]==minM
            colors[1]=c[3]
            colors[2]=c[2]
        else
            colors[2]=c[3]
            colors[1]=c[2]
        end
    end
    return colors
end

function makeCircleVals(r,center=[0,0,0])
    #(x,y,z)=(rcos(theta)sin(phi),rsin(theta)sin(phi),rcos(phi))
    xOffset=center[1]
    yOffset=center[2]
    zOffset=center[3]
    xVals=[]
    yVals=[]
    zVals=[]
    for i=0:pi/128:2*pi
        for j=0:pi/64:pi
            x=r*cos(i)*sin(j)+xOffset
            y=r*sin(i)*sin(j)+yOffset
            z=r*cos(j)+zOffset
            push!(xVals,x)
            push!(yVals,y)
            push!(zVals,z)
        end
    end
    return xVals,yVals,zVals
end

plotData,t,m,rad=getInteresting3Body(15)
c=[:DodgerBlue,:Gold,:Tomato] #most massive to least massive, also roughly corresponds to temp
colors=getColors(m,c)
#adding fake stars
numStars=2500
starsX=zeros(numStars)
starsY=zeros(numStars)
starsZ=zeros(numStars)
for i=1:numStars
    num=rand(-200:200,3) #box size is 70 AU but we need some extra padding for movement
    starsX[i]=num[1]
    starsY[i]=num[2]
    starsZ[i]=num[3]
end

#this new way runs significantly faster (~2x improvement over @anim)
#Downside is it spams folder with png images of every frame and must manually compile with ffmpeg
#Comment out and use older way (after this below) if performance/specific formatting is not an issue

plotLoadPath="/home/kirk/Documents/3Body/tmpPlots/"
threeBodyAnim=Animation(plotLoadPath,String[])
for i=1:333:length(t) #this makes animation scale ~1 sec/year with other conditions
    GR.inline("png") #added to eneable cron/jobber compatibility, also this makes frames generate WAY faster? Prior to adding this when run from cron/jobber frames would stop generating at 408 for some reason.
    gr(legendfontcolor = plot_color(:white)) #legendfontcolor=:white plot arg broken right now (at least in this backend)
    print("$(@sprintf("%.2f",i/length(t)*100)) % complete\r") #output percent tracker
    pos=[plotData[1][i],plotData[2][i],plotData[3][i],plotData[4][i],plotData[5][i],plotData[6][i],plotData[7][i],plotData[8][i],plotData[9][i]] #current pos
    limx,limy,limz=getLims(pos./1.5e11,10) #convert to AU, 10 AU padding
    p=plot3d(plotData[1][1:33:i]./1.5e11,plotData[2][1:33:i]./1.5e11,plotData[3][1:33:i]./1.5e11,label="",linecolor=colors[1],linealpha=max.((1:33:i) .+ 10000 .- i,2500)/10000) #plot orbits up to i
    p=plot3d!(plotData[4][1:33:i]./1.5e11,plotData[5][1:33:i]./1.5e11,plotData[6][1:33:i]./1.5e11,label="",linecolor=colors[2],linealpha=max.((1:33:i) .+ 10000 .- i,2500)/10000) #linealpha argument causes lines to decay
    p=plot3d!(plotData[7][1:33:i]./1.5e11,plotData[8][1:33:i]./1.5e11,plotData[9][1:33:i]./1.5e11,label="",linecolor=colors[3],linealpha=max.((1:33:i) .+ 10000 .- i,2500)/10000) #example: alpha=max.((1:i) .+ 100 .- i,0) causes only last 100 to be visible
    p=scatter3d!(starsX,starsY,starsZ,markercolor=:white,markersize=:1,label="") #fake background stars
    star1=makeCircleVals(rad[1],[plotData[1][i],plotData[2][i],plotData[3][i]]) #generate spheres with appropriate sizes for each star
    star2=makeCircleVals(rad[2],[plotData[4][i],plotData[5][i],plotData[6][i]]) #at current positions
    star3=makeCircleVals(rad[3],[plotData[7][i],plotData[8][i],plotData[9][i]])
    p=plot3d!(star1[1]./1.5e11,star1[2]./1.5e11,star1[3]./1.5e11,label="$(@sprintf("%.1f", m[1]./2e30))",color=colors[1],fill=true) #plot star circles with labels
    p=plot3d!(star2[1]./1.5e11,star2[2]./1.5e11,star2[3]./1.5e11,label="$(@sprintf("%.1f", m[2]./2e30))",color=colors[2],fill=true)
    p=plot3d!(star3[1]./1.5e11,star3[2]./1.5e11,star3[3]./1.5e11,label="$(@sprintf("%.1f", m[3]./2e30))",color=colors[3],fill=true)
    p=plot3d!(background_color=:black,background_color_legend=:transparent,foreground_color_legend=:transparent,
        background_color_outside=:white,aspect_ratio=:equal,legendtitlefontcolor=:white) #formatting for plot frame
    p=plot3d!(title="Random Three-Body Problem\nt: $(@sprintf("%0.2f",t[i]/365/24/3600)) years after start",
        legend=:best,xaxis=("x: AU",(limx[1],limx[2]),font(9,"Courier")),yaxis=("y: AU",(limy[1],limy[2]),font(9,"Courier")),zaxis=("z: AU",(limz[1],limz[2]),font(9,"Courier")),
        gridalpha=0.5,gridcolor=:white,titlefont=font(14,"Courier"),size=(720,721),legendfontsize=8,legendtitle="Mass (in solar masses)",legendtitlefontsize=8) #add in axes/title/legend with formatting
    frame(threeBodyAnim,p) #generate the frame
    closeall() #close plots
end

#threeBodyFile="3Body_fps30.mp4"

#crf is compression value (17 or 18 "visually lossless"), pix_fmt is for twitter specific vid req, -b:v specifies target bitrate, -vcodec specifies codec (h264 in this case) -y says overwrite existing file
#run( `ffmpeg -framerate 30 -i $plotLoadPath"%06d.png" -vcodec libx264 -pix_fmt yuv420p -profile:v high -b:v 2048K -y -vf "scale=720:720,setdar=1/1" $threeBodyFile` ) #-vf scale=720:72 -crf 25
#run( `ffmpeg -framerate 30 -i $plotLoadPath"%06d.png" -c:v libx264 -preset slow -coder 1 -movflags +faststart -g 15 -crf 18 -pix_fmt yuv420p -profile:v high -y -bf 2 -vf "scale=720:720,setdar=1/1" $threeBodyFile` ) #all this bullshit to hopefully satisfy twitter requirements

#NOTE: moved ffmpeg commands to shell script

#old (simpler) way of generating animation
#uncomment and use this way if you just want a simple animation saved and don't
#care about performance/specific formatting of video.

# threeBodyAnim=@animate for i=1:length(t)
#     gr(legendfontcolor = plot_color(:white)) #plot arg broken right now in Julia
#     print("$(@sprintf("%.2f",i/length(t)*100)) % complete\r") #output percent tracker
#     pos=[plotData[1][i],plotData[2][i],plotData[3][i],plotData[4][i],plotData[5][i],plotData[6][i]] #current pos
#     limx,limy=getLims(pos./1.5e11,5) #convert to AU, 5 AU padding
#     plot(plotData[1][1:i]./1.5e11,plotData[2][1:i]./1.5e11,label="",linecolor=colors[1])
#     plot!(plotData[3][1:i]./1.5e11,plotData[4][1:i]./1.5e11,label="",linecolor=colors[2])
#     plot!(plotData[5][1:i]./1.5e11,plotData[6][1:i]./1.5e11,label="",linecolor=colors[3])
#     scatter!(starsX,starsY,markercolor=:white,markersize=:1,label="") #fake background stars
#     star1=makeCircleVals(rad[1],[plotData[1][i],plotData[2][i]])
#     star2=makeCircleVals(rad[2],[plotData[3][i],plotData[4][i]])
#     star3=makeCircleVals(rad[3],[plotData[5][i],plotData[6][i]])
#     plot!(star1[1]./1.5e11,star1[2]./1.5e11,label="$(@sprintf("%.1f", m[1]./2e30))",color=colors[1],fill=true)
#     plot!(star2[1]./1.5e11,star2[2]./1.5e11,label="$(@sprintf("%.1f", m[2]./2e30))",color=colors[2],fill=true)
#     plot!(star3[1]./1.5e11,star3[2]./1.5e11,label="$(@sprintf("%.1f", m[3]./2e30))",color=colors[3],fill=true)
#     plot!(background_color=:black,background_color_legend=:transparent,background_color_outside=:white,aspect_ratio=:equal,legendtitlefontcolor=:white) #legendfontcolor=:white
#     plot!(xlabel="x: AU",ylabel="y: AU",title="Random Three Body Problem\nt: $(@sprintf("%0.2f",t[i]/365/24/3600)) yrs after start",
#         legend=:best,xaxis=("x: AU",(limx[1],limx[2]),font(12,"Courier")),yaxis=("y: AU",(limy[1],limy[2]),font(12,"Courier")),
#         grid=false,titlefont=font(24,"Courier"),size=(720,720),legendfontsize=12,legendtitle="Mass (in solar masses)",legendtitlefontsize=14)
#     end every 25

#mp4(threeBodyAnim,"3Body_fps30.mp4",fps=30)
#OR
#gif(threeBodyAnim,"3Body_fps30.gif",fps=30)


text:
function test_bisection_search()
    n = 500
    x0 = fill(0.0, n)
    x1 = fill(1.0, n)
    tol = 1e-5
    t = clamp.(rand(n), 1e-3, 1.0 - 1e-3)
    @info "test_bisection_search" n tol
    for e in [10.0^i for i in -1:0.1:1]
        println("f(x) = x.^$(e)")
        f(x) = x.^e
        x = SimTools.bisection_search(f, x0, x1, t; verbose = false)
        @test !isnothing(x) && all(isapprox.(f(x), t; atol = tol))
    end
end
test_bisection_search()


text:
# This file was generated, do not modify it. # hide
lm = LR()
mlm = machine(lm, select(X, :Horsepower), y)
fit!(mlm, rows=train)
rms(MLJ.predict(mlm, rows=test), y[test])^2

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
export ARC_Cub
function ARC_Cub(h :: AbstractNLPModel;
                 t₀ :: Float64 = h.meta.x0[1],
                 tₘ :: Float64 = 100.0,
                 tol :: Float64 = 1e-7,
                 maxiter :: Int = 50,
                 verbose :: Bool = false)

    (length(h.meta.x0) > 1) && warn("Not a 1-D problem ")
    # ARC parameters
    eps1 = 0.2; eps2 = 0.8
    red = 0.5; aug = 20
    Δ = 1.0
    t = t₀

    iter = 0;
    fₖ = obj(h,[t])[1]
    gₖ = grad(h, [t])[1]

    secₖ = 1.0
    dN=-gₖ #pour la première iter
    A=0.5
    B=0.0
    q(d) = fₖ + gₖ * d + A * d^2 + B * d^3

    verbose && @printf(" iter  t         gₖ          Δ        pred         ared\n")
    verbose && @printf(" %4d %7.2e  %7.2e  %7.2e \n", iter, t, gₖ, Δ)

    while ((abs(gₖ[1])>tol) & (iter < maxiter)) || (iter == 0)

      #Step computation
      Quad(t) = fₖ + gₖ * t + A * t^2 + B * t^3 + (1 / (4 * Δ)) * t^4
      dQuad(t) = gₖ + 2 * A * t + 3 * B * t^2 + (1 / Δ) * t^3

      #dR=PolynomialRoots.roots([gₖ,2*A,3*B,(1/Δ)])
      #TECHNIQUEMENT LE CALCUL DES RACINES PAR PolynomialRoots À LA MÊME PRÉCISION QUE SCILAB, MAIS IL DONNE UN SEG FAULT (comme pour TR_Cub)
      p = Poly([gₖ, 2 * A, 3 * B, (1 / Δ)])
      dR = roots(p)
      vmin=Inf
      dN=0
      for i=1:length(dR)
        rr = dR[i]
        if isreal(rr)
          rr = real(rr)
          vact = Quad(rr)
          if rr * gₖ[1] < 0.0
            if vact < vmin
              dN = rr
              vmin = vact
            end
          end
        end
      end

      d = dN

      if d == 0.0
        return (t, fₖ, norm(gₖ, Inf), iter + 1, false, false, :FailedDirectionComputation, h.counters.neval_obj, h.counters.neval_grad, h.counters.neval_hess)
        #L'OUTIL DE CALCUL DES RACINES DU PACKAGE Polynomials ARRONDIE À 0 CERTAINES RACINES D'OÙ LA NOUVELLE CONDITION
      end

      ftestTR=obj(h, [t + d])[1]
      gtestTR=grad(h, [t + d])[1]

      pred = gₖ[1] * d + A *d ^2 + B * d^3
      if pred > -1e-10
        ared=(gₖ[1] + gtestTR[1]) * d / 2
      else
        ared = ftestTR - fₖ
      end

      ratio = ared / pred
      if ratio < eps1
        Δ = red * Δ
      else
        #two point memory
        tpred = t
        gkm1 = gₖ
        fkm1 = fₖ

        t = t + d
        gₖ = gtestTR
        fₖ = ftestTR

        #Cubic two points interpolation
        s = t - tpred
        y = gₖ - gkm1
        α = -s
        z = gₖ + gkm1 + 3 *( fₖ - fkm1) / α
        discr= z^2 - gₖ * gkm1
        denom = gₖ + gkm1 + 2 * z
        B = (1 / 3) * (gₖ + gkm1 + 2 * z) / ( α * α)
        A = -(gₖ + z) / α
        if ratio > eps2
          Δ = aug * Δ
        end
      end
      iter += 1
      verbose && @printf(" %4d %7.2e  %7.2e  %7.2e %7.2e %7.2e\n", iter, t, gₖ, Δ, pred, ared)
    end
    status = :NotSolved
    (abs(gₖ) < tol) && (status = :Optimal)
    (iter >= maxiter) && (status = :Tired)
    tired = iter > maxiter
    optimal = abs(gₖ) < tol
    return (t, fₖ, norm(gₖ, Inf), iter, optimal, tired, status, h.counters.neval_obj, h.counters.neval_grad, h.counters.neval_hess)

end


text:
using ExternalBallistics
using Test

@testset "ExternalBallistics.jl" begin
    # Write your tests here.
end


text:
"""
    $SIGNATURES

Calculate the proposals and products on `destvert` using `factors` in factor graph `dfg`.

Notes
- Returns tuple of product and whether full dimensional (=true) or partial (=false).
- `N` determines the number of samples to draw from the marginal.
- `dens` can contain mixed full and partial dimension `ManifoldKernelDensity` beliefs

Related

[`approxConv`](@ref), [`proposalbeliefs!`](@ref), [`AMP.manifoldProduct`](@ref)
"""
function propagateBelief( dfg::AbstractDFG,
                          destvar::DFGVariable,
                          factors::AbstractVector; #{<:DFGFactor};
                          solveKey::Symbol = :default,
                          dens::AbstractVector{<:ManifoldKernelDensity} = Vector{ManifoldKernelDensity}(), # TODO, abstract is unstable
                          N::Integer = getSolverParams(dfg).N,
                          needFreshMeasurements::Bool = true,
                          dbg::Bool = false,
                          logger=ConsoleLogger()  )
  #

  # get proposal beliefs
  destlbl = getLabel(destvar)
  ipc = proposalbeliefs!(dfg, destlbl, factors, dens, solveKey=solveKey, N=N, dbg=dbg)

  # @show dens[1].manifold

  # make sure oldpts has right number of points
  oldBel = getBelief(dfg, destlbl, solveKey)
  oldpts = if Npts(oldBel) == N
    getPoints(oldBel)
  else
    sample(oldBel, N)[1]
  end

  # few more data requirements
  varType = getVariableType(destvar)
  M = getManifold(varType)
  # @info "BUILDING MKD" varType M isPartial.(dens)

  # take the product
  mkd = AMP.manifoldProduct(dens, M, Niter=1, oldPoints=oldpts, N=N, u0=getPointDefault(varType) )

  # @info "GOT" mkd.manifold
  return mkd, ipc
end

"""
    $SIGNATURES

This is an old function name that will be replaced by [`propagateBelief`](@ref).
"""
function predictbelief( dfg::AbstractDFG,
                        destvert::DFGVariable,
                        factors::AbstractVector; #{<:DFGFactor};
                        asPartial::Bool=false,
                        kw... )
  #
  # new
  mkd, ifd = propagateBelief(dfg,destvert,factors;kw...)

  # legacy interface
  return getPoints(mkd, asPartial), ifd
end

predictbelief(dfg::AbstractDFG,
              destlbl::Symbol,
              fctlbls::AbstractVector{Symbol};
              kw... ) = predictbelief(dfg, getVariable(dfg, destlbl), map(x->getFactor(dfg, x), fctlbls); kw... )
#

predictbelief(dfg::AbstractDFG,
              destlbl::Symbol,
              ::Colon;
              kw... ) = predictbelief(dfg, destlbl, getNeighbors(dfg, destlbl); kw... )
#

"""
    $(SIGNATURES)

Using factor graph object `dfg`, project belief through connected factors
(convolution with likelihood) to variable `sym` followed by a approximate functional product.

Return: product belief, full proposals, partial dimension proposals, labels
"""
function localProduct(dfg::AbstractDFG,
                      sym::Symbol;
                      solveKey::Symbol=:default,
                      N::Int=getSolverParams(dfg).N, #maximum([length(getPoints(getBelief(dfg, sym, solveKey))); getSolverParams(dfg).N]),
                      dbg::Bool=false,
                      logger=ConsoleLogger() )
  #
  # vector of all neighbors as Symbols
  lb = getNeighbors(dfg, sym)

  # store proposal beliefs, TODO replace Abstract with concrete type
  dens = Vector{ManifoldKernelDensity}()

  fcts = map(x->getFactor(dfg, x), lb)
  mkd, sinfd = propagateBelief(dfg, getVariable(dfg, sym), fcts; solveKey=solveKey, logger=logger, dens=dens, N=N )

  return mkd, dens, lb, sinfd
end
localProduct(dfg::AbstractDFG, lbl::AbstractString; kw...) = localProduct(dfg, Symbol(lbl); kw...)

"""
    $SIGNATURES

Basic wrapper to take local product and then set the value of `sym` in `dfg`.

Notes
- returns `::Tuple{ManifoldKernelDensity, Float64, Vector{Symbol}}`

DevNotes:
- Unknown issue first occurred here near IIF v0.8.4 tag, recorded case at 2020-01-17T15:26:17.673
"""
function localProductAndUpdate!(dfg::AbstractDFG,
                                sym::Symbol,
                                setkde::Bool=true,
                                logger=ConsoleLogger();
                                solveKey::Symbol=:default )
  #
  # calculate new points for sym using existing structure around sym in dfg
  newPts, dens, lbl, ipc = localProduct(dfg, sym, solveKey=solveKey, N=getSolverParams(dfg).N, logger=logger)
  # maybe update dfg sym with newly calculated points
  setkde && 0 < length(getPoints(newPts)) ? setValKDE!(dfg, sym, newPts, false, ipc, solveKey=solveKey) : nothing

  return newPts, ipc, lbl
end

#


text:
# # Motivation

using InteractiveUtils # hide
using InteractiveUtils: subtypes # hide

# Before going into details about Julia type system, we will spend a few minutes motivating
# the two main roles of the type system, which are:
#
# 1. Structuring the code
# 2. Communicating to the compiler how your type will be used.
#
# The first aspect is important for the convenience of the programmer and enables abstractions
# in the language, the latter aspect is important for the speed of the generated code.
#
# Type systems according to [Wikipedia](https://en.wikipedia.org/wiki/Data_type):
# * In computer science and computer programming, a **data type** or simply **type** is an attribute of data which tells the compiler or interpreter how the programmer intends to use the data.*
# * A **type system** is a logical system comprising a set of rules that assigns a property called a type to the various constructs of a computer program, such as variables, expressions, functions or modules. These types formalize and enforce the otherwise implicit categories the programmer uses for algebraic data types, data structures, or other components.*
#
# ## Structuring the code / enforcing the categories
# The role of **structuring** the code and imposing semantic restriction
# means that the type system allows you to logically divide your program,
# and to prevent certain types of errors.
# Consider for example two types, `Wolf` and `Sheep` which share the same
# definition but the types have different names.

struct Wolf
  name::String
  energy::Int
end

struct Sheep
  name::String
  energy::Int
end

# This allows us to define functions applicable only to the corresponding type
#
howl(wolf::Wolf) = println(wolf.name, " has howled.")
baa(sheep::Sheep) = println(sheep.name, " has baaed.")
nothing # hide
#
# Therefore the compiler (or interpretter) **enforces** that a wolf can only `howl`
# and never `baa` and vice versa a sheep can only `baa`. In this sense, it ensures
# that `howl(sheep)` and `baa(wolf)` never happen.
# For comparison, consider an alternative definition as follows
#
howl(animal) = println(animal.name, " has howled.")
baa(animal)  = println(animal.name, " has baaed.")
#
# in which case the burden of ensuring that a wolf will never baa rests upon the
# programmer which inevitably leads to errors (note that severely constrained
# type systems are difficult to use).

# ## Intention of use and restrictions on compilers
# The *intention of use* in types is related to how efficient code a compiler can
# produce for that given intention. As an example, consider a following two
# alternatives to represent a set of animals:
a = [Wolf("1", 1), Wolf("2", 2), Sheep("3", 3)]
b = (Wolf("1", 1), Wolf("2", 2), Sheep("3", 3))
# where `a` is an array which can contain arbitrary types and have arbitrary length
# whereas `b` is a `Tuple` which has fixed length in which the first two items are of type `Wolf`
# and the third item is of type `Sheep`. Moreover, consider a function which calculates the
# energy of all animals as
energy(animals) = mapreduce(x -> x.energy, +, animals)
nothing # hide
# A good compiler makes use of the information provided by the type system to generate effiecint code
# which we can verify by inspecting the compiled code using `@code_native` macro
@code_native energy(a)
#
@code_native energy(b)
# one observes the second version produces more optimal code. Why is that?
# * In the first representation, `a`, the animals are stored in an `Array` which can have arbitrary size and can contain arbitrary animals. This means that the compiler has to compile `energy(a)` such that it works on such arrays.
# * In the second representation, `b`, the animals are stored in a `Tuple`, which specializes for lengths and types of items. This means that the compiler knows the number of animals and the type of each animal on each position within the tuple, which allows it to specialize.

# This difference will indeed have an impact on the time of code execution.
# On my i5-8279U CPU, the difference (as measured by BenchmarkTools) is
using BenchmarkTools
@benchmark energy(a)
@benchmark energy(b)
# Which nicely demonstrates that the choice of types affects performance. Does it mean that we should always use `Tuples` instead of `Arrays`? Surely not, it is  just that each is better for different use-cases. Using Tuples means that the compiler will compile a special function for each length of tuple and each combination of types of items it contains, which is clearly wasteful.

# # Julia's type system

# ## Julia is dynamicaly typed
# Julia's type system is dynamic, which means that all types are resolved during runtime. **But**, if the compiler can infer types of all variables of the called function, it can specialize the function for that given type of variables which leads to efficient code. Consider a modified example where we represent two wolfpacks:
wolfpack_a =  [Wolf("1", 1), Wolf("2", 2), Wolf("3", 3)]
wolfpack_b =  Any[Wolf("1", 1), Wolf("2", 2), Wolf("3", 3)]
# `wolfpack_a` carries a type `Vector{Wolf}` while `wolfpack_b` has the type `Vector{Any}`. This means that in the first case, the compiler knows that all items are of the type `Wolf`and it can specialize functions using this information. In case of `wolfpack_b`, it does not know which animal it will encounter (although all are of the same type), and therefore it needs to dynamically resolve the type of each item upon its use. This ultimately leads to less performant code.
@btime energy(wolfpack_a)
@btime energy(wolfpack_b)
nothing # hide

# To conclude, julia is indeed a dynamically typed language, **but** if the compiler can infer
# all types in a called function in advance, it does not have to perform the type resolution
# during execution, which produces performant code.

# ## Classes of types
# Julia divides types into three classes: primitive, composite, and abstract.

# ### Primitive types
# Citing the [documentation](https://docs.julialang.org/en/v1/manual/types/#Primitive-Types):  *A primitive type is a concrete type whose data consists of plain old bits. Classic examples  of primitive types are integers and floating-point values. Unlike most languages, Julia lets you declare your own primitive types, rather than providing only a fixed set of built-in ones. In fact, the standard primitive types are all defined in the language itself.*

# The definition of primitive types look as follows
# ```julia
# primitive type Float16 <: AbstractFloat 16 end
# primitive type Float32 <: AbstractFloat 32 end
# primitive type Float64 <: AbstractFloat 64 end
# ```
# and they are mainly used to jump-start julia's type system. It is rarely needed to
# define a special primitive type, as it makes sense only if you define special functions
# operating on its bits. This is almost excusively used for exposing special operations
# provided by the underlying CPU / LLVM compiler. For example `+` for `Int32` is different
# from `+` for `Float32` as they call a different intrinsic operations. You can inspect this
# jump-starting of the type system yourself by looking at Julia's source.
# ```julia
# julia> @which +(1,2)
# +(x::T, y::T) where T<:Union{Int128, Int16, Int32, Int64, Int8, UInt128, UInt16, UInt32, UInt64, UInt8} in Base at int.jl:87
# ```

# At `int.jl:87`
# ```julia
# (+)(x::T, y::T) where {T<:BitInteger} = add_int(x, y)
# ```
# we see that `+` of integers is calling the function `add_int(x, y)`, which is defined in the core
# part of the compiler in `Intrinsics.cpp` (yes, in C++).

# From Julia docs: *Core is the module that contains all identifiers considered "built in" to
# the language, i.e. part of the core language and not libraries. Every module implicitly
# specifies using Core, since you can't do anything without those definitions.*

# Primitive types are rarely used, and they will not be used in this course. We mention them
# for the sake of completeness and refer the reader to the official Documentation (and source code
# of Julia).

# ### Abstract Type
#
# An abstract type can be viewed as a set of concrete types. For example, an
# `AbstractFloat` represents the set of concrete types `(BigFloat,Float64,Float32,Float16)`.
# This is used mainly to define general methods for sets of types for which we expect the same behavior (recall the Julia design motivation: *if it quacks like a duck, waddles like a duck and looks like a duck, chances are it's a duck*). Abstract types are defined with `abstract type TypeName end`. For example the following set of abstract types defines part of julia's number system.
# ```julia
# abstract type Number end
# abstract type Real          <: Number end
# abstract type Complex       <: Number end
# abstract type AbstractFloat <: Real end
# abstract type Integer       <: Real end
# abstract type Signed        <: Integer end
# abstract type Unsigned      <: Integer end
# ```
# where `<:` means "is a subtype of" and it is used in declarations where the right-hand is an immediate sypertype of a given type (`Integer` has the immediate supertype `Real`.) If the supertype is not supplied, it is considered to be Any, therefore in the above defition `Number` has the supertype `Any`. Children of a particular type can be viewed as
using AbstractTrees
function AbstractTrees.children(t::Type)
    t === Function ? Vector{Type}() : filter!(x -> x !== Any,subtypes(t))
end
AbstractTrees.printnode(io::IO,t::Type) = print(io,t)
print_tree(Number)
# As was mentioned, abstract types allows as to define functions that can be applied variebles of types with a given abstract type as a supertype. For example we can define a `sgn` function for **all** real numbers as
sgn(x::Real) = x > 0 ? 1 : x < 0 ? -1 : 0
nothing # hide
# and we know it would be correct for all real numbers. This means that if anyone creates
# a new subtype of `Real`, the above function can be used. This also means that
# **it is expected** that comparison operations are defined for any real number. Also notice that
# `Complex` numbers are excluded, since they do not have a total order.

# For unsigned numbers, the `sgn` can be simplified, as it is sufficient to verify if they are different (greater) than zero, therefore the function can read
sgn(x::Unsigned) = x > 0 ? 1 : 0
nothing # hide
# and again, it applies to all numbers derived from `Unsigned`. Recall that
# `Unsigned <: Integer <: Real,` how does Julia decide,
# which version of the function `sgn` to use for `UInt8(0)`? It chooses the most
# specific version, and thus for `sgn(UInt8(0))` it will use `sgn(x::Unsinged)`.
# If the compiler cannot decide, typically it encounters an ambiguity, it throws an error
# and recommends which function you should define to resolve it.

# The above behavior allows to define default "fallback" implementations and while allowing
# to specialize for sub-types. A great example is matrix multiplication, which has a
# generic (and slow) implementation with many specializations, which can take advantage
# of structure (sparse, banded), or use optimized implementations (e.g. blas implementation
# for dense matrices with eltype `Float32` and `Float64`).

# Again, Julia does not make a difference between abstract types defined in `Base`
# libraries shipped with the language and those defined by you (the user). All are treated
# the same.

# (![From Julia documentation](https://docs.julialang.org/en/v1/manual/types/#man-abstract-types))
# Abstract types cannot be instantiated, which means that we cannot create a variable that
# would have an abstract type (try `typeof(Number(1f0))`). Also, abstract types cannot have
# any fields, therefore there is no composition (there are lengty discussions of why this is so,
# one of the most definite arguments of creators is that abstract types with fields frequently lead
# to children types not using some fields (consider circle vs. ellipse)).

# ### [Composite types](@id composite_types)
# Composite types are similar to `struct` in C (they even have the same memory layout) as they logically join together other types. It is not a great idea to think about them as objects (in OOP sense), because objects tie together *data* and *functions* on owned data. Contrary in Julia (as in C), functions operate on data of structures, but are not tied to them and they are defined outside them. Composite types are workhorses of Julia's type system, as user-defined types are mostly composite (or abstract).

# Composite types are defined using `struct TypeName [fields] end`. To define a position of an animal on the Euclidean plane as a type, we would write
struct PositionF64
  x::Float64
  y::Float64
end
# which defines a structure with two fields `x` and `y` of type `Float64`. Julia's compiler creates a default constructor, where both (but generally all) arguments are converted using `(convert(Float64, x), convert(Float64, y)` to the correct type. This means that we can construct a PositionF64 with numbers of different type that are convertable to Float64, e.g. `PositionF64(1,1//2)` but we cannot construct `PositionF64` where the fields would be of different type (e.g. `Int`, `Float32`, etc.) or they are not trivially convertable (e.g. `String`).

# Fields in composite types do not have to have a specified type.  We can define a `VaguePosition` without specifying the type
struct VaguePosition
  x
  y
end
# This works as the definition above except that the arguments are not converted to `Float64` now. One can store different values in `x` and `y`, for example `String` (e.g. VaguePosition("Hello","world")). Although the above definition might be convenient, it limits the compiler's ability to specialize, as the type  `VaguePosition` does not carry information about type of `x` and `y`, which has a negative impact on the performance. For example
using BenchmarkTools
move(a::T,b::T) where T = T(a.x + b.x, a.y + b.y)
x = [PositionF64(rand(), rand()) for _ in 1:100]
y = [VaguePosition(rand(), rand()) for _ in 1:100]
@benchmark reduce(move, x)
@benchmark reduce(move, y)
# Giving fields of a composite type an abstract type does not really solve the problem of the compiler not knowing the type. In this example, it still does not know, if it should use instructions for `Float64` or `Int8`.
struct LessVaguePosition
  x::Real
  y::Real
end
z = [LessVaguePosition(rand(), rand()) for _ in 1:100];
@benchmark reduce(move, z)
# From the perspective of generating optimal code, both definitions are equally uninformative to the compiler as it cannot assume anything about the code. However, the  `LessVaguePosition` will ensure that the position will contain only numbers, hence catching trivial errors like instantiating `VaguePosition` with non-numeric types for which arithmetic operators will not be defined (recall the discussion on the  beginning of the lecture).

# All structs defined above are immutable (as we have seen above in the case of `Tuple`), which means that one cannot change a field (unless the struct wraps a container, like and array, which allows that). For example this raises an error
a = LessVaguePosition(1,2)
a.x = 2

# If one needs to make a struct mutable, use the keyword `mutable` before the keyword `struct` as
mutable struct MutablePosition
  x::Float64
  y::Float64
end
# In mutable structures, we can change the values of fields.
a = MutablePosition(1e0, 2e0)
a.x = 2
a
# Note, that the memory layout of mutable structures is different, as fields now contain references to memory locations, where the actual values are stored.

# ### Parametric Types
# So far, we had to trade-off flexibility for generality in type definitions. Can we have both? The answer is affirmative. The way to achieve this  **flexibility** in definitions of the type while being  able to generate optimal code is to  **parametrize** the type definition. This is achieved by replacing types with a parameter (typically a single uppercase character) and decorating in definition by specifying different type in curly brackets. For example
struct PositionT{T}
  x::T
  y::T
end
u = [PositionT(rand(), rand()) for _ in 1:100];
@btime reduce(move, u);
# Notice that the compiler can take advantage of specializing for different types (which does not have an effect here as in modern processors addition of `Float` and `Int` takes the same time).
v = [PositionT(rand(1:100), rand(1:100)) for _ in 1:100];
@btime reduce(move, v);
# The above definition suffers the same problem as `VaguePosition`, which is that it allows us to instantiate the `PositionT` with non-numeric types, e.g. `String`. We solve this by restricting the types `T` to be children of some supertype, in this case `Real`
struct Position{T<:Real}
  x::T
  y::T
end
# which will throw an error if we try to initialize it with `Position("1.0", "2.0")`.
#
# Naturally, fields in structures can be of different types, as is in the below pointless example.
struct PositionXY{X<:Real, Y<:Real}
  x::X
  y::Y
end

# ### Abstract parametric types
# Like Composite types, Abstract types can also have parameters. These parameters define types that are common for all child types. A very good example is Julia's definition of arrays of arbitrary dimension `N` and type `T` of its items as
# ```julia
# abstract type AbstractArray{T,N} end
# ```
# Different `T` and `N` give rise to different variants of `AbstractArrays`,
# therefore `AbstractArray{Float32,2}` is different from `AbstractArray{Float64,2}`
# and from `AbstractArray{Float64,1}.` Note that these are still `Abstract` types,
# which means you cannot instantiate them. Their purpose is
# * to allow to define operations for broad class of concrete types
# * to inform the compiler about constant values, which can be used
# Notice in the above example that parameters of types do not have to be types, but can also be values of primitive types, as in the above example of `AbstractArray` `N` is the number of dimensions which is an integer value.
#
# For convenience, it is common to give some important partially instantiated Abstract types an **alias**, for example `AbstractVector` as
# ```julia
# const AbstractVector{T} = AbstractArray{T,1}
# ```
# is defined in `array.jl:23` (in Julia 1.6.2), which allows us to define for example general prescription for the `dot` product of two abstract vectors as
function dot(a::AbstractVector, b::AbstractVector)
  @assert length(a) == length(b)
  mapreduce(*, +, a, b)
end
# You can verify that the above general function can be compiled to performant code if
# specialized for particular arguments.
@code_native mapreduce(*,+, [1,2,3], [1,2,3])

# ## More on the use of types in function definitions
# ### Terminology
# * A *function* refers to a set of "methods" for a different combination of type parameters (the term function can be therefore considered as refering to a mere **name**). *Methods* define different behavior for different types of arguments for a given function. For example
move(a::Position, b::Position) = Position(a.x + b.x, a.y + b.y)
move(a::Vector{<:Position}, b::Vector{<:Position}) = move.(a,b)
# `move` refers to a function with methods `move(a::Position, b::Position)` and `move(a::Vector{<:Position}, b::Vector{<:Position})`. When different behavior on different types is defined by a programmer, as shown above, it is also called *implementation specialization*. There is another type of specialization, called *compiler specialization*, which occurs when the compiler generates different functions for you from a single method. For example for
move(Position(1,1), Position(2,2))
move(Position(1.0,1.0), Position(2.0,2.0))
# the compiler generates two methods, one for `Position{Int64}` and the other for `Position{Float64}`. Notice that inside generated functions, the compiler needs to use different intrinsic operations, which can be viewed from
@code_native move(Position(1,1), Position(2,2))
# and
@code_native move(Position(1.0,1.0), Position(2.0,2.0))

# ## Intermezzo: How does the Julia compiler work?
# Let's walk through an example. Consider the following definitions
move(a::Position, by::Position) = Position(a.x + by.x, a.y + by.y)
move(a::T, by::T) where {T<:Position} = Position(a.x + by.x, a.y + by.y)
move(a::Position{Float64}, by::Position{Float64}) = Position(a.x + by.x, a.y + by.y)
move(a::Vector{<:Position}, by::Vector{<:Position}) = move.(a, by)
move(a::Vector{<:Position}, by::Position) = move.(a, by)
nothing # hide
# and a function call
a = Position(1.0, 1.0)
by = Position(2.0, 2.0)
move(a, by)
# 1. The compiler knows that you call the function `move`.
# 2. The compiler infers the type of the arguments. You can view the result with
(typeof(a),typeof(by))
# 3. The compiler identifies all `move`-methods with arguments of type `(Position{Float64}, Position{Float64})`:
Base.method_instances(move, (typeof(a), typeof(by)))
# (For demonstration we can store the method we want in a variable `m`)
m = Base.method_instances(move, (typeof(a), typeof(by))) |> first
nothing # hide
# 4a. If the method has been specialized (compiled), then the arguments are prepared and the method is invoked. The compiled specialization can be seen from
m.cache
# 4b. If the method has not been specialized (compiled), the method is compiled for the given type of arguments and continues as in step 4a.
# A compiled function is therefore  a "blob" of **native code** living in a particular memory location. When Julia calls a function, it needs to pick the right block corresponding to a function with particular type of parameters.
#
# If the compiler cannot narrow the types of arguments to concrete types, it has to perform the above procedure inside the called function, which has negative effects on performance, as the type resulution and identification of the methods can be slow, especially for methods with many arguments (e.g. 30ns for a method with one argument,
# 100 ns for method with two arguements).
# Recall the above example
wolfpack_a =  [Wolf("1", 1), Wolf("2", 2), Wolf("3", 3)]
@benchmark energy(wolfpack_a)
# and
wolfpack_b =  Any[Wolf("1", 1), Wolf("2", 2), Wolf("3", 3)]
@benchmark energy(wolfpack_b)
# An interesting intermediate between fully abstract and fully concrete type happens, when the compiler knows that arguments have abstract type, which is composed of a small number of concrete types. This case  called Union-Splitting, which happens when there is just a little bit of uncertainty. Julia will do something like
# ```julia
# argtypes = typeof(args)
# push!(execution_stack, args)
# if T == Tuple{Int, Bool}
#   @goto compiled_blob_1234
# else # the only other option is Tuple{Float64, Bool}
#   @goto compiled_blob_1236
# end
# ```
# For example
const WolfOrSheep = Union{Wolf, Sheep}
wolfpack_c =  WolfOrSheep[Wolf("1", 1), Wolf("2", 2), Wolf("3", 3)]
@benchmark energy(wolfpack_c)
# Thanks to union splitting, Julia is able to have performant operations on arrays with undefined / missing values for example
[1, 2, 3, missing] |> typeof

# ### More on matching methods and arguments
# In the above process, the step, where Julia looks for a method instance with corresponding parameters can be very confusing. The rest of this lecture will focus on this. For those who want to have a formal background, we recommend [talk of  Francesco Zappa Nardelli](https://www.youtube.com/watch?v=Y95fAipREHQ) and / or the one of [Jan Vitek](https://www.youtube.com/watch?v=LT4AP7CUMAw).
#
# When Julia needs to specialize a method instance, it needs to find it among multiple definitions. A single function can have many method instances, see for example `methods(+)` which  lists all method instances of the `+`-function. How does Julia select the proper one?
# 1. It finds all methods where the type of arguments match or are subtypes of restrictions on arguments in the method definition.
# 2a. If there are multiple matches, the compiler selects the most specific definition.
# 2b. If the compiler cannot decide, which method instance to choose, it throws an error.
confused_move(a::Position{Float64}, by) = Position(a.x + by.x, a.y + by.y)
confused_move(a, by::Position{Float64}) = Position(a.x + by.x, a.y + by.y)
confused_move(Position(1.0,2.0), Position(1.0,2.0))
# 2c. If it cannot find a suitable method, it throws an error.
move(Position(1,2), VaguePosition("hello","world"))

# Some examples
# Consider following definitions
move(a::Position, by::Position) = Position(a.x + by.x, a.y + by.y)
move(a::T, by::T) where {T<:Position} = T(a.x + by.x, a.y + by.y)
move(a::Position{Float64}, by::Position{Float64}) = Position(a.x + by.x, a.y + by.y)
move(a::Vector{<:Position}, by::Vector{<:Position}) = move.(a, by)
move(a::Vector{T}, by::Vector{T}) where {T<:Position} = move.(a, by)
move(a::Vector{<:Position}, by::Position) = move.(a, by)
# Which method will compiler select for
move(Position(1.0,2.0), Position(1.0,2.0))
# The first three methods match the types of argumens, but the compiler will select the third one, since it is the most specific.
#
# Which method will compiler select for
move(Position(1,2), Position(1,2))
# Again, the first and second method definitions match the argument, but the second is the most specific.
#
# Which method will the compiler select for
move([Position(1,2)], [Position(1,2)])
# Again, the fourth and fifth method definitions match the argument, but the fifth is the most specific.
move([Position(1,2), Position(1.0,2.0)], [Position(1,2), Position(1.0,2.0)])

# ### Frequent problems
# 1. Why does the following fail?
foo(a::Vector{Real}) = println("Vector{Real}")
foo([1.0,2,3])
# Julia's type system is **invariant**, which means that `Vector{Real}` is different from `Vector{Float64}` and from `Vector{Float32}`, even though `Float64` and `Float32` are sub-types of `Real`. Therefore `typeof([1.0,2,3])` isa `Vector{Float64}` which is not subtype of `Vector{Real}.` For **covariant** languages, this would be true. For more information on variance in computer languages, [see here](https://en.wikipedia.org/wiki/Covariance_and_contravariance_(computer_science)). If de above definition of `foo` should be applicable to all vectors which has elements of subtype of `Real` we have define it as
foo(a::Vector{T}) where {T<:Real} = println("Vector{T} where {T<:Real}")
# or equivalently but more tersely as
foo(a::Vector{<:Real}) = println("Vector{T} where {T<:Real}")
# 2. Diagonal rule says that a repeated type in a method signature has to be a concrete type. Consider for example the function below
move(a::T, b::T) where {T<:Position}
# we cannot call it with `move(Position(1.0,2.0), Position(1,2))`, since in this case `Position(1.0,2.0)` is of type `Position{Float64}` while `Position(1,2)` is of type `Position{Int64}`.
# 3. When debugging why arguments do not match a particular method definition, it is useful to use `typeof`, `isa`, and `<:` commands. For example
typeof(Position(1.0,2.0))
#
typeof(Position(1,2))
#
Position(1,2) isa Position{Float64}
#
Position(1,2) isa Position{Real}
#
Position(1,2) isa Position{<:Real}
#
typeof(Position(1,2)) <: Position{<:Float64}
#
typeof(Position(1,2)) <: Position{<:Real}

# ### A bizzare definition which you can encounter
# The following definition of a one-hot matrix is taken from [Flux.jl](https://github.com/FluxML/Flux.jl/blob/1a0b51938b9a3d679c6950eece214cd18108395f/src/onehot.jl#L10-L12)
struct OneHotArray{T<:Integer, L, N, var"N+1", I<:Union{T,AbstractArray{T, N}}} <: AbstractArray{Bool, var"N+1"}
  indices::I
end
# The parameters of the type carry information about the type used to encode the position of `one` in each column in `T`, the dimension of one-hot vectors in `L`, the dimension of the storage of `indices` in `N` (which is zero for `OneHotVector` and one for `OneHotMatrix`), number of dimensions of the `OneHotArray` in `var"N+1"` and the type of underlying storage of indicies `I`.


text:
# Implements the Hillis--Steele algorithm using global memory
# See algorithm 1 at https://en.wikipedia.org/wiki/Prefix_sum#Parallel_algorithm

# TODO: features
# - init::Some
# - CuMatrix
# - pairwise

# TODO: performance
# - shared memory / shuffle (see CUDAnative.jl/examples/scan)

function Base._accumulate!(op::Function, vout::CuVector{T}, v::CuVector, dims::Int,
                           init::Nothing) where {T}
    if dims != 1
        return copyto!(vout, v)
    end

    return Base._accumulate!(op::Function, vout::CuVector{T}, v::CuVector, nothing, nothing)
end

function Base._accumulate!(op::Function, vout::CuVector{T}, v::CuVector, dims::Nothing,
                           init::Nothing) where {T}
    vin = T.(v)  # convert to vector with eltype T

    Δ = 1   # Δ = 2^d
    n = ceil(Int, log2(length(v)))

    num_threads = 256
    num_blocks = ceil(Int, length(v) / num_threads)

    for d in 0:n   # passes through data
        @cuda blocks=num_blocks threads=num_threads _partial_accumulate!(op, vout, vin, Δ)

        vin, vout = vout, vin
        Δ *= 2
    end

    return vin
end

function _partial_accumulate!(op, vout, vin, Δ)
    @inbounds begin
        k = threadIdx().x + (blockIdx().x - 1) * blockDim().x

        if k <= length(vin)
            if k > Δ
                vout[k] = op(vin[k - Δ], vin[k])
            else
                vout[k] = vin[k]
            end
        end
    end

    return
end

Base.accumulate_pairwise!(op, result::CuVector, v::CuVector) = accumulate!(op, result, v)


text:
# ---------------------------------------------------------------------------- #
#                                     REGEX                                    #
# ---------------------------------------------------------------------------- #
# ---------------------------------- markup ---------------------------------- #
const OPEN_TAG_REGEX = r"\[[a-zA-Z _0-9. ,()#]+[^/\[]\]"
const CLOSE_TAG_REGEX = r"\[\/[a-zA-Z _0-9. ,()#]+[^/\[]\]"
const GENERIC_CLOSER_REGEX = r"\[\/\]"

"""
    remove_markup(input_text::AbstractString)::AbstractString

Remove all markup tags from a string of text.
"""
function remove_markup(input_text::AbstractString)::AbstractString
    text = replace_double_brackets(input_text)

    # remove extra closing tags
    text = remove_markup_open(text)
    text = replace(text, GENERIC_CLOSER_REGEX => "")
    text = replace(text, CLOSE_TAG_REGEX => "")

    return reinsert_double_brackets(text)
end

"""
    remove_markup_open(text::AbstractString)

Remove all opening markup tags from a string of text
"""
remove_markup_open(text::AbstractString)::AbstractString =
    replace(text, OPEN_TAG_REGEX => "")

# ----------------------------------- ansi ----------------------------------- #
const ANSI_REGEXEs = [r"\e\[[0-9]*m", r"\e\[[0-9;]*m"]

"""
    remove_ansi(str::AbstractString)::AbstractString

Remove all ANSI tags from a string of text
"""
function remove_ansi(str::AbstractString)::AbstractString
    for regex in ANSI_REGEXEs
        str = replace(str, regex => "")
    end
    return str
end

"""
    cleantext(str::AbstractString)

Remove all style information from a string.
"""
cleantext(str::AbstractString) = (remove_ansi ∘ remove_markup_open ∘ remove_markup)(str)

# --------------------------------- brackets --------------------------------- #
const brackets_regexes = [(r"\[", "[["), (r"\]", "]]")]

"""
    remove_ansi(str::AbstractString)::AbstractString

Replace each squared bracket with a double copy of itself
"""
function escape_brackets(text::AbstractString)::AbstractString
    for (regex, replacement) in brackets_regexes
        text = replace(text, regex => replacement)
    end
    return text
end

const remove_brackets_regexes = [(r"\[\[", "["), (r"\]\]", "]")]

"""
    unescape_brackets(text::AbstractString)::AbstractString

Replece every double squared parenthesis with a single copy of itself
"""
function unescape_brackets(text::AbstractString)::AbstractString
    for (regex, replacement) in remove_brackets_regexes
        text = replace(text, regex => replacement)
    end
    return text
end

"""
    replace_double_brackets(text::AbstractString)::AbstractString

Replace double brackets with %% and ±± to avoid them being picked up by markup extraction
"""
function replace_double_brackets(text::AbstractString)::AbstractString
    text = replace(text, "[[" => "%%")
    text = replace(text, "]]" => "±±")
    return text
end

"""
    reinsert_double_brackets(text::AbstractString)::AbstractString

Insert previously replaced double brackets
"""
function reinsert_double_brackets(text::AbstractString)::AbstractString
    text = replace(text, "%%" => "[[")
    text = replace(text, "±±" => "]]")
    return text
end

# ---------------------------------------------------------------------------- #
#                                      I/O                                     #
# ---------------------------------------------------------------------------- #
"""
    read_file_lines(path::String, start::Int, stop::Int)

Read a file and select only lines in range `start` -> `stop`.

Returns a vector of tuples with the line number and line content.
"""
function read_file_lines(path::AbstractString, start::Int, stop::Int)
    start = start < 1 ? 1 : start
    stop = stop >= countlines(path) ? countlines(path) : stop
    lines = readlines(path; keep = true)
    return collect(enumerate(lines))[start:stop]
end

# ---------------------------------------------------------------------------- #
#                                     MISC                                     #
# ---------------------------------------------------------------------------- #

"""
    replace_text(text::AbstractString, start::Int, stop::Int, replace::AbstractString)

Replace a section of a `text` between `start` and `stop` with `replace`.
"""
function replace_text(text::AbstractString, start::Int, stop::Int, replace::AbstractString)
    if start == 1
        return replace * text[stop:end]
    elseif stop == length(text)
        return text[1:start] * replace
    else
        return text[1:start] * replace * text[stop:end]
    end
end

"""
    replace_text(text::AbstractString, start::Int, stop::Int, char::Char='_')

Replace a section of a `text`  between `start` and `stop` with another string composed of repeats of a given character `char`.
"""
function replace_text(text::AbstractString, start::Int, stop::Int, char::Char = '_')
    replacement = char^(stop - start - 1)
    return replace_text(text, start, stop, replacement)
end

"""
    nospaces(text::AbstractString)

Remove all spaces from a string.
"""
nospaces(text::AbstractString) = replace(text, " " => "")

"""
    remove_brackets(text::AbstractString)

Remove all () brackets from a string.
"""
remove_brackets(text::AbstractString) = replace(replace(text, "(" => ""), ")" => "")

"""
    square_to_round_brackets(text::AbstractString)

Replace square brackets with round ones.
"""
function square_to_round_brackets(text::AbstractString)
    return replace(replace(text, "[" => "("), "]" => ")")
end

"""
    unspace_commas(text::AbstractString)

Remove spaces after commas.
"""
unspace_commas(text::AbstractString) = replace(replace(text, ", " => ","), ". " => ".")

"""
Split a string into a vector of Chars.
"""
chars(text::AbstractString)::Vector{Char} = [x for x in text]

"""
    join_lines(lines)

Merge a vector of strings in a single string.
"""
join_lines(lines) = join(lines, "\n")

"""
    split_lines(text::AbstractString)

Split a string into its composing lines.
"""
function split_lines(text::AbstractString)
    return split(text, "\n")
end

"""
    split_lines(renderable)

Split a renderable's text.
"""
function split_lines(renderable)
    if string(typeof(renderable)) == "Segment"
        return split_lines(renderable.text)
    else
        [s.text for s in renderable.segments]
    end
end

# ------------------------------- reshape text ------------------------------- #
"""
    fillin(text::AbstractString)::AbstractString

Ensure that each line in a multi-line text has the same width.
"""
function fillin(text::AbstractString)::AbstractString
    lines = split_lines(text)
    length(lines) == 1 && return text

    widths = map(textlen, lines)
    w = max(widths...)
    filled::Vector{AbstractString} = []
    for (i, ln) in enumerate(lines)
        push!(filled, ln * " "^(w - widths[i]))
    end
    return join_lines(filled)
end

"""
    truncate(text::AbstractString, width::Int)

Shorten a string of text to a target width
"""
function truncate(text::AbstractString, width::Int)
    if textlen(text) <= width
        return text
    end

    cut = get_last_valid_str_idx(text, width - 3)
    return text[1:cut] * "..."
end

"""
    get_valid_chars!(valid_chars::Vector{Int}, tag, δ::Int)

Recursively extract valid characters (i.e. not in markup tags) from a string.
"""
function get_valid_chars!(valid_chars::Vector{Int}, tag, δ::Int)
    # get correct start/stop positions
    s1, e1 = δ + tag.open.start, δ + tag.open.stop
    s2, e2 = δ + tag.close.start, δ + tag.close.stop

    # do nested tags
    for inner in tag.inner_tags
        get_valid_chars!(valid_chars, inner, e1)
    end

    if s2 > length(valid_chars)
        @debug "How can tag close after valid chars?" tag tag.open tag.close length(
            valid_chars
        ) tag.text
    else
        valid_chars[s1:e1] .= 0
        valid_chars[s2:e2] .= 0
    end

    return valid_chars
end

"""
    textlen(x::AbstractString)

Get length of text after all style information is removed.
"""
textlen(x::AbstractString) = (textwidth ∘ remove_markup ∘ remove_markup_open ∘ remove_ansi)(x)

"""
    reshape_text(text::AbstractString, width::Int)

Reshape `text` to have a given `width`.

When `text` is longer than `width`, it gets cut into multiple lines.
This is done carefully to preserve style information by: avoiding
cutting inside style markup and copying markup tags over to new lines
so that the style is correctly applied.
"""
function reshape_text(text::AbstractString, width::Int)::AbstractString

    # check if no work is required
    if textlen(text) <= width
        return text
    end

    # extract tag and mark valid characters and "cutting" places
    tags = extract_markup(text)
    valid_chars = ones(Int, length(text))
    spaces = isspace.(chars(text))
    for tag in tags
        get_valid_chars!(valid_chars, tag, 0)
    end

    # create lines with splitted tex
    lines::Vector{AbstractString} = []
    j = 1
    while textlen(text) > width
        # get a cutting index not in a tag's markup
        condition = (cumsum(valid_chars[j:end]) .<= width) .& valid_chars[j:end] .== 1 .& spaces[j:end] .== 1
        cut = findlast(condition)

        if isnothing(cut)
            # couldnt find a cut in the spaces, try without
            condition = (cumsum(valid_chars[j:end]) .<= width) .& valid_chars[j:end] .== 1
            cut = findlast(condition)
        end

        if cut + j > length(valid_chars)
            @warn "ops not valid" j cut length(valid_chars)
            cut = length(valid_chars) - j
        end

        # prep line
        try
            newline = strip(text[1:cut])

            # pad new line with spaces to ensure it has the right lengt
            newline *= " ".^(width - textlen(newline))

            push!(lines, newline)
            text = lstrip(text[(cut + 1):end])
            j += cut
        catch err
            throw("Failed to reshape text: $err - target width: $width")
        end
    end

    # add what's left of the text
    if length(text) > 0
        push!(lines, text * " "^(width - textlen(text)))
    end

    # do checks and pad line
    for (n, line) in enumerate(lines)
        h = remove_markup(line)
        # @assert length(remove_markup(line)) <= width

        ll = length(remove_markup(line))
        if ll < width
            lines[n] = line * " "^(width - ll)
        end
    end

    return join_lines(pairup_tags(lines))
end

# ------------------------------------ end ----------------------------------- #

"""
    do_by_line(fn::Function, text::AbstractString)

Apply `fn` to each line in the `text`.

The function `fn` should accept a single `::AbstractString` argument.
"""
function do_by_line(fn::Function, text::AbstractString)
    lines::Vector{AbstractString} = []
    for line in split_lines(text)
        push!(lines, fn(line))
    end
    return join_lines(lines)
end

do_by_line(fn, text::Vector) = do_by_line(fn, join_lines(text))

"""
    get_last_valid_str_idx(str::AbstractString, idx::Int)

Get valid index to cut a string at.

When indexing a string, the number of indices is given by the
the sum of the `ncodeunits` of each `Char`, but some indices
will not be valid. This function ensures that given a (potentially)
not valid index, the last valid one is elected.
"""
function get_last_valid_str_idx(str::AbstractString, idx::Int)
    while !isvalid(str, idx)
        idx -= 1

        if idx <= 0
            throw("Failed to find a valid index for $str starting at $idx")
        end
    end
    return idx
end

function get_last_valid_str_idx(str::AbstractString, idx::Int, valid_places::Vector{Int64})
    while !isvalid(str, idx) || valid_places[idx] == 0
        idx -= 1

        if idx == 0
            break
        end
    end
    return idx
end

"""
get_next_valid_str_idx(str::AbstractString, idx::Int)

Get valid index to cut a string at.

When indexing a string, the number of indices is given by the
the sum of the `ncodeunits` of each `Char`, but some indices
will not be valid. This function ensures that given a (potentially)
not valid index, the next valid one is elected.
"""
function get_next_valid_str_idx(str::AbstractString, idx::Int)
    while !isvalid(str, idx)
        idx += 1

        if idx >= length(str)
            throw("Failed to find a valid index for $str starting at $idx")
        end
    end
    return idx
end