text:
# This file is a part of BATBase.jl, licensed under the MIT License (MIT).

"""
    abstract type AbstractDensity

Subtypes of `AbstractDensity` (e.g. `MyDensity`) must implement the methods

* `BATBase.eval_logval_unchecked(density::MyDensity)`

To query the logarithm of the density for a given point/variate `v`, use
the function [`logdensityof`](@ref).

For likelihood densities, implementing `BATBase.eval_logval_unchecked is
typically sufficient, since shape, and variate bounds can be inferred from
the prior.

Densities with a known variate shape may also implement

* `ValueShapes.varshape(density::MyDensity)`

Densities with known variate bounds may also implement

* `BAT.var_bounds(density::MyDensity)`

!!! note

    Only implement `BAT.var_bounds` in special cases. `BAT.var_bounds` is not
    part of the stable public BATBase.jl and BAT.jl APIs and subject to change
    without deprecation.
"""
abstract type AbstractDensity end
export AbstractDensity

Base.convert(::Type{AbstractDensity}, d::ContinuousDistribution) = convert(AbstractDistributionDensity, d)

"""
    BATBase.eval_logval_unchecked(density::AbstractDensity, v::Any)

Compute log of the value of a multivariate density function for the given
variate/parameter-values.

Input:

* `density`: density function
* `v`: argument, i.e. variate / parameter-values

Note: If `eval_logval_unchecked` is called with an argument that is out of bounds,
the behavior is undefined. The result for arguments that are not within
bounds is *implicitly* `-Inf`, but it is the caller's responsibility to handle
these cases.
"""
function eval_logval_unchecked end

"""
    ValueShapes.totalndof(density::AbstractDensity)::Union{Int,Missing}

Get the number of degrees of freedom of the variates of `density`. May return
`missing`, if the shape of the variates is not fixed.
"""
@inline function ValueShapes.totalndof(density::AbstractDensity)
    shape = varshape(density)
    ismissing(shape) ? missing : ValueShapes.totalndof(shape)
end

"""
    ValueShapes.varshape(
        density::AbstractDensity
    )::Union{ValueShapes.AbstractValueShape,Missing}

    ValueShapes.varshape(
        density::DistLikeDensity
    )::ValueShapes.AbstractValueShape

Get the shapes of the variates of `density`.

For prior densities, the result must not be `missing`, but may be `nothing` if
the prior only supports unshaped variate/parameter vectors.
"""
ValueShapes.varshape(density::AbstractDensity) = missing

"""
    logdensityof(density::AbstractDensity, v)::Real
    logdensityof(density::AbstractDensity)::Function

Computes the logarithmic value of `density` at a given point, resp. returns a
function that does so:

```julia
logdensityof(density, v) == logdensityof(density)(v)
```

The function returned by `logdensityof(density)` supports

```julia
(- logdensityof(density))(v) == - logdensityof(density)(v)
```

Note: This function should *not* be specialized for custom density types,
implement [`BATBase.eval_logval_unchecked`](@ref) instead.

!!! note

    This functionality is implemented in BAT.jl, not in BATBase.jl.
"""
function logdensityof end
export logdensityof

# Implemented in BAT.jl:
# logdensityof(density::AbstractDensity, v::Any) = ...
# logdensityof(density::AbstractDensity) = ...

# ToDo, move to BAT.jl:
#=

logdensityof(density::AbstractDensity, v::Any) = eval_logval(density, v, default_dlt())

logdensityof(density::AbstractDensity) = LogDensityOf(density)

struct LogDensityOf{D<:AbstractDensity} <: Function
    density::D
end

@inline (lvd::LogDensityOf)(v::Any) = logdensityof(lvd.density, v)

(Base.:-)(lvd::LogDensityOf) = NegLogDensityOf(lvd.density)

ValueShapes.varshape(lvd::LogDensityOf) = varshape(lvd.density)
ValueShapes.unshaped(lvd::LogDensityOf) = LogDensityOf(unshaped(lvd.density))

struct NegLogDensityOf{D<:AbstractDensity} <: Function
    density::D
end

@inline (lvd::NegLogDensityOf)(v::Any) = - logdensityof(lvd.density, v)

(Base.:-)(lvd::NegLogDensityOf) = LogDensityOf(lvd.density)

ValueShapes.varshape(lvd::NegLogDensityOf) = varshape(lvd.density)
ValueShapes.unshaped(lvd::NegLogDensityOf) = NegLogDensityOf(unshaped(lvd.density))
=#

"""
    abstract type DistLikeDensity <: AbstractDensity

A density that implements part of the `Distributions.Distribution` interface.
Such densities are suitable for use as a priors.

In contrast to a distribution, the integral of a `DistLikeDensity` is not
required to be normalized to one (but will be, in many cases).

Typically, custom priors should be implemented as subtypes of
`Distributions.Distribution`. BAT will automatically wrap them in a subtype of
[`AbstractDistributionDensity`](@ref), which is a subtype of
`DistLikeDensity`.

Subtypes of `DistLikeDensity` are required to support more functionality
than an [`AbstractDensity`](@ref), but less than a
`Distribution{Multivariate,Continuous}`.

A `d::ContinuousDistribution` can be converted into (wrapped
in) an suitable subtype of [`AbstractDistributionDensity`](@ref) via
`conv(DistLikeDensity, d)`.

The following functions must be implemented for subtypes (e.g.
`MyDistLikeDensity`), in addition to the mandatory [`AbstractDensity`](@ref)
interface:

* `ValueShapes.varshape(density::MyDistLikeDensity)`

* `Distributions.sampler(density::MyDistLikeDensity)`

!!! note

    The function `BAT.var_bounds` is not part of the stable public BATBase.jl
    and BAT.jl APIs and subject to change without deprecation.
"""
abstract type DistLikeDensity <: AbstractDensity end
export DistLikeDensity

Base.convert(::Type{DistLikeDensity}, d::ContinuousDistribution) = convert(AbstractDistributionDensity, d)

"""
    abstract type AbstractDistributionDensity <: DistLikeDensity

A density that can be converted efficiently to a `Distributions.Distribution`.

Subtypes (e.g. `MyDistributionDensity`) support, in addition to the
[`DistLikeDensity`](@ref) interface, support:

* `Base.convert(Distributions.Distribution, density::MyDistributionDensity)`

A `d::ContinuousDistribution` can be converted into (resp. wrapped in) an
suitable subtype of `DistLikeDensity` via `conv(AbstractDistributionDensity, d)`.

!!! note

    This functionality is implemented in BAT.jl, not in BATBase.jl.
"""
abstract type AbstractDistributionDensity <: DistLikeDensity end
export AbstractDistributionDensity

# Implemented in BAT.jl:
# Base.convert(::Type{AbstractDistributionDensity}, d::ContinuousDistribution) = ...

"""
    BAT.AnyDensityLike = Union{...}

Union of all types that BAT will accept as a probability density, resp. that
`convert(AbstractDensity, d)` supports:

* [`AbstractDensity`](@ref)
* `Distributions.Distribution`
"""
const AnyDensityLike = Union{
    AbstractDensity,
    Distributions.ContinuousDistribution
}
export AnyDensityLike

"""
    BAT.AnySampleable = Union{...}

Union of all types that BAT can sample from:

* [`AbstractDensity`](@ref)
* [`AbstractDensitySampleVector`](@ref)
* `Distributions.Distribution`
"""
const AnySampleable = Union{
    AbstractDensity,
    AbstractDensitySampleVector,
    Distributions.Distribution
}
export AnySampleable

"""
    BAT.AnyIIDSampleable = Union{...}

Union of all distribution/density-like types that BAT can draw i.i.d.
(independent and identically distributed) samples from:

* [`DistLikeDensity`](@ref)
* `Distributions.Distribution`
"""
const AnyIIDSampleable = Union{
    DistLikeDensity,
    Distributions.Distribution
}
export AnyIIDSampleable


text:
type FiveNearUniformArms <: BanditEnvironment
end

function arms(env::FiveNearUniformArms, trial::Int64)
  [
    BernoulliArm(0.48),
    BernoulliArm(0.49),
    BernoulliArm(0.51),
    BernoulliArm(0.50),
    BernoulliArm(0.52)
   ]
end

function n_arms(env::FiveNearUniformArms)
  5
end

ext: jl

max_stars_repo_name: rmminusrslash/BanditsBook

max_issues_repo_name: rmminusrslash/BanditsBook

max_forks_repo_name: rmminusrslash/BanditsBook


text:
@testset "exceptions" begin

    @test_broken @interval("[+infinity]") === emptyinterval()

    @test interval(+Inf,-Inf) === emptyinterval()

    @test_broken interval(nai()) === emptyinterval()

    @test_broken @interval("[1.0000000000000001, 1.0000000000000002]") === Interval(1.0, 0x1.0000000000001p+0)

end


text:
using EffectSizes
using Test

@testset "EffectSizes.jl" begin
    include("confint.jl")
    include("effectsize.jl")
end


text:
@testset "Generated test set" begin

    dir = "integration_test_files/"

    # Iterate over a list of small problems and solutions
    for prob in eachline(dir * "list.dat")

        # Ignore blank lines
        (length(strip(prob)) == 0) && continue

        dname, sname = split(prob)

        # Data file

        open(dir * dname, "r") do fp

            global N = parse(Int, readline(fp))

            global data = readdlm(fp)[:, [1, 2]]

        end

        # Solution file
        fsol = open(dir * sname, "r")

        # Number of parameters
        n = Meta.parse(readline(fsol))

        # Solution vector
        answer = eval(Meta.parse(readline(fsol)))

        # Model function to fit data
        model = eval(Meta.parse(readline(fsol)))

        close(fsol)

        # Call raff
        rout = raff(model, data, n)

        @test rout.solution ≈ answer atol=1.0e-2

    end

end


text:
# Before running this, please make sure to activate and instantiate the environment
# corresponding to [this `Project.toml`](https://raw.githubusercontent.com/alan-turing-institute/DataScienceTutorials.jl/master/Project.toml) and [this `Manifest.toml`](https://raw.githubusercontent.com/alan-turing-institute/DataScienceTutorials.jl/master/Manifest.toml)
# so that you get an environment which matches the one used to generate the tutorials:
#
# ```julia
# cd("DataScienceTutorials") # cd to folder with the *.toml
# using Pkg; Pkg.activate("."); Pkg.instantiate()
# ```

using MLJ

@load LinearRegressor pkg=MLJLinearModels

import RDatasets: dataset
import DataFrames: describe, select, Not, rename!
boston = dataset("MASS", "Boston")
first(boston, 3)

describe(boston, :mean, :std, :eltype)

data = coerce(boston, autotype(boston, :discrete_to_continuous));

y = data.MedV
X = select(data, Not(:MedV));

mdl = LinearRegressor()

X_uni = select(X, :LStat) # only a single feature
mach_uni = machine(mdl, X_uni, y)
fit!(mach_uni)

fp = fitted_params(mach_uni)
@show fp.coefs
@show fp.intercept

using PyPlot

figure(figsize=(8,6))
plot(X.LStat, y, ls="none", marker="o")
Xnew = (LStat = collect(range(extrema(X.LStat)..., length=100)),)
plot(Xnew.LStat, predict(mach_uni, Xnew))

mach = machine(mdl, X, y)
fit!(mach)

fp = fitted_params(mach)
coefs = fp.coefs
intercept = fp.intercept
for (name, val) in coefs
    println("$(rpad(name, 8)):  $(round(val, sigdigits=3))")
end
println("Intercept: $(round(intercept, sigdigits=3))")

ŷ = predict(mach, X)
round(rms(ŷ, y), sigdigits=4)

figure(figsize=(8,6))
res = ŷ .- y
stem(res)

figure(figsize=(8,6))
hist(res, density=true)

X2 = hcat(X, X.LStat .* X.Age);

rename!(X2, :x1 => :interaction);

mach = machine(mdl, X2, y)
fit!(mach)
ŷ = predict(mach, X2)
round(rms(ŷ, y), sigdigits=4)

X3 = hcat(X.LStat, X.LStat.^2)
mach = machine(mdl, X3, y)
fit!(mach)
ŷ = predict(mach, X3)
round(rms(ŷ, y), sigdigits=4)

Xnew = (LStat = Xnew.LStat, LStat2 = Xnew.LStat.^2)

figure(figsize=(8,6))
plot(X.LStat, y, ls="none", marker="o")
plot(Xnew.LStat, predict(mach, Xnew))

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
import Base: getindex, setindex!, copy

mutable struct Board
    data::Array{Int}
    score::Int
    level::Int
    lines_to_goal::Int
    tile::Tile
    location::Vector{Int}
    orientation::Int
    nexttiles::Vector{Tile}
    holdtile::Tile
    allowhold::Bool
end
Board(i=0) = Board(ones(Int, 20, 10)*i, 0, 1+i, 5,
            rand(Tiles)(), [4,1], 0,
            [rand(Tiles)() for i in 1:3], rand(Tiles)(), true)

copy(b::Board) = Board(copy(b.data), b.score, b.level, b.lines_to_goal,
                    b.tile, b.location, b.orientation,
                    copy(b.nexttiles), b.holdtile, b.allowhold)

function getindex(b::Board)
    dy,dx = size(rotatedtile(b)) .-1
    x,y = b.location
    return b.data[y:y+dy, x:x+dx]
end

function setindex!(b::Board, s::AbstractArray)
    dy,dx = size(rotatedtile(b)) .-1
    x,y = b.location
    b.data[y:y+dy, x:x+dx] = s
end

rotatedtile(b::Board) = rotr90(data(b.tile), b.orientation)

function nexttile!(b::Board)
    push!(b.nexttiles, rand(Tiles)())
    b.tile = popfirst!(b.nexttiles)
end

function add_tile!(b::Board)
    nexttile!(b)
    b.location = start_location(b.tile)
    b.orientation = 0
    if all(b[] .== 0)
        oldboard = copy(b)
        b[] += data(b.tile)
        update_board!(oldboard, b)
        return true
    end
    false
end

function delete_lines!(b::Board)
    oldboard = copy(b)
    nr_lines = 0
    for i in 1:20
        if all(b.data[i, :] .!= 0)
            b.data[2:i, :] = b.data[1:i-1, :]
            b.data[1,:] .= 0
            nr_lines += 1
        end
    end
    b.lines_to_goal -= nr_lines
    b.score += [0 1 3 5 8][nr_lines+1] * b.level * 100
    if b.lines_to_goal ≤ 0
        b.level += 1
        b.lines_to_goal += b.level*5
    end
    update_board!(oldboard, b)
end

max_stars_repo_name: ChristianKurz/REPLTetris

max_issues_repo_name: ChristianKurz/REPLTetris

max_forks_repo_name: ChristianKurz/REPLTetris


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using AbaqusReader

"""
    abaqus_read_mesh(fn::String)

Read and parse ABAQUS `.inp` file.

`fn` (filename) is the name of the file to parse.
"""
function abaqus_read_mesh(fn::String)
    m = AbaqusReader.abaqus_read_mesh(fn)
    return Mesh(m)
end

"""
    create_surface_elements(mesh::Mesh, surface_name::Symbol)

Create a set of surface elements from solid elements.

Notation follow what is defined in ABAQUS. For example, if solid elements
are Tet10, surface elements will be Tri6 and they can be used to define
boundary conditions.
"""
function create_surface_elements(mesh::Mesh, surface_name::Symbol)
    elements = Element[]
    for (elid, elsi) in mesh.surface_sets[surface_name]
        elty = mesh.element_types[elid]
        elco = mesh.elements[elid]
        chel, chcon = AbaqusReader.create_surface_element(elty, elsi, elco)
        ch = Element(getfield(JuliaFEM, chel), chcon)
        push!(elements, ch)
    end
    update!(elements, "geometry", mesh.nodes)
    return elements
end

"""
    create_surface_elements(mesh::Mesh, surface_name::String)

Create a set of surface elements from solid elements.

Notation follow what is defined in ABAQUS. For example, if solid elements
are Tet10, surface elements will be Tri6 and they can be used to define
boundary conditions.
"""
function create_surface_elements(mesh::Mesh, surface_name::String)
    return create_surface_elements(mesh, Symbol(surface_name))
end

"""
    create_nodal_elements(mesh::Mesh, node_set_name::String)

Create a set of "nodal" elements from node set.

Nodal elements are of type Poi1. In principle they don't have volume and it's
not possible to integrate over them. It is however possible to add boundary
conditions to them.
"""
function create_nodal_elements(mesh::Mesh, node_set_name::String)
    node_ids = mesh.node_sets[Symbol(node_set_name)]
    elements = [Element(Poi1, [j]) for j in node_ids]
    update!(elements, "geometry", mesh.nodes)
    return elements
end


text:
__precompile__()

module SparseRelaxation

export gauss_seidel, gauss_seidel!, weighted_jacobi!, weighted_jacobi

# TODO: how to improve locality
# order is different than it should be?
"""
    gauss_seidel(A::SparseMatrixCSC, x::Vector, b::Vector; backwards::Bool = false, iterations = 1)

Return a vector that is the result of applying one iteration of Gauss Seidel to
`Ax=b`. This iteration is equivalent to `L(b - Ux)` where `L` is the lower
triangular part of `A` and `U` is the strictly upper triangular part of `A`. If
`backwards` is true, then a backwards sweep is performed. `A` is assumed to be symmetric. If not, then `gauss_seidel(A)` is really Gauss-Seidel applied to `A^T`.
"""
function gauss_seidel{T}( A :: SparseMatrixCSC{T}
                        , x :: Vector{T}
                        , b :: Vector{T}
                        ; backwards :: Bool = false
                        , iterations :: Int = 1
                        ) :: Vector{T}
    x_new = copy(x)
    gauss_seidel!(A, x_new, b, backwards=backwards, iterations=iterations)
    x_new
end

"""
    gauss_seidel!(A::SparseMatrixCSC, x::Vector, b::Vector, backwards::Bool = false)

In place, more efficient version of `gauss_seidel(A, x, b)`.
"""
function gauss_seidel!{T}( A :: SparseMatrixCSC{T}
                         , x :: Vector{T}
                         , b :: Vector{T}
                         ; backwards :: Bool = false
                         , iterations :: Int = 1
                         ) :: Vector{T}
    rows = rowvals(A)
    vals = nonzeros(A)
    m, n = size(A)
    for col in 1:n
        colsum :: T = 0
        diag :: T = 0

        js = nzrange(A, col)
        # hack to avoid allocating
        start = backwards ? js.stop : js.start
        stop = backwards ? js.start : js.stop
        step = backwards ? -1 : 1
        @inbounds for j in start:step:stop
            row = rows[j]
            val = vals[j]
            if row != col
                colsum += val * x[row]
            else
                diag = val
            end
        end

        @inbounds x[col] = (b[col] - colsum) / diag
    end

    x
end

"""
In place version of `weighted_jacobi`.
"""
function weighted_jacobi!{T}( A :: SparseMatrixCSC{T}
                            , x :: Vector{T}
                            , b :: Vector{T}
                            ; weight :: T = 2/3
                            , iterations :: Int = 1
                            ) :: Vector{T}
    Di = spdiagm(map(x -> 1/x, diag(A)))
    R = A - spdiagm(diag(A))
    for i in 1:iterations
        x[:] = weight * Di * (b - R*x) + (1 - weight) * x
    end

    x
end

"""
    weighted_jacobi(A::SparseMatrixCSC, x::Vector, b::Vector; weight = 2/3, iterations = 1)

Return a vector that is the result of applying one iteration of Jacobi smoothing to
`Ax=b`. This iteration is equivalent to `wD(b - R*x) + (1-w)x` where `D` is the inverse diagonal and `R` is `A` without its diagonal.
"""
function weighted_jacobi{T}( A :: SparseMatrixCSC{T}
                           , x :: Vector{T}
                           , b :: Vector{T}
                           ; weight :: T = 2/3
                           , iterations :: Int = 1
                           ) :: Vector{T}
    x_new = copy(x)
    weighted_jacobi!(A, x_new, b, weight=weight, iterations=iterations)
    x_new
end

end


text:
using DrWatson
@quickactivate "FractalDimension"

# Import all needed functions for benchmarks from packages and by direct definition.
using BenchmarkTools, DynamicalSystems, Neighborhood, Entropies, DataFrames, CSV

function correlationsum_new(X, ε; q = 2, metric = Euclidean(), w = 0)
    ntype = WithinRange(maximum(ε))
    tree = KDTree(X)
    qf = q == 2 ? identity : x -> float(x)^(q-1)
    _correlationsum_new(tree, X, ε, w, qf)
end

function _correlationsum_new(tree, X, ε, w, qf)
    theiler = Neighborhood.Theiler(w)
    searchtype = Neighborhood.WithinRange(maximum(ε))
    N = length(X)
    Cq = similar(ε, Float64)
    @inbounds for i in 1:N
        idxs = Neighborhood.isearch(tree, X[i], searchtype, theiler(i))
        isempty(idxs) && continue
        dist = [Neighborhood.evaluate(tree.metric, X[i], X[j]) for j in idxs]
        sort!(dist)
        for k in 1:length(ε)-1
            j = searchsortedfirst(dist, ε[k])
            Cq[k] += qf(j-1)
        end
        Cq[end] += qf(length(dist))
    end
    return Cq ./ (N*qf(N-w-1))
end

# Define a function to access data sets.
"""
	sample_data(data, N_sample)
Sample `N_sample` unique random points out of `data`.
"""
function sample_data(data, N_sample)
    N_data = length(data)
    possible_indices = 1:N_data
    indices = unique(rand(possible_indices, 2N_sample))[1:N_sample]
    data[indices] |> Dataset
end

# Define Benchmark functions.
# All functions can output probabilities so the time is benchmarked for probabilities.
function benchmark_correlation(data)
    bm = @benchmark begin
        ε0 = min_pairwise_distance($data)[2]
        ε1 = minimum(maxima($data) - minima($data))
        εs = 10 .^ range(log10(ε0), log10(ε1), length = 12)
        [correlationsum($data, ε) for ε in εs]
    end
    minimum(bm.times)
end

function benchmark_new_correlation(data)
    bm = @benchmark begin
        ε0 = min_pairwise_distance($data)[2]
        ε1 = minimum(maxima($data) - minima($data))
        εs = 10 .^ range(log10(ε0), log10(ε1), length = 12)
        correlationsum_new($data, εs)
    end
    minimum(bm.times)
end

function benchmark_boxed_correlation(data::Dataset{D,T}, prismdim = D) where {D,T}
    bm = @benchmark begin
        ε0 = min_pairwise_distance($data)[2]
        ε1 = estimate_r0_buenoorovio($data, $prismdim)
        εs = 10 .^ range(log10(ε0), log10(ε1), length = 12)
        boxed_correlationsum($data, εs, ε1, M = $prismdim)
    end
    minimum(bm.times)
end

function benchmark_molteno(data)
    bm = @benchmark molteno_boxing($data)
    minimum(bm.times)
end

function benchmark_genentropy(data)
    bm = @benchmark begin
        ε0 = min_pairwise_distance($data)[2]
        ε1 = minimum(maxima($data) - minima($data))
        εs = 10 .^ range(log10(ε0), log10(ε1), length = 12)
		[Entropies.genentropy($data, ε) for ε in εs]
	end
	minimum(bm.times)
end

# Define size benchmark with logarithmically spread values.
# `base = 1` allows non-logarithmic values
function benchmark_N(
	lower_number,
	upper_number,
	number_of_numbers;
	base = 10,
	feedback = false
)
	result = DataFrame(N = Int[], Model = [], Method = [], Time = [])
	numbers = floor.(Int, base .^ range(lower_number, upper_number, length = number_of_numbers))

	# Define Systems.
	henon_sys = Systems.henon(rand(2))
	lorenz_sys = Systems.lorenz(randn(3))
	N_traj = 5 * floor(base ^ upper_number)

	# Integrate trajectories.
	henon = trajectory(henon_sys, N_traj, Ttr = 100)
	lorenz = trajectory(lorenz_sys, N_traj * 0.1, dt = 0.1, Ttr = 10)

	trajectories = [(:henon, henon), (:lorenz, lorenz)]

	for (name, traj) in trajectories
		feedback && println("Now using the $name system.")
		for N in numbers
			feedback && println("Currently calculating with $N random data points.")
			data = sample_data(traj, N)
			feedback && println("Benchmarking the classic correlation…")
			push!(
				result,
				(N, name, "Correlation",benchmark_correlation(data))
			)
			feedback && println("Benchmarking the new correlation…")
			push!(
				result,
				(N, name, "New Correlation", benchmark_new_correlation(data))
			)
			feedback && println("Benchmarking the boxed correlation…")
			push!(
				result,
				(N, name, "Boxed Correlation", benchmark_boxed_correlation(data))
			)
			feedback && println("Benchmarking the generalized entropy with factor 2 boxing…")
			push!(
				result,
				(N, name, "Molteno", benchmark_molteno(data))
			)
			feedback && println("Benchmarking the generalized entropy with standard boxing…")
			push!(
				result,
				(N, name, "Generalized Dimension", benchmark_genentropy(data))
			)
		end
	end
	wsave(
		datadir(
			"benchmark",
			savename(
				"benchmark_number",
				Dict(
					:lowerN => lower_number,
					:upperN => upper_number,
					:NofN => number_of_numbers,
					:base => base,
				),
				"csv",
			),
		),
		result,
	)
	result
end

# Define Dimension Benchmark.
function benchmark_dimension(
	lower_dimension::Int = 4,
	upper_dimension::Int = 16,
	distance_dimension::Int = 1;
	N_sample = 10000,
	feedback = false,
)
	result = DataFrame(Dimension = Int[], Method = [], Time = [])
	dimensions = lower_dimension:distance_dimension:upper_dimension
	for dimension in dimensions
		feedback && println("The current dimension is $dimension.")
		# The current trajectroy is longer than the data sample taken from the trajectory.
		traj = trajectory(
			Systems.lorenz96(dimension, rand(dimension), F = 8.0),
			N_sample,
			Ttr = 100,
			dt = 0.1,
		)
		# Only 10000 unique data points are sampled from the 100000 trajectory points
		data = sample_data(traj, N_sample)
		# Benchmark the classic correlation.
		feedback && println("Now benchmarking the classic correlation…")
		push!(
			result,
			(dimension, "Classic Correlation", benchmark_correlation(data)),
		)
		# Benchmark the new correlation that uses KDTrees.
		feedback && println("Now benchmarking the new correlation…")
		push!(
			result,
			(dimension, "New Correlation", benchmark_new_correlation(data)),
		)
		# Benchmark the Boxed Correlation without accessing the prism property.
		feedback && println("Now benchmarking the boxed correlation…")
		push!(
			result,
			(dimension, "Boxed Correlation", benchmark_boxed_correlation(data)),
		)
		feedback && println("Now benchmarking the prism correlation…")
		# Benchmark the Prism Correlation with Prism Dimension by Theiler.
		embedding_dim =
			0.75log2(N_sample) < dimension ?
			(0.5log2(N_sample) |> ceil |> Int) :
			dimension
		push!(
			result,
			(dimension, "Prism Correlation", benchmark_boxed_correlation(data, embedding_dim)),
		)
		feedback && println("Now benchmarking the generalized entropy by Molteno…")
		# Benchmark the generalized entropy method by Molteno.
		push!(
			result,
			(dimension, "Molteno", benchmark_molteno(data))
		)
		feedback && println("Now benchmarking the classic generalized entropy…")
		# Benchmark the classic generalized entropy.
		push!(
			result,
			(dimension, "Generalized Entropy", benchmark_genentropy(data))
		)
	end
	wsave(
		datadir(
			"benchmark",
			savename(
				"benchmark_dim",
				Dict(
					:N => N_sample,
					:lowerdim => lower_dimension,
					:upperdim => upper_dimension,
					:ddim => distance_dimension,
				),
				"csv",
			),
		),
		result,
	)
	result
end

# df_dim = benchmark_dimension()
#

max_stars_repo_name: Datseris/FractalDimension

max_issues_repo_name: Datseris/FractalDimension

max_forks_repo_name: Datseris/FractalDimension


text:
# NOTE: This file contains machine-generated code by translate_exported_functions.jl.

# sam_hdr_t *sam_hdr_init(void);
@defun sam_hdr_init()::Ptr{sam_hdr_t}

# sam_hdr_t *bam_hdr_read(BGZF *fp);
@defun bam_hdr_read(fp::Ptr{BGZF})::Ptr{sam_hdr_t}

# int bam_hdr_write(BGZF *fp, const sam_hdr_t *h) HTS_RESULT_USED;
@defun bam_hdr_write(fp::Ptr{BGZF}, h::Ptr{sam_hdr_t})::Cint

# void sam_hdr_destroy(sam_hdr_t *h);
@defun sam_hdr_destroy(h::Ptr{sam_hdr_t})::Cvoid

# sam_hdr_t *sam_hdr_dup(const sam_hdr_t *h0);
@defun sam_hdr_dup(h0::Ptr{sam_hdr_t})::Ptr{sam_hdr_t}

# sam_hdr_t *sam_hdr_parse(size_t l_text, const char *text);
@defun sam_hdr_parse(l_text::Csize_t, text::Cstring)::Ptr{sam_hdr_t}

# sam_hdr_t *sam_hdr_read(samFile *fp);
@defun sam_hdr_read(fp::Ptr{samFile})::Ptr{sam_hdr_t}

# int sam_hdr_write(samFile *fp, const sam_hdr_t *h) HTS_RESULT_USED;
@defun sam_hdr_write(fp::Ptr{samFile}, h::Ptr{sam_hdr_t})::Cint

# size_t sam_hdr_length(sam_hdr_t *h);
@defun sam_hdr_length(h::Ptr{sam_hdr_t})::Csize_t

# const char *sam_hdr_str(sam_hdr_t *h);
@defun sam_hdr_str(h::Ptr{sam_hdr_t})::Cstring

# int sam_hdr_nref(const sam_hdr_t *h);
@defun sam_hdr_nref(h::Ptr{sam_hdr_t})::Cint

# int sam_hdr_add_lines(sam_hdr_t *h, const char *lines, size_t len);
@defun sam_hdr_add_lines(h::Ptr{sam_hdr_t}, lines::Cstring, len::Csize_t)::Cint

# int sam_hdr_add_line(sam_hdr_t *h, const char *type, ...);
# ErrorException("unknown token: '.'")

# int sam_hdr_find_line_id(sam_hdr_t *h, const char *type, const char *ID_key, const char *ID_val, kstring_t *ks);
@defun sam_hdr_find_line_id(h::Ptr{sam_hdr_t}, type::Cstring, ID_key::Cstring, ID_val::Cstring, ks::Ptr{kstring_t})::Cint

# int sam_hdr_find_line_pos(sam_hdr_t *h, const char *type, int pos, kstring_t *ks);
@defun sam_hdr_find_line_pos(h::Ptr{sam_hdr_t}, type::Cstring, pos::Cint, ks::Ptr{kstring_t})::Cint

# int sam_hdr_remove_line_id(sam_hdr_t *h, const char *type, const char *ID_key, const char *ID_value);
@defun sam_hdr_remove_line_id(h::Ptr{sam_hdr_t}, type::Cstring, ID_key::Cstring, ID_value::Cstring)::Cint

# int sam_hdr_remove_line_pos(sam_hdr_t *h, const char *type, int position);
@defun sam_hdr_remove_line_pos(h::Ptr{sam_hdr_t}, type::Cstring, position::Cint)::Cint

# int sam_hdr_update_line(sam_hdr_t *h, const char *type, const char *ID_key, const char *ID_value, ...);
# ErrorException("unknown token: '.'")

# int sam_hdr_remove_except(sam_hdr_t *h, const char *type, const char *ID_key, const char *ID_value);
@defun sam_hdr_remove_except(h::Ptr{sam_hdr_t}, type::Cstring, ID_key::Cstring, ID_value::Cstring)::Cint

# int sam_hdr_remove_lines(sam_hdr_t *h, const char *type, const char *id, void *rh);
@defun sam_hdr_remove_lines(h::Ptr{sam_hdr_t}, type::Cstring, id::Cstring, rh::Ptr{Cvoid})::Cint

# int sam_hdr_count_lines(sam_hdr_t *h, const char *type);
@defun sam_hdr_count_lines(h::Ptr{sam_hdr_t}, type::Cstring)::Cint

# int sam_hdr_line_index(sam_hdr_t *bh, const char *type, const char *key);
@defun sam_hdr_line_index(bh::Ptr{sam_hdr_t}, type::Cstring, key::Cstring)::Cint

# const char *sam_hdr_line_name(sam_hdr_t *bh, const char *type, int pos);
@defun sam_hdr_line_name(bh::Ptr{sam_hdr_t}, type::Cstring, pos::Cint)::Cstring

# int sam_hdr_find_tag_id(sam_hdr_t *h, const char *type, const char *ID_key, const char *ID_value, const char *key, kstring_t *ks);
@defun sam_hdr_find_tag_id(h::Ptr{sam_hdr_t}, type::Cstring, ID_key::Cstring, ID_value::Cstring, key::Cstring, ks::Ptr{kstring_t})::Cint

# int sam_hdr_find_tag_pos(sam_hdr_t *h, const char *type, int pos, const char *key, kstring_t *ks);
@defun sam_hdr_find_tag_pos(h::Ptr{sam_hdr_t}, type::Cstring, pos::Cint, key::Cstring, ks::Ptr{kstring_t})::Cint

# int sam_hdr_remove_tag_id(sam_hdr_t *h, const char *type, const char *ID_key, const char *ID_value, const char *key);
@defun sam_hdr_remove_tag_id(h::Ptr{sam_hdr_t}, type::Cstring, ID_key::Cstring, ID_value::Cstring, key::Cstring)::Cint

# int sam_hdr_name2tid(sam_hdr_t *h, const char *ref);
@defun sam_hdr_name2tid(h::Ptr{sam_hdr_t}, ref::Cstring)::Cint

# const char *sam_hdr_tid2name(const sam_hdr_t *h, int tid);
@defun sam_hdr_tid2name(h::Ptr{sam_hdr_t}, tid::Cint)::Cstring

# hts_pos_t sam_hdr_tid2len(const sam_hdr_t *h, int tid);
@defun sam_hdr_tid2len(h::Ptr{sam_hdr_t}, tid::Cint)::hts_pos_t

# const char *sam_hdr_pg_id(sam_hdr_t *h, const char *name);
@defun sam_hdr_pg_id(h::Ptr{sam_hdr_t}, name::Cstring)::Cstring

# int sam_hdr_add_pg(sam_hdr_t *h, const char *name, ...);
# ErrorException("unknown token: '.'")

# char *stringify_argv(int argc, char *argv[]);
@defun stringify_argv(argc::Cint, argv::Ptr{Ptr{Cchar}})::Ptr{Cchar}

# void sam_hdr_incr_ref(sam_hdr_t *h);
@defun sam_hdr_incr_ref(h::Ptr{sam_hdr_t})::Cvoid

# bam1_t *bam_init1(void);
@defun bam_init1()::Ptr{bam1_t}

# void bam_destroy1(bam1_t *b);
@defun bam_destroy1(b::Ptr{bam1_t})::Cvoid

# int bam_read1(BGZF *fp, bam1_t *b) HTS_RESULT_USED;
@defun bam_read1(fp::Ptr{BGZF}, b::Ptr{bam1_t})::Cint

# int bam_write1(BGZF *fp, const bam1_t *b) HTS_RESULT_USED;
@defun bam_write1(fp::Ptr{BGZF}, b::Ptr{bam1_t})::Cint

# bam1_t *bam_copy1(bam1_t *bdst, const bam1_t *bsrc) HTS_RESULT_USED;
@defun bam_copy1(bdst::Ptr{bam1_t}, bsrc::Ptr{bam1_t})::Ptr{bam1_t}

# bam1_t *bam_dup1(const bam1_t *bsrc);
@defun bam_dup1(bsrc::Ptr{bam1_t})::Ptr{bam1_t}

# hts_pos_t bam_cigar2qlen(int n_cigar, const uint32_t *cigar);
@defun bam_cigar2qlen(n_cigar::Cint, cigar::Ptr{UInt32})::hts_pos_t

# hts_pos_t bam_cigar2rlen(int n_cigar, const uint32_t *cigar);
@defun bam_cigar2rlen(n_cigar::Cint, cigar::Ptr{UInt32})::hts_pos_t

# hts_pos_t bam_endpos(const bam1_t *b);
@defun bam_endpos(b::Ptr{bam1_t})::hts_pos_t

# int bam_str2flag(const char *str);
@defun bam_str2flag(str::Cstring)::Cint

# char *bam_flag2str(int flag);
@defun bam_flag2str(flag::Cint)::Ptr{Cchar}

# int bam_set_qname(bam1_t *b, const char *qname);
@defun bam_set_qname(b::Ptr{bam1_t}, qname::Cstring)::Cint

# int sam_idx_init(htsFile *fp, sam_hdr_t *h, int min_shift, const char *fnidx);
@defun sam_idx_init(fp::Ptr{htsFile}, h::Ptr{sam_hdr_t}, min_shift::Cint, fnidx::Cstring)::Cint

# int sam_idx_save(htsFile *fp) HTS_RESULT_USED;
@defun sam_idx_save(fp::Ptr{htsFile})::Cint

# hts_idx_t *sam_index_load(htsFile *fp, const char *fn);
@defun sam_index_load(fp::Ptr{htsFile}, fn::Cstring)::Ptr{hts_idx_t}

# hts_idx_t *sam_index_load2(htsFile *fp, const char *fn, const char *fnidx);
@defun sam_index_load2(fp::Ptr{htsFile}, fn::Cstring, fnidx::Cstring)::Ptr{hts_idx_t}

# hts_idx_t *sam_index_load3(htsFile *fp, const char *fn, const char *fnidx, int flags);
@defun sam_index_load3(fp::Ptr{htsFile}, fn::Cstring, fnidx::Cstring, flags::Cint)::Ptr{hts_idx_t}

# int sam_index_build(const char *fn, int min_shift) HTS_RESULT_USED;
@defun sam_index_build(fn::Cstring, min_shift::Cint)::Cint

# int sam_index_build2(const char *fn, const char *fnidx, int min_shift) HTS_RESULT_USED;
@defun sam_index_build2(fn::Cstring, fnidx::Cstring, min_shift::Cint)::Cint

# int sam_index_build3(const char *fn, const char *fnidx, int min_shift, int nthreads) HTS_RESULT_USED;
@defun sam_index_build3(fn::Cstring, fnidx::Cstring, min_shift::Cint, nthreads::Cint)::Cint

# hts_itr_t *sam_itr_queryi(const hts_idx_t *idx, int tid, hts_pos_t beg, hts_pos_t end);
@defun sam_itr_queryi(idx::Ptr{hts_idx_t}, tid::Cint, beg::hts_pos_t, end_::hts_pos_t)::Ptr{hts_itr_t}

# hts_itr_t *sam_itr_querys(const hts_idx_t *idx, sam_hdr_t *hdr, const char *region);
@defun sam_itr_querys(idx::Ptr{hts_idx_t}, hdr::Ptr{sam_hdr_t}, region::Cstring)::Ptr{hts_itr_t}

# hts_itr_t *sam_itr_regions(const hts_idx_t *idx, sam_hdr_t *hdr, hts_reglist_t *reglist, unsigned int regcount);
@defun sam_itr_regions(idx::Ptr{hts_idx_t}, hdr::Ptr{sam_hdr_t}, reglist::Ptr{hts_reglist_t}, regcount::Cuint)::Ptr{hts_itr_t}

# hts_itr_t *sam_itr_regarray(const hts_idx_t *idx, sam_hdr_t *hdr, char **regarray, unsigned int regcount);
@defun sam_itr_regarray(idx::Ptr{hts_idx_t}, hdr::Ptr{sam_hdr_t}, regarray::Ptr{Ptr{Cchar}}, regcount::Cuint)::Ptr{hts_itr_t}

# const char *sam_parse_region(sam_hdr_t *h, const char *s, int *tid, hts_pos_t *beg, hts_pos_t *end, int flags);
@defun sam_parse_region(h::Ptr{sam_hdr_t}, s::Cstring, tid::Ptr{Cint}, beg::Ptr{hts_pos_t}, end_::Ptr{hts_pos_t}, flags::Cint)::Cstring

# int sam_open_mode(char *mode, const char *fn, const char *format);
@defun sam_open_mode(mode::Ptr{Cchar}, fn::Cstring, format::Cstring)::Cint

# char *sam_open_mode_opts(const char *fn, const char *mode, const char *format);
@defun sam_open_mode_opts(fn::Cstring, mode::Cstring, format::Cstring)::Ptr{Cchar}

# int sam_hdr_change_HD(sam_hdr_t *h, const char *key, const char *val);
@defun sam_hdr_change_HD(h::Ptr{sam_hdr_t}, key::Cstring, val::Cstring)::Cint

# int sam_parse1(kstring_t *s, sam_hdr_t *h, bam1_t *b) HTS_RESULT_USED;
@defun sam_parse1(s::Ptr{kstring_t}, h::Ptr{sam_hdr_t}, b::Ptr{bam1_t})::Cint

# int sam_format1(const sam_hdr_t *h, const bam1_t *b, kstring_t *str) HTS_RESULT_USED;
@defun sam_format1(h::Ptr{sam_hdr_t}, b::Ptr{bam1_t}, str::Ptr{kstring_t})::Cint

# int sam_read1(samFile *fp, sam_hdr_t *h, bam1_t *b) HTS_RESULT_USED;
@defun sam_read1(fp::Ptr{samFile}, h::Ptr{sam_hdr_t}, b::Ptr{bam1_t})::Cint

# int sam_write1(samFile *fp, const sam_hdr_t *h, const bam1_t *b) HTS_RESULT_USED;
@defun sam_write1(fp::Ptr{samFile}, h::Ptr{sam_hdr_t}, b::Ptr{bam1_t})::Cint

# uint8_t *bam_aux_get(const bam1_t *b, const char tag[2]);
@defun bam_aux_get(b::Ptr{bam1_t}, tag::Ptr{Cchar})::Ptr{UInt8}

# int64_t bam_aux2i(const uint8_t *s);
@defun bam_aux2i(s::Ptr{UInt8})::Int64

# double bam_aux2f(const uint8_t *s);
@defun bam_aux2f(s::Ptr{UInt8})::Cdouble

# char bam_aux2A(const uint8_t *s);
@defun bam_aux2A(s::Ptr{UInt8})::Cchar

# char *bam_aux2Z(const uint8_t *s);
@defun bam_aux2Z(s::Ptr{UInt8})::Ptr{Cchar}

# uint32_t bam_auxB_len(const uint8_t *s);
@defun bam_auxB_len(s::Ptr{UInt8})::UInt32

# int64_t bam_auxB2i(const uint8_t *s, uint32_t idx);
@defun bam_auxB2i(s::Ptr{UInt8}, idx::UInt32)::Int64

# double bam_auxB2f(const uint8_t *s, uint32_t idx);
@defun bam_auxB2f(s::Ptr{UInt8}, idx::UInt32)::Cdouble

# int bam_aux_append(bam1_t *b, const char tag[2], char type, int len, const uint8_t *data);
@defun bam_aux_append(b::Ptr{bam1_t}, tag::Ptr{Cchar}, type::Cchar, len::Cint, data::Ptr{UInt8})::Cint

# int bam_aux_del(bam1_t *b, uint8_t *s);
@defun bam_aux_del(b::Ptr{bam1_t}, s::Ptr{UInt8})::Cint

# int bam_aux_update_str(bam1_t *b, const char tag[2], int len, const char *data);
@defun bam_aux_update_str(b::Ptr{bam1_t}, tag::Ptr{Cchar}, len::Cint, data::Cstring)::Cint

# int bam_aux_update_int(bam1_t *b, const char tag[2], int64_t val);
@defun bam_aux_update_int(b::Ptr{bam1_t}, tag::Ptr{Cchar}, val::Int64)::Cint

# int bam_aux_update_float(bam1_t *b, const char tag[2], float val);
@defun bam_aux_update_float(b::Ptr{bam1_t}, tag::Ptr{Cchar}, val::Cfloat)::Cint

# int bam_aux_update_array(bam1_t *b, const char tag[2], uint8_t type, uint32_t items, void *data);
@defun bam_aux_update_array(b::Ptr{bam1_t}, tag::Ptr{Cchar}, type::UInt8, items::UInt32, data::Ptr{Cvoid})::Cint

# bam_plp_t bam_plp_init(bam_plp_auto_f func, void *data);
@defun bam_plp_init(func::bam_plp_auto_f, data::Ptr{Cvoid})::bam_plp_t

# void bam_plp_destroy(bam_plp_t iter);
@defun bam_plp_destroy(iter::bam_plp_t)::Cvoid

# int bam_plp_push(bam_plp_t iter, const bam1_t *b);
@defun bam_plp_push(iter::bam_plp_t, b::Ptr{bam1_t})::Cint

# const bam_pileup1_t *bam_plp_next(bam_plp_t iter, int *_tid, int *_pos, int *_n_plp);
@defun bam_plp_next(iter::bam_plp_t, _tid::Ptr{Cint}, _pos::Ptr{Cint}, _n_plp::Ptr{Cint})::Ptr{bam_pileup1_t}

# const bam_pileup1_t *bam_plp_auto(bam_plp_t iter, int *_tid, int *_pos, int *_n_plp);
@defun bam_plp_auto(iter::bam_plp_t, _tid::Ptr{Cint}, _pos::Ptr{Cint}, _n_plp::Ptr{Cint})::Ptr{bam_pileup1_t}

# const bam_pileup1_t *bam_plp64_next(bam_plp_t iter, int *_tid, hts_pos_t *_pos, int *_n_plp);
@defun bam_plp64_next(iter::bam_plp_t, _tid::Ptr{Cint}, _pos::Ptr{hts_pos_t}, _n_plp::Ptr{Cint})::Ptr{bam_pileup1_t}

# const bam_pileup1_t *bam_plp64_auto(bam_plp_t iter, int *_tid, hts_pos_t *_pos, int *_n_plp);
@defun bam_plp64_auto(iter::bam_plp_t, _tid::Ptr{Cint}, _pos::Ptr{hts_pos_t}, _n_plp::Ptr{Cint})::Ptr{bam_pileup1_t}

# void bam_plp_set_maxcnt(bam_plp_t iter, int maxcnt);
@defun bam_plp_set_maxcnt(iter::bam_plp_t, maxcnt::Cint)::Cvoid

# void bam_plp_reset(bam_plp_t iter);
@defun bam_plp_reset(iter::bam_plp_t)::Cvoid

# void bam_plp_constructor(bam_plp_t plp, int (*func)(void *data, const bam1_t *b, bam_pileup_cd *cd));
# ErrorException("function pointer?")

# void bam_plp_destructor(bam_plp_t plp, int (*func)(void *data, const bam1_t *b, bam_pileup_cd *cd));
# ErrorException("function pointer?")

# int bam_plp_insertion(const bam_pileup1_t *p, kstring_t *ins, int *del_len) HTS_RESULT_USED;
@defun bam_plp_insertion(p::Ptr{bam_pileup1_t}, ins::Ptr{kstring_t}, del_len::Ptr{Cint})::Cint

# bam_mplp_t bam_mplp_init(int n, bam_plp_auto_f func, void **data);
@defun bam_mplp_init(n::Cint, func::bam_plp_auto_f, data::Ptr{Ptr{Cvoid}})::bam_mplp_t

# int bam_mplp_init_overlaps(bam_mplp_t iter);
@defun bam_mplp_init_overlaps(iter::bam_mplp_t)::Cint

# void bam_mplp_destroy(bam_mplp_t iter);
@defun bam_mplp_destroy(iter::bam_mplp_t)::Cvoid

# void bam_mplp_set_maxcnt(bam_mplp_t iter, int maxcnt);
@defun bam_mplp_set_maxcnt(iter::bam_mplp_t, maxcnt::Cint)::Cvoid

# int bam_mplp_auto(bam_mplp_t iter, int *_tid, int *_pos, int *n_plp, const bam_pileup1_t **plp);
@defun bam_mplp_auto(iter::bam_mplp_t, _tid::Ptr{Cint}, _pos::Ptr{Cint}, n_plp::Ptr{Cint}, plp::Ptr{Ptr{bam_pileup1_t}})::Cint

# int bam_mplp64_auto(bam_mplp_t iter, int *_tid, hts_pos_t *_pos, int *n_plp, const bam_pileup1_t **plp);
@defun bam_mplp64_auto(iter::bam_mplp_t, _tid::Ptr{Cint}, _pos::Ptr{hts_pos_t}, n_plp::Ptr{Cint}, plp::Ptr{Ptr{bam_pileup1_t}})::Cint

# void bam_mplp_reset(bam_mplp_t iter);
@defun bam_mplp_reset(iter::bam_mplp_t)::Cvoid

# void bam_mplp_constructor(bam_mplp_t iter, int (*func)(void *data, const bam1_t *b, bam_pileup_cd *cd));
# ErrorException("function pointer?")

# void bam_mplp_destructor(bam_mplp_t iter, int (*func)(void *data, const bam1_t *b, bam_pileup_cd *cd));
# ErrorException("function pointer?")

# int sam_cap_mapq(bam1_t *b, const char *ref, hts_pos_t ref_len, int thres);
@defun sam_cap_mapq(b::Ptr{bam1_t}, ref::Cstring, ref_len::hts_pos_t, thres::Cint)::Cint

# int sam_prob_realn(bam1_t *b, const char *ref, hts_pos_t ref_len, int flag);
@defun sam_prob_realn(b::Ptr{bam1_t}, ref::Cstring, ref_len::hts_pos_t, flag::Cint)::Cint


text:
export Lowrank

struct Lowrank{T <: Real} <: ConvexSet
    k::T
end
```
This function projects the matrix y onto the closest matrix of rank k or less.
```
function project!(s::Lowrank{<:Real}, v::Matrix{T}, y::Matrix{T}) where {T <: Real}
    U, S, V = LinearAlgebra.svd(y)
    Z = zeros(size(y))
    for j = 1:s.k
      Z = Z .+ S[j]*U[:,j]*V[:,j]'
    end
    v = Z
end

function project(s::Lowrank{<:Real}, y::Matrix{<:Real})
    project!(s, similar(y), y)
end


text:
@testset "208.implement-trie-prefix-tree.jl" begin
    @testset "208 case 1" begin
        node = PrefixNode()
        insert_node!(node, "apple")
        @test search_word(node, "apple")
        @test !search_word(node, "app")
        @test starts_with(node, "app")
        insert_node!(node, "app")
        @test search_word(node, "app")
    end
end


text:
try
    verstring = read(`genhtml -v`, String)
    m = match(r"^.*(\d.\d+)$", verstring)
    if m ≡ nothing
        @warn "Could not parse `genhtml -v` version information, please open an issue."
    else
        ver = VersionNumber(m[1])
        OLDVERWARN = """
        The installed version of `genhtml` is $ver, HTML generation may not work.
        See https://github.com/tpapp/LocalCoverage.jl/issues/1 for workarounds.
        """
        ver ≥ v"1.12" || warn(OLDVERWARN)
    end
catch
    @warn "Could not find `genhtml`, you need to install it to generate HTML reports."
end


text:
using Test: @testset, @test

using Agents: allagents, run!, dummystep
using DataFrames: DataFrame
using Statistics: mean
using LightGraphs: SimpleGraph, barabasi_albert, add_edge!
using GraphPlot: gplot

include("../src/agent.jl")
include("../src/model.jl")
include("../src/exec.jl")

println("Julia $VERSION")
println("Thread count: $(Threads.nthreads())")

G = barabasi_albert(10^3, 2)
println(G)
println("2 * ne(G) / nv(G) = $(2 * ne(G) / nv(G))")

@testset "agent.Player" begin
    @test Agent.Player(Agent.C) isa Agent.Player
    @test Agent.Player(1, Agent.D) isa Agent.Player
    @test Agent.Player(1, 0.5) isa Agent.Player
end

@testset "model.build_model" begin
    model = Model.build_model(G = barabasi_albert(10^3, 4), cost_model = :fixed_cost, r = 2.0, C_rate = 0.7)
    @test length(allagents(model)) == 10^3
    @test 0.65 < mean([a.strategy for a in allagents(model)]) < 0.75
end

@testset "model.model_step!" begin
    @testset "executable" begin
        G = barabasi_albert(10^3, 4)
        model = Model.build_model(G = G, r = calc_r(G, 0.8), cost_model = :fixed_cost, C_rate = 0.5)
        agent_df, model_df = run!(model, dummystep, Model.model_step!, 10; adata = [(:strategy, mean)])
        @test size(agent_df) == (11, 2)
        @test size(model_df) == (0, 0)
    end

    @testset "collect logic" begin
        # make a double-star graph
        G = SimpleGraph(16)
        [add_edge!(G, edge) for edge in [
            (1, 10),
            (1, 2), (1, 3), (1, 4), (1, 5), (1, 6), (1, 7), (1, 8), (1, 9),
            (10, 11), (10, 12), (10, 13), (10, 14), (10, 15), (10, 16)
        ]]
        model = Model.build_model(;G = G, r = 9.0, cost_model = :fixed_cost, C_rate = 1.0)
        model[1].strategy = Agent.D

        @testset "status 1" begin
            @test model[1].strategy == Agent.D
            for id in 2:10
                @test model[id].strategy == Agent.C
            end
        end

        # status 1 -> 2
        Model.model_step!(model)

        @testset "status 2" begin
            for id in 1:9
                @test model[id].strategy == Agent.D
            end
            for id in 10:16
                @test model[id].strategy == Agent.C
            end
        end

        # status 2 -> 3
        Model.model_step!(model)

        @testset "status 3" begin
            @test model[1].strategy == Agent.C
            for id in 2:9
                @test model[id].strategy == Agent.D
            end
            for id in 10:16
                @test model[id].strategy == Agent.C
            end
        end

        # status 3 -> 4
        Model.model_step!(model)

        @testset "status 4" begin
            for id in 1:16
                @test model[id].strategy == Agent.C
            end
        end
    end
end

@testset "exec.calc_r" begin
    G = barabasi_albert(10^3, 2)
    @test r = [round(calc_r(G, η), digits = 1) for η in 0.2:0.1:0.8] == [1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0]
end

@testset "exec.run_simulation" begin
    for cost_model in [:fixed_cost, :variable_cost]
        c_rate = run_simulation(:scale_free, cost_model, 0.7; N_gen = 10^2, N_sim = 3)
        println("[$(cost_model)] c_rate = $(c_rate)")
        @test c_rate isa Float64
    end
end


text:
using FileIO, GeometryBasics
using Test
const tf = joinpath(dirname(@__FILE__), "testfiles")
using MeshIO

function test_face_indices(mesh)
    for face in faces(mesh)
        for index in face
            pass = firstindex(coordinates(mesh)) <= index <= lastindex(coordinates(mesh))
            pass || return false
        end
    end
    return true
end

@testset "MeshIO" begin
    dirlen = 1f0
    baselen = 0.02f0
    mesh = [
        FRect3D(Vec3f0(baselen), Vec3f0(dirlen, baselen, baselen)),
        FRect3D(Vec3f0(baselen), Vec3f0(baselen, dirlen, baselen)),
        FRect3D(Vec3f0(baselen), Vec3f0(baselen, baselen, dirlen))
    ]
    mesh = merge(map(triangle_mesh, mesh))

    mktempdir() do tmpdir
        for ext in ["2dm", "off"]
            @testset "load save $ext" begin
                save(joinpath(tmpdir, "test.$ext"), mesh)
                mesh_loaded = load(joinpath(tmpdir, "test.$ext"))
                @test mesh_loaded == mesh
            end
        end
        @testset "PLY ascii and binary" begin
            f = File(format"PLY_ASCII", joinpath(tmpdir, "test.ply"))
            save(f, mesh)
            mesh_loaded = load(joinpath(tmpdir, "test.ply"))
            @test mesh_loaded == mesh
            save(File(format"PLY_BINARY", joinpath(tmpdir, "test.ply")), mesh)
        end
        @testset "STL ascii and binary" begin
            save(File(format"STL_ASCII", joinpath(tmpdir, "test.stl")), mesh)
            mesh_loaded = load(joinpath(tmpdir, "test.stl"))
            @test Set(mesh.position) == Set(mesh_loaded.position)
            save(File(format"STL_BINARY", joinpath(tmpdir, "test.stl")), mesh)
            mesh_loaded = load(joinpath(tmpdir, "test.stl"))
            @test Set(mesh.position) == Set(mesh_loaded.position)
        end
    end
    @testset "Real world files" begin

        @testset "STL" begin
            msh = load(joinpath(tf, "ascii.stl"))
            @test length(faces(msh)) == 12
            @test length(coordinates(msh)) == 36
            @test length(normals(msh)) == 36
            @test test_face_indices(msh)

            msh = load(joinpath(tf, "binary.stl"))
            @test msh isa GLNormalMesh
            @test length(faces(msh)) == 828
            @test length(coordinates(msh)) == 2484
            @test length(msh.normals) == 2484
            @test test_face_indices(msh)

            mktempdir() do tmpdir
                save(File(format"STL_BINARY", joinpath(tmpdir, "test.stl")), msh)
                msh1 = load(joinpath(tmpdir, "test.stl"))
                @test msh1 isa GLNormalMesh
                @test faces(msh) == faces(msh1)
                @test coordinates(msh) == coordinates(msh1)
                @test msh.normals == msh1.normals
            end

            msh = load(joinpath(tf, "binary_stl_from_solidworks.STL"))
            @test msh isa GLNormalMesh
            @test length(faces(msh)) == 12
            @test length(coordinates(msh)) == 36
            @test test_face_indices(msh)

            # STL Import
            msh = load(joinpath(tf, "cube_binary.stl"))
            @test length(coordinates(msh)) == 36
            @test length(faces(msh)) == 12
            @test test_face_indices(msh)

            msh = load(joinpath(tf, "cube.stl"))
            @test length(coordinates(msh)) == 36
            @test length(faces(msh)) == 12
            @test test_face_indices(msh)

        end
        @testset "PLY" begin
            msh = load(joinpath(tf, "ascii.ply"))
            @test length(faces(msh)) == 36
            @test test_face_indices(msh)
            @test length(coordinates(msh)) == 72
            # msh = load(joinpath(tf, "binary.ply")) # still missing
            #@test msh isa GLNormalMesh
            msh = load(joinpath(tf, "cube.ply")) # quads
            @test length(coordinates(msh)) == 24
            @test length(faces(msh)) == 12
            @test test_face_indices(msh)
        end
        @testset "OFF" begin
            msh = load(joinpath(tf, "test.off"))
            @test length(faces(msh)) == 28
            @test length(coordinates(msh)) == 20
            @test test_face_indices(msh)

            msh = load(joinpath(tf, "test2.off"))
            @test length(faces(msh)) == 810
            @test length(coordinates(msh)) == 405
            @test test_face_indices(msh)

            msh = load(joinpath(tf, "cube.off"))
            @test length(faces(msh)) == 12
            @test length(coordinates(msh)) == 8
            @test test_face_indices(msh)
        end
        @testset "OBJ" begin
            msh = load(joinpath(tf, "test.obj"))
            @test length(faces(msh)) == 3954
            @test length(coordinates(msh)) == 2520
            @test length(normals(msh)) == 2520
            @test test_face_indices(msh)

            msh = load(joinpath(tf, "cube.obj")) # quads
            @test length(faces(msh)) == 12
            @test length(coordinates(msh)) == 8
            @test test_face_indices(msh)

            msh = load(joinpath(tf, "polygonal_face.obj"))
            @test length(faces(msh)) == 4
            @test length(coordinates(msh)) == 6
            @test test_face_indices(msh)

            msh = load(joinpath(tf, "test_face_normal.obj"))
            @test length(faces(msh)) == 1
            @test length(coordinates(msh)) == 3
            @test test_face_indices(msh)
        end
        @testset "2DM" begin
            msh = load(joinpath(tf, "test.2dm"))
            @test test_face_indices(msh)
        end
        @testset "GMSH" begin
            msh = load(joinpath(tf, "cube.msh"))
            @test length(faces(msh)) == 24
            @test length(coordinates(msh)) == 14
            @test test_face_indices(msh)
        end
    end
end


text:
using Test, FileIO
include(joinpath(dirname(@__FILE__), "../../src/includeall.jl"))

rp = joinpath(dirname(@__FILE__), "../reference/SLM")
in_inc = load(joinpath(rp, "find_solution_input.jld2"))
in_dec = load(joinpath(rp, "find_solution_input_decrease.jld2"))
in_sin = load(joinpath(rp, "find_solution_input_sine.jld2"))
in_sinnoineq = load(joinpath(rp, "find_solution_input_sinenoineq.jld2"))
out_inc = load(joinpath(rp, "slm_rescale.jld2"))
out_dec = load(joinpath(rp, "slm_rescale_decrease.jld2"))
out_sin = load(joinpath(rp, "slm_rescale_sine.jld2"))
out_sinnoineq = load(joinpath(rp, "slm_rescale_sinenoineq.jld2"))
cu_int  = [π 2*π; 3*π 4*π]
cd_int  = [0 π; 2*π 3*π]
inc_int = [0 π/2; 3*π/2 5*π/2; 7*π/2 4*pi]
dec_int = [π/2 3*π/2; 5*π/2 7*π/2]
test_inc = load(joinpath(rp, "solve_slm_system.jld2"))
scale_inc = load(joinpath(rp, "scaleproblem.jld2"))

@testset "SLM estimation" begin
    for (i, in_data, out_data) in zip(1:4, [in_inc, in_dec, in_sin, in_sinnoineq], [out_inc, out_dec, out_sin, out_sinnoineq])
        if i == 1
            slm = SLM(vec(in_data["w"]), vec(in_data["p_sol"]); knots = Int(in_data["knots"]),
                      increasing = true, concave_down = true, left_value = in_data["p0_norun"],
                      right_value = in_data["p1_norun"], min_value = in_data["p0_norun"] - 1e-3,
                      max_value = in_data["p1_norun"] + 1e-3)
        elseif i == 2
            slm = SLM(vec(in_data["w"]), vec(in_data["rev_p_sol"]); knots = Int(in_data["knots"]),
                      decreasing = true, concave_up = true, right_value = in_data["p0_norun"],
                      left_value = in_data["p1_norun"], min_value = in_data["p0_norun"] - 1e-3,
                      max_value = in_data["p1_norun"] + 1e-3)
        elseif i == 3
            slm = SLM(vec(in_data["x"]), vec(in_data["y"]); knots = Int(in_data["knots"]),
                      increasing_intervals = inc_int, decreasing_intervals = dec_int,
                      concave_up_intervals = cu_int, concave_down_intervals = cd_int,
                      right_value = 0., left_value = 0.)
        else
            slm = SLM(vec(in_data["x"]), vec(in_data["y"]); knots = Int(in_data["knots"]),
                      right_value = 0., left_value = 0.)
        end
        if i == 1 || i == 2
            @test maximum(abs.(get_coef(slm) - out_data["rescale_coef"])) < 5e-4
        else
            @test @test_matrix_approx_eq get_coef(slm) out_data["rescale_coef"]
        end
    end
end

#=
println("No sparse matrices")
@btime begin
    SLM(vec(in_sin["x"]), vec(in_sin["y"]); knots = Int(in_sin["knots"]),
        increasing_intervals = inc_int, decreasing_intervals = dec_int,
        concave_up_intervals = cu_int, concave_down_intervals = cd_int,
        right_value = 0., left_value = 0.)
end
println("Sparse matrices")
@btime begin
    SLM(vec(in_sin["x"]), vec(in_sin["y"]); knots = Int(in_sin["knots"]),
        increasing_intervals = inc_int, decreasing_intervals = dec_int,
        concave_up_intervals = cu_int, concave_down_intervals = cd_int,
        right_value = 0., left_value = 0., use_sparse = true)
end
=#


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using FMIImport.FMICore: fmi2VariableNamingConventionStructured

myFMU = fmi2Load("SpringFrictionPendulum1D", ENV["EXPORTINGTOOL"], ENV["EXPORTINGVERSION"])
myFMU.executionConfig.assertOnWarning = true

@test fmi2GetVersion(myFMU) == "2.0"
@test fmi2GetTypesPlatform(myFMU) == "default"

@test fmi2GetModelName(myFMU) == "SpringFrictionPendulum1D"
@test fmi2GetVariableNamingConvention(myFMU) == fmi2VariableNamingConventionStructured
@test fmi2IsCoSimulation(myFMU) == true
@test fmi2IsModelExchange(myFMU) == true

@test fmi2GetGUID(myFMU) == "{df491d8d-0598-4495-913e-5b025e54d7f2}"
@test fmi2GetGenerationTool(myFMU) == "Dymola Version 2022x (64-bit), 2021-10-08"
@test fmi2GetGenerationDateAndTime(myFMU) == "2022-03-03T15:09:18Z"
@test fmi2GetNumberOfEventIndicators(myFMU) == 24
@test fmi2CanGetSetState(myFMU) == true
@test fmi2CanSerializeFMUstate(myFMU) == true
@test fmi2ProvidesDirectionalDerivative(myFMU) == true

@test fmi2GetDefaultStartTime(myFMU.modelDescription) ≈ 0.0
@test fmi2GetDefaultStopTime(myFMU.modelDescription) ≈ 1.0
@test fmi2GetDefaultTolerance(myFMU.modelDescription) ≈ 1e-4
@test fmi2GetDefaultStepSize(myFMU.modelDescription) === nothing

fmi2Unload(myFMU)


text:
"""Julia program to implement Fibonacci Search algorithm.
Fibonacci Search is a Divide and Conquer strategy based algorithm that works for sorted arrays
Fibonacci Search avoids CPU costly division operation and uses addition and subtraction opertations
This algorithm uses Fibonacci series as an underlying principle, where every element of the series is defined by:

    fib(n) = fib(n-1) + fib(n-2), where fib(0) = 0, fib(1) = 1
"""

function fibonacci_search(arr, n, ele)

    # We start the series from 1, 1 as the Julia language is 1-based indexing
    fibn2 = 1
    fibn1 = 1
    fibn = fibn1 + fibn2
    while (fibn < n)
        fibn2 = fibn1
        fibn1 = fibn
        fibn = fibn1 + fibn2
    end

    # To denote the eliminted range from left.
    offset = 0

    while (fibn > 1)
        i = min(offset + fibn2, n - 1)

        """Compares the element in i'th index with ele and decide the next fibonacci element. If the element doesn't
            turns out to be greater than or less than the element at i'th index then we find the required element"""
        if (arr[i] < ele)
            fibn = fibn1;
            fibn1 = fibn2;
            fibn2 = fibn - fibn1;
            offset = i;
        elseif (arr[i] > ele)
            fibn = fibn2;
            fibn1 = fibn1 - fibn2;
            fibn2 = fibn - fibn1;
        else
            return true
        end
        if (arr[offset + 1] == ele)
            return true
        end
    end
    return false
end

print("How many numbers are present in the array? ")
n = readline()
n = parse(Int, n)
if (n <= 0)
    println("Array is Empty!!!")
    exit()
end
arr = Int[]
print("Enter the numbers: ")
arr = [parse(Int, num) for num in split(readline())]
print("Which number do you want to search in the array? ")
ele = readline()
ele = parse(Int, ele)
# Sort the array in ascending order
arr = sort(arr)
res = fibonacci_search(arr, n, ele)
if (res == 0)
    print("The number $ele is not present in the array")
else
    print("The number $ele is present in the array.")
end

"""
Time Complexity - O(log(n)), where 'n' is the size of the array
Space Complexity - O(n)

SAMPLE INPUT AND OUTPUT

SAMPLE I

How many numbers are present in the array? 5
Enter the numbers: 1 2 3 4 5
Which number do you want to search in the array? 6
The number 6 is not present in the array

SAMPLE II

How many numbers are present in the array? 3
Enter the numbers: 3 1 2
Which number do you want to search in the array? 2
The number 2 is present in the array.

"""

max_stars_repo_name: TechSpiritSS/NeoAlgo


text:
# This script solves an elliptic partial differential equation

using Printf
using LinearAlgebra

# Define the paths
CURRENT_DIR = @__DIR__
ROOT_DIR = basename(CURRENT_DIR) == "scripts" ? dirname(CURRENT_DIR) : CURRENT_DIR
TABLES = joinpath(ROOT_DIR, "tables")

# Prepare the output directories
mkpath(TABLES)

# Define the problem
l_x = 1
l_y = 1
p(x, y) = 1
q(x, y) = 1
u(x, y) = x^3 * y + y^2 * x
f(x, y) = -(6 * x * y + 2 * x)

# Define the grids
N = [5, 10, 20]
M = copy(N)

# Set the expected precision of the results
ε = 1e-4
digits = 4

"Compute the value of the differential operator L"
function L(U, x, y, h_x, h_y, i, j)::Float64
    return Λ₁(U, x, y, h_x, i, j) + Λ₂(U, x, y, h_y, i, j)
end

"Compute the value of the differential operator Λ₁"
function Λ₁(U, x, y, h_x, i, j)::Float64
    return p(x[i] + h_x / 2, y[j]) * (U[i+1, j] - U[i, j]) / h_x^2 -
           p(x[i] - h_x / 2, y[j]) * (U[i, j] - U[i-1, j]) / h_x^2
end

"Compute the value of the differential operator Λ₂"
function Λ₂(U, x, y, h_y, i, j)::Float64
    return q(x[i], y[j] + h_y / 2) * (U[i, j+1] - U[i, j]) / h_y^2 -
           q(x[i], y[j] - h_y / 2) * (U[i, j] - U[i, j-1]) / h_y^2
end

"Compute ξ, which is needed to compute an optimal number of iterations"
function compute_xi(l_x, l_y, h_x, h_y)::Float64
    rx = 0:0.001:l_x
    ry = 0:0.001:l_y
    pairs = [(x, y) for x in rx, y in ry]
    c₁ = minimum(args -> p(args...), pairs)
    c₂ = maximum(args -> p(args...), pairs)
    d₁ = minimum(args -> q(args...), pairs)
    d₂ = maximum(args -> q(args...), pairs)
    k₁ = 4 / h_x^2
    k₂ = 4 / h_y^2
    arg₁ = (π * h_x) / (2 * l_x)
    arg₂ = (π * h_y) / (2 * l_y)
    δ = c₁ * k₁ * sin(arg₁)^2 +
        d₁ * k₂ * sin(arg₂)^2
    Δ = c₂ * k₁ * cos(arg₁)^2 +
        d₂ * k₂ * cos(arg₂)^2
    return δ / Δ
end

"Initialize a new matrix and compute the boundary values"
function boundary_values(N, M, x, y)::Matrix{Float64}
    # Prepare a matrix for the solution
    U = zeros(N + 1, M + 1)
    # Compute the boundary values
    U[:, 1] .= u.(x, 0)
    U[:, M+1] .= u.(x, l_y)
    U[1, 2:M] .= u.(0, y[2:M])
    U[N+1, 2:M] .= u.(l_x, y[2:M])
    return U
end

"Iterate using the simple iteration method"
function simple_iteration(Uₖ₋₁, ξ, ε, N, M, h_x, h_y, x, y)::Tuple{Matrix{Float64},Int}
    # Prepare a new matrix
    Uₖ = copy(Uₖ₋₁)
    # Compute the optimal number of iterations
    m = ceil(Int, log(1 / ε) / (2ξ))
    # Iterate enough times to achieve desired precision
    for _ in 1:m
        # Compute values at the inner nodes
        for i = 2:N, j = 2:M
            k₁ = p(x[i] - h_x / 2, y[j]) / h_x^2
            k₂ = p(x[i] + h_x / 2, y[j]) / h_x^2
            k₃ = q(x[i], y[j] - h_y / 2) / h_y^2
            k₄ = q(x[i], y[j] + h_y / 2) / h_y^2
            Uₖ[i, j] = (k₁ * Uₖ₋₁[i-1, j] +
                        k₂ * Uₖ₋₁[i+1, j] +
                        k₃ * Uₖ₋₁[i, j-1] +
                        k₄ * Uₖ₋₁[i, j+1] +
                        f(x[i], y[j])) /
                       (k₁ + k₂ + k₃ + k₄)
        end
        # Reassign the current iteration as the previous one
        Uₖ₋₁ = Uₖ
    end
    return Uₖ, m
end

"Iterate using the Seidel's method"
function seidel(Uₖ₋₁, ξ, ε, N, M, h_x, h_y, x, y)::Tuple{Matrix{Float64},Int}
    # Prepare a new matrix
    Uₖ = copy(Uₖ₋₁)
    # Compute the optimal number of iterations
    m = ceil(Int, log(1 / ε) / (4ξ))
    # Iterate enough times to achieve desired precision
    for _ in 1:m
        # Compute values at the inner nodes
        for i = N:-1:2, j = M:-1:2
            k₁ = p(x[i] - h_x / 2, y[j]) / h_x^2
            k₂ = p(x[i] + h_x / 2, y[j]) / h_x^2
            k₃ = q(x[i], y[j] - h_y / 2) / h_y^2
            k₄ = q(x[i], y[j] + h_y / 2) / h_y^2
            Uₖ[i, j] = (k₁ * Uₖ[i-1, j] +
                        k₂ * Uₖ₋₁[i+1, j] +
                        k₃ * Uₖ[i, j-1] +
                        k₄ * Uₖ₋₁[i, j+1] +
                        f(x[i], y[j])) /
                       (k₁ + k₂ + k₃ + k₄)
        end
        # Reassign the current iteration as the previous one
        Uₖ₋₁ = Uₖ
    end
    return Uₖ, m
end

"Iterate using the upper relaxation method"
function upper_relaxation(Uₖ₋₁, ξ, ε, N, M, h_x, h_y, x, y)::Tuple{Matrix{Float64},Int}
    # Prepare a new matrix
    Uₖ = copy(Uₖ₋₁)
    # Compute the optimal number of iterations
    m = ceil(Int, log(1 / ε) / sqrt(ξ))
    # Define ω in (0,2), which affects the speed of convergence
    ω = 1.0
    # Iterate enough times to achieve desired precision
    for _ in 1:m
        # Compute values at the inner nodes
        for i = N:-1:2, j = M:-1:2
            k₁ = p(x[i] + h_x / 2, y[j]) / h_x^2
            k₂ = p(x[i] - h_x / 2, y[j]) / h_x^2
            k₃ = q(x[i], y[j] + h_y / 2) / h_y^2
            k₄ = q(x[i], y[j] - h_y / 2) / h_y^2
            Uₖ[i, j] = Uₖ₋₁[i, j] +
                       ω * (k₁ * (Uₖ₋₁[i+1, j] - Uₖ₋₁[i, j]) -
                            k₂ * (Uₖ₋₁[i, j] - Uₖ[i-1, j]) +
                            k₃ * (Uₖ₋₁[i, j+1] - Uₖ₋₁[i, j]) -
                            k₄ * (Uₖ₋₁[i, j] - Uₖ[i, j-1]) +
                            f(x[i], y[j])) /
                       (k₁ + k₂ + k₃ + k₄)
        end
        # Reassign the current iteration as the previous one
        Uₖ₋₁ = Uₖ
    end
    return Uₖ, m
end

"Iterate, using the alternately triangular iterative method"
function triangular(Uₖ₋₁, ε, N, M, h_x, h_y, x, y)::Tuple{Matrix{Float64},Int}
    # Prepare a new matrix
    Uₖ = copy(Uₖ₋₁)
    # Compute the coefficients
    pairs = [(x, y) for x in 0:0.001:l_x, y in 0:0.001:l_y]
    c₁ = minimum(args -> p(args...), pairs)
    c₂ = maximum(args -> p(args...), pairs)
    d₁ = minimum(args -> q(args...), pairs)
    d₂ = maximum(args -> q(args...), pairs)
    k₁ = 4 / h_x^2
    k₂ = 4 / h_y^2
    δ = c₁ * k₁ * sin((π * h_x) / (2 * l_x))^2 +
        d₁ * k₂ * sin((π * h_y) / (2 * l_y))^2
    Δ = c₂ * k₁ + d₂ * k₂
    ω = 2 / sqrt(δ * Δ)
    η = δ / Δ
    γ₁ = δ / (2 + 2 * sqrt(η))
    γ₂ = δ / (4 * sqrt(η))
    ξ = γ₁ / γ₂
    κ₁ = ω / h_x^2
    κ₂ = ω / h_y^2
    τ = 2 / (γ₁ + γ₂)
    # Compute the optimal number of iterations
    m = ceil(Int, log(1 / ε) / log((1 + ξ) / (1 - ξ)))
    # Iterate enough times to achieve desired precision
    for _ in 1:m
        # Prepare intermediate matrices
        Ũ = zeros(size(Uₖ₋₁)...)
        U̅ = zeros(size(Uₖ₋₁)...)
        # Compute values at the inner nodes of the first intermediate matrix
        for i = 2:N, j = 2:M
            k₁ = κ₁ * p(x[i] - h_x / 2, y[j])
            k₂ = κ₂ * q(x[i], y[j] - h_y / 2)
            Ũ[i, j] = (k₁ * Ũ[i-1, j] + k₂ * Ũ[i, j-1] +
                        L(Uₖ₋₁, x, y, h_x, h_y, i, j) + f(x[i], y[j])) /
                       (1 + k₁ + k₂)
        end
        # Compute values at the inner nodes of the second intermediate matrix
        for i = N:-1:2, j = M:-1:2
            k₁ = κ₁ * p(x[i] + h_x / 2, y[j])
            k₂ = κ₂ * q(x[i], y[j] + h_y / 2)
            U̅[i, j] = (k₁ * U̅[i+1, j] + k₂ * U̅[i, j+1] + Ũ[i, j]) /
                       (1 + k₁ + k₂)
        end
        # Add up the matrices
        Uₖ[2:N, 2:M] .= Uₖ₋₁[2:N, 2:M] .+ τ .* U̅[2:N, 2:M]
        # Reassign the current iteration as the previous one
        Uₖ₋₁ = Uₖ
    end
    return Uₖ, m
end

"Iterate using the variable directions method"
function variable_directions(Uₖ₋₁, ε, N, M, h_x, h_y, x, y)::Tuple{Matrix{Float64},Int}
    # Prepare a new matrix and an intermediate matrix
    Uₖ = copy(Uₖ₋₁)
    U̅ = copy(Uₖ₋₁)
    # Compute the optimal number of iterations
    m = ceil(Int, N / (2 * π) * log(1 / ε))
    # Compute the coefficients
    pairs = [(x, y) for x in 0:0.001:l_x, y in 0:0.001:l_y]
    c₁ = minimum(args -> p(args...), pairs)
    c₂ = maximum(args -> p(args...), pairs)
    d₁ = minimum(args -> q(args...), pairs)
    d₂ = maximum(args -> q(args...), pairs)
    k₁ = 4 / h_x^2
    k₂ = 4 / h_y^2
    arg₁ = (π * h_x) / (2 * l_x)
    arg₂ = (π * h_y) / (2 * l_y)
    δ₁ = c₁ * k₁ * sin(arg₁)^2
    δ₂ = d₁ * k₂ * sin(arg₂)^2
    Δ₁ = c₂ * k₁ * cos(arg₁)^2
    Δ₂ = d₂ * k₁ * cos(arg₂)^2
    δ = min(δ₁, δ₂)
    Δ = max(Δ₁, Δ₂)
    τ = 2 / sqrt(δ * Δ)
    # Iterate enough times to achieve desired precision
    for _ in 1:m
        # Compute values at the inner nodes of the intermediate matrix
        for j = 2:M
            # Compute the linear system's matrix
            Ũ = Tridiagonal(
                # A's
                [repeat([τ / (2 * h_x^2)], N - 1); 0],
                # B's
                [1; repeat([-τ / h_x^2 + 1], N - 1); 1],
                # C's
                [0; repeat([τ / (2 * h_x^2)], N - 1)],
            )
            # Compute the linear system's right-hand side vector
            g = [
                Uₖ₋₁[1, j]
                [-Uₖ₋₁[i, j] -
                 τ / 2 * (Λ₂(Uₖ₋₁, x, y, h_y, i, j) +
                          f(x[i], y[j]))
                 for i in 2:N
                ]
                Uₖ₋₁[N+1, j]]
            # Compute the solution of the linear system
            U̅[:, j] .= Ũ \ g
        end
        # Compute values at the inner nodes of the input matrix
        for i = 2:N
            # Compute the linear system's matrix
            Ũ = Tridiagonal(
                # A's
                [repeat([τ / (2 * h_y^2)], N - 1); 0],
                # B's
                [1; repeat([-τ / h_y^2 + 1], N - 1); 1],
                # C's
                [0; repeat([τ / (2 * h_y^2)], N - 1)],
            )
            # Compute the linear system's right-hand side vector
            g = [
                Uₖ₋₁[i, 1]
                [-U̅[i, j] -
                 τ / 2 * (Λ₁(U̅, x, y, h_x, i, j) +
                          f(x[i], y[j]))
                 for j in 2:M
                ]
                Uₖ₋₁[i, M+1]]
            # Compute the solution of the linear system
            Uₖ[i, :] .= Ũ \ g
        end
        # Reassign the current iteration as the previous one
        Uₖ₋₁ = Uₖ
    end
    return Uₖ, m
end

"Create the TeX tables and write them to disk"
function tables(dir, U, N, M; digits = 8)
    # Prepare the output directory
    OUTPUT_DIR = joinpath(TABLES, dir)
    mkpath(OUTPUT_DIR)
    # Prepare the ranges
    ri = Int.([1:N/5:N+1]...)
    rj = Int.([1:M/5:M+1]...)
    # Create and write the table with the solution grid
    open(joinpath(OUTPUT_DIR, "$N, $M.tex"), "w") do io
        for i in ri
            s = ""
            for j in rj
                number = round(U[i, j]; digits)
                s = "$(s)& \$ $(Printf.format(Printf.Format("%.0$(digits)f"), number == 0 ? 0 : number)) \$ "
            end
            s = "$(s)\\\\"
            println(io, s)
        end
    end
end

"""
Compare the approximate solution to the exact one,
find the biggest absolute difference between the nodes
"""
function max_difference(U, x, y)::Float64
    M1, N1 = size(U)
    # Find the maximum
    Δu = 0
    for i = 1:N1, j = 1:M1
        _Δu = abs(u(x[i], y[j]) - U[i, j])
        if _Δu > Δu
            Δu = _Δu
        end
    end
    return Δu
end

# For each grid
for l in eachindex(N)
    # Compute the steps
    h_x = l_x / N[l]
    h_y = l_y / M[l]
    # Print the numbers of nodes
    println(
        '\n',
        " "^5, "N: ", N[l], '\n',
        " "^5, "M: ", M[l]
    )
    # Compute the nodes
    x = collect(0:h_x:1.0)
    y = collect(0:h_y:1.0)
    # Compute ξ
    ξ = compute_xi(l_x, l_y, h_x, h_y)
    # Create a new matrix and compute the boundary values
    U = boundary_values(N[l], M[l], x, y)
    # Iterate, using the simple iteration method
    U, m = simple_iteration(U, ξ, ε, N[l], M[l], h_x, h_y, x, y)
    # Create and write the TeX tables
    tables("simple_iteration", U, N[l], M[l]; digits)
    # Print the info
    println(
        '\n',
        " "^5, "> Simple iteration:", '\n',
        " "^5, "  m = ", m, '\n',
        " "^5, "  Δu = ", max_difference(U, x, y)
    )
    # Iterate, using the Seidel's method
    U, m = seidel(U, ξ, ε, N[l], M[l], h_x, h_y, x, y)
    # Create and write the TeX tables
    tables("seidel", U, N[l], M[l]; digits)
    # Print the info
    println(
        '\n',
        " "^5, "> Seidel:", '\n',
        " "^5, "  m = ", m, '\n',
        " "^5, "  Δu = ", max_difference(U, x, y)
    )
    # Iterate, using the upper relaxation method
    U, m = upper_relaxation(U, ξ, ε, N[l], M[l], h_x, h_y, x, y)
    # Create and write the TeX tables
    tables("upper_relaxation", U, N[l], M[l]; digits)
    # Print the info
    println(
        '\n',
        " "^5, "> Upper relaxation:", '\n',
        " "^5, "  m = ", m, '\n',
        " "^5, "  Δu = ", max_difference(U, x, y)
    )
    # Iterate, using the alternately triangular iterative method
    U, m = triangular(U, ε, N[l], M[l], h_x, h_y, x, y)
    # Create and write the TeX tables
    tables("triangular", U, N[l], M[l]; digits)
    # Print the info
    println(
        '\n',
        " "^5, "> Alternately triangular iterative method:", '\n',
        " "^5, "  m = ", m, '\n',
        " "^5, "  Δu = ", max_difference(U, x, y)
    )
    # Iterate, using the variable directions method
    U, m = variable_directions(U, ε, N[l], M[l], h_x, h_y, x, y)
    # Create and write the TeX tables
    tables("variable_directions", U, N[l], M[l]; digits)
    # Print the info
    println(
        '\n',
        " "^5, "> Variable directions method:", '\n',
        " "^5, "  m = ", m, '\n',
        " "^5, "  Δu = ", max_difference(U, x, y)
    )
end

println()


text:
const ContactImplicitMPC = Main
include(joinpath(@__DIR__, "..", "src/dynamics", "quadruped", "visuals.jl"))
T = Float64
vis = Visualizer()
open(vis)
render(vis)

s_sim = get_simulation("quadruped", "sine1_2D_lc", "sinusoidal")
s = get_simulation("quadruped", "flat_2D_lc", "flat")
model = s.model
env = s.env

nq = model.dim.q

s_sim.model.μ_world
s.model.μ_world

ref_traj_ = deepcopy(ContactImplicitMPC.get_trajectory(s.model, s.env,
    joinpath(module_dir(), "src/dynamics/quadruped/gaits/gait2.jld2"),
    load_type = :split_traj_alt))
ref_traj = deepcopy(ref_traj_)

# time
H = ref_traj.H
h = ref_traj.h
N_sample = 5
H_mpc = 10
h_sim = h / N_sample
H_sim = 5000 #4000 #5000

# barrier parameter
κ_mpc = 1.0e-4

obj = TrackingObjective(model, env, H_mpc,
    q = [Diagonal(1e-2 * [10; 0.02; 0.25; 0.25 * ones(nq-3)]) for t = 1:H_mpc],
    u = [Diagonal(3e-2 * ones(model.dim.u)) for t = 1:H_mpc],
    γ = [Diagonal(1.0e-100 * ones(model.dim.c)) for t = 1:H_mpc],
    b = [Diagonal(1.0e-100 * ones(model.dim.c * friction_dim(env))) for t = 1:H_mpc])

p = linearized_mpc_policy(ref_traj, s, obj,
    H_mpc = H_mpc,
    N_sample = N_sample,
    κ_mpc = κ_mpc,
    n_opts = NewtonOptions(
        r_tol = 3e-4,
        max_iter = 5,
        # verbose = true,
        ),
    mpc_opts = LinearizedMPCOptions(
        # live_plotting=true,
        altitude_update = true,
        altitude_impact_threshold = 0.05,
        altitude_verbose = true,
        )
    )

q1_ref = copy(ref_traj.q[2])
q0_ref = copy(ref_traj.q[1])
q1_sim = SVector{model.dim.q}(q1_ref)
q0_sim = SVector{model.dim.q}(copy(q1_sim - (q1_ref - q0_ref) / N_sample))
@assert norm((q1_sim - q0_sim) / h_sim - (q1_ref - q0_ref) / h) < 1.0e-8

w_amp = [+0.02, -0.20]
sim = simulator(s_sim, q0_sim, q1_sim, h_sim, H_sim,
    p = p,
    # d = open_loop_disturbances([rand(model.dim.w) .* w_amp for i=1:H_sim]),
    ip_opts = InteriorPointOptions(
        r_tol = 1.0e-8,
        κ_init = 1.0e-6,
        κ_tol = 2.0e-6),
    sim_opts = SimulatorOptions(warmstart = true)
    )

@time status = simulate!(sim)

# plt = plot(layout=(3,1), legend=false)
# plot!(plt[1,1], hcat(Vector.(vcat([fill(ref_traj.q[i], N_sample) for i=1:H]...))...)',
#     color=:red, linewidth=3.0)
# plot!(plt[1,1], hcat(Vector.(sim.traj.q)...)', color=:blue, linewidth=1.0)
# plot!(plt[2,1], hcat(Vector.(vcat([fill(ref_traj.u[i][1:nu], N_sample) for i=1:H]...))...)',
#     color=:red, linewidth=3.0)
# plot!(plt[2,1], hcat(Vector.([u[1:nu] for u in sim.traj.u]*N_sample)...)', color=:blue, linewidth=1.0)
# plot!(plt[3,1], hcat(Vector.([γ[1:nc] for γ in sim.traj.γ]*N_sample)...)', color=:blue, linewidth=1.0)

plot_surface!(vis, s_sim.env, ylims=[0.3, -0.05])
plot_lines!(vis, model, sim.traj.q[1:5:end])
anim = visualize_meshrobot!(vis, model, sim.traj, sample=5)
# anim = visualize_robot!(vis, model, sim.traj, anim=anim)
anim = visualize_force!(vis, model, sim.traj, anim=anim, h=h_sim)

# Display ghosts
t_ghosts = [1, 1333, 2666]
mvis_ghosts = []
for (i,t) in enumerate(t_ghosts)
    α = i/(length(t_ghosts)+1)
    name = Symbol("ghost$i")
    mvis = build_meshrobot!(vis, model, name=name, α=α)
    push!(mvis_ghosts, mvis)
end

for (i,t) in enumerate(t_ghosts)
    name = Symbol("ghost$i")
    set_meshrobot!(vis, mvis_ghosts[i], model, sim.traj.q[t], name=name)
end