text:
##########################
# Approximation problems
##########################

using BasisFunctions: Dictionary, TensorProductDict
using DomainSets: Domain

"""
    abstract type ApproximationProblem end

Approximation problem types group the arguments that a user supplies
to the `Fun` constructor. This may be a dictionary or a platform, for example.
By passing around approximation problem objects, functions can easily
implement the same interface as `Fun`.

All functionality is eventually delegated back to the underlying dictionary or
platform. Thus, an approximation problem is nothing but an empty intermediate
layer that passes on the questions to the right object regardless of how the
question was asked.

Some of the properties of an approximation problem may be cached in order to
improve efficiency.
"""
abstract type ApproximationProblem end

# Caching of properties
cache(ap::ApproximationProblem) = ap.cache
cache(ap::ApproximationProblem, property) = ap.cache[property]
cache!(ap::ApproximationProblem, property, value) = ap.cache[property] = value

ap_hasproperty(ap::ApproximationProblem, property) = haskey(cache(ap), property)

ap_fun(ap::ApproximationProblem) = cache(ap, ap_fun)
ap_fun!(ap::ApproximationProblem, f) = cache!(ap, ap_fun, f)

# For convenience
measure(ap::ApproximationProblem) = measure(platform(ap))
coefficienttype(ap::ApproximationProblem) = coefficienttype(platform(ap))

approximation_measure(ap::ApproximationProblem) = approximation_measure(SamplingStyle(ap), ap)
approximation_measure(ss::DiscreteStyle, ap::ApproximationProblem) = discretemeasure(ap)
approximation_measure(ss::ProjectionStyle, ap::ApproximationProblem) = measure(ap)
approximation_measure(ss::ProductSamplingStyle, ap::ApproximationProblem) = approximation_measure(factor(ss,1),ap)

"A `PlatformApproximation` stores a platform and a platform parameter."
mutable struct PlatformApproximation <: ApproximationProblem
    platform    ::  Platform
    plt_par         # platform parameter
    cache

    PlatformApproximation(platform, plt_par) = new(platform, plt_par, Dict{Any,Any}())
    function PlatformApproximation(platform, plt_par, smpl_par)
        cache = Dict{Any,Any}()
        cache[samplingparameter] = smpl_par
        new(platform, plt_par, cache)
    end
end

platform(ap::PlatformApproximation) = ap.platform

export platformparameter
"Return the platform parameter of the approximation problem."
platformparameter(ap::PlatformApproximation) = ap.plt_par
@deprecate parameter(ap::ApproximationProblem) platformparameter(ap)

function dictionary(ap::PlatformApproximation)
    if ap_hasproperty(ap, dictionary)
        cache(ap, dictionary)
    else
        Φ = dictionary(platform(ap), platformparameter(ap))
        cache!(ap, dictionary, Φ)
    end
end

function dualdictionary(ap::PlatformApproximation)
    if ap_hasproperty(ap, dualdictionary)
        cache(ap, dualdictionary)
    else
        Φ = dualdictionary(platform(ap), platformparameter(ap), approximation_measure(ap))
        cache!(ap, dualdictionary, Φ)
    end
end

coefficienttype(ap::PlatformApproximation) = coefficienttype(dictionary(ap))

# specialized functions for when the platform is a product platform

productparameter(ap::PlatformApproximation) = productparameter(platform(ap), platformparameter(ap))

function factors(ap::PlatformApproximation)
    if ap_hasproperty(ap, samplingparameter)
        smpl_par = samplingparameter(ap)
        if length(smpl_par)==length(productparameter(ap))
            map(approximationproblem, factors(platform(ap)), productparameter(ap), smpl_par)
        else
            error("sampling parameter should contain $(length(productparameter(ap))) elements")
        end
    else
        map(approximationproblem, components(platform(ap)), productparameter(ap))
    end
end

unsafe_factors(ap::PlatformApproximation) =
    map(approximationproblem, components(platform(ap)), productparameter(ap))

components(ap::PlatformApproximation) = factors(ap)

"""
    struct AdaptiveApproximation <: ApproximationProblem

An `AdaptiveApproximation` only stores a platform, with which to compute adaptive
approximations. It has no single platform parameter.
"""
struct AdaptiveApproximation <: ApproximationProblem
    platform    ::  Platform
    cache
    AdaptiveApproximation(platform::Platform) = new(platform, Dict{Any,Any}())
end

platform(ap::AdaptiveApproximation) = ap.platform

dictionary(ap::AdaptiveApproximation) =
    error("An adaptive approximation problem does not have a single dictionary.")
dictionary(ap::AdaptiveApproximation, plt_par) = dictionary(platform(ap), plt_par)

platformparameter(ap::AdaptiveApproximation) =
    error("An adaptive approximation problem does not contain a platform parameter.")

## Now on to the actual interface:

to_platform(dict::Dictionary) = (platform(dict), platformparameter(dict))
to_platform(dict::Dictionary, smpl_par) =
    (platform(dict), platformparameter(dict), smpl_par)
function to_platform(dict::Dictionary, domain::Domain, args...)
    if domain == support(dict)
        to_platform(dict, args...)
    else
        to_platform(extensionframe(domain, dict), args...)
    end
end

export approximationproblem
"""
    approximationproblem(args...)

# examples

## Create an approximationproblem from a dictionary
```jldocs
julia> approximationproblem(Fourier(10))
```

## Create an approximationproblem from a dictionary and a domain (extensionframe)
```jldocs
julia> approximationproblem(Fourier(10),0.0..0.5)
```

## Create an approximationproblem from a platform
```jldocs
julia> approximationproblem(FourierPlatform())
```

## Create an approximationproblem from a platform and a platform parameter
```jldocs
julia> approximationproblem(FourierPlatform(), 10)
```
"""
approximationproblem(dict::Dictionary, args...) =
    approximationproblem(to_platform(dict, args...)...)

# 2. An approximation problem from a platform and a concrete parameter value
function approximationproblem(platform::Platform, plt_par)
    @assert correctparamformat(platform, plt_par)
    PlatformApproximation(platform, plt_par)
end
function approximationproblem(platform::Platform, plt_par, smpl_par)
    @assert correctparamformat(platform, plt_par)
    PlatformApproximation(platform, plt_par, smpl_par)
end

# 3. An adaptive approximation problem from a platform
approximationproblem(platform::Platform) = AdaptiveApproximation(platform)

function copy_cache(src, dest)
    for k in keys(cache(src))
        cache!(dest, k, cache(src, k))
    end
end

# From the adaptive approximation problem we can create a concrete platform
# approximation by specifying a parameter value.
function approximationproblem(ap::AdaptiveApproximation, plt_par)
    ap2 = approximationproblem(platform(ap), plt_par)
    copy_cache(ap, ap2)
    ap2
end
function approximationproblem(ap::AdaptiveApproximation, plt_par, smpl_par)
    ap2 = approximationproblem(platform(ap), plt_par, smpl_par)
    copy_cache(ap, ap2)
    ap2
end

# 4. Enable options
function approximationproblem(args...; options...)
    ap = approximationproblem(args...)
    ap_process_options(ap; options...)
    ap
end

# 5. Allow a function to be specified as first argument
approximationproblem(fun, dict::Dictionary, args...) =
    approximationproblem(fun, approximationproblem(dict, args...))
approximationproblem(fun, platform::Platform, args...) =
    approximationproblem(fun, approximationproblem(platform, args...))
function approximationproblem(fun, ap::ApproximationProblem)
    ap_fun!(ap, fun)
    ap
end

export _ap
const _ap = approximationproblem


text:
"""
    o_mincut(graph, values; interdiction = 0)
Compute the capacity of a cut (determined by the state of the solver) with a possible `interdiction` on the highest capacited links.
"""
function o_mincut(graph, values; interdiction = 0)
    capacity = Vector{Float64}()
    n = size(graph, 1)
    for i in 1:n, j in 1:n
        (values[i] < values[n + 1] < values[j]) && push!(capacity, graph[i, j])
    end
    return sum(sort!(capacity)[1:(end - interdiction)])
end

ext: jl


text:
function farloi01()
#
# Farmer-Loizou
#
    p = conv([1.0;1;2],[1;1;3]);
    p = conv(p,p); p = conv(p,p);
    z = [-0.50000000000000 + 1.65831239517770im
        -0.50000000000000 - 1.65831239517770im
        -0.50000000000000 + 1.32287565553230im
        -0.50000000000000 - 1.32287565553230im];
    z = [z 4*ones(4)];

    p, PolyZeros(z)
end


text:
# Geometry which displays points at given (x, y) positions.
immutable PointGeometry <: Gadfly.GeometryElement
end

const point = PointGeometry

function element_aesthetics(::PointGeometry)
    [:x, :y, :size, :color]
end

# Generate a form for a point geometry.
#
# Args:
#   geom: point geometry.
#   theme: the plot's theme.
#   aes: aesthetics.
#
# Returns:
#   A compose Form.
#
function render(geom::PointGeometry, theme::Gadfly.Theme,
                aes::Gadfly.Aesthetics, data::Gadfly.Data,
                scales::Dict{Symbol, ScaleElement})
    Gadfly.assert_aesthetics_defined("Geom.point", aes, :x, :y)
    Gadfly.assert_aesthetics_equal_length("Geom.point", aes,
                                          element_aesthetics(geom)...)

    default_aes = Gadfly.Aesthetics()
    default_aes.color = PooledDataArray(RGB{Float32}[theme.default_color])
    default_aes.size = Measure[theme.default_point_size]
    aes = inherit(aes, default_aes)

    lw_hover_scale = 10
    lw_ratio = theme.line_width / aes.size[1]

    if !all(map(isconcrete, aes.x)) || !all(map(isconcrete, aes.y))
        aes_x = Array(eltype(aes.x), 0)
        aes_y = Array(eltype(aes.y), 0)
        for (x, y) in zip(aes.x, aes.y)
            if !isconcrete(x) || !isconcrete(y)
                continue
            end
            push!(aes_x, x)
            push!(aes_y, y)
        end
    else
        aes_x, aes_y = aes.x, aes.y
    end

    ctx = compose!(
        context(),
        circle(aes.x, aes.y, aes.size),
        fill(aes.color),
        linewidth(theme.highlight_width))

    if aes.color_key_continuous != nothing && aes.color_key_continuous
        compose!(ctx,
            stroke(map(theme.continuous_highlight_color, aes.color)))
    else
        compose!(ctx,
            stroke(map(theme.discrete_highlight_color, aes.color)),
            svgclass([svg_color_class_from_label(escape_id(aes.color_label([c])[1]))
                      for c in aes.color]))
    end

    return compose!(context(order=4), svgclass("geometry"), ctx)
end


text:
using Term
import Term: chars

prt(pan) = begin
    print(" " * hLine(pan.measure.w; style="red"))
    print(vLine(pan.measure.h; style="red") * pan)
    println(pan.measure, "   ", length(pan.segments) )
end

print("\n"^5)

# check that unreasonable widhts are ignored
@time prt(TextBox(
    "nofit"^25;
    width=1000
))

@time prt(TextBox(
    "truncate"^25;
    width=100,
    fit=:truncate
))

@time prt(TextBox(
    "truncate"^25;
    width=100,
))

@time prt(TextBox(
    "truncate"^8;
    fit=:fit
))

@time prt(TextBox(
    "[red]truncate[/red]"^8;
    fit=:fit
))

@time prt(TextBox(
    "[red]tru\nncate[/red]test"^1;
    fit=:fit
))

@time prt(
    TextBox(
    "[red]truncate[/red]test"^8;
    fit=:fit,
    justify=:left
)
)


text:
module EEGTest

# package code goes here
export AnalogData, DigitalData, Spectrogram

include("eegtypes.jl")

export ad_equals, normalize_data, threshold_01, lowpass, highpass,
down_sample, toDecibels, debounce_discrete_signal, dd_equals,
truncate_by_index, truncate_by_value

include("eegutilities.jl")

end # module


text:
using Test
using LinearAlgebra
using WavePropBase
using WavePropBase.Geometry
using WavePropBase.Integration

@testset "Trapezoidal quadrature" begin
    q       = Trapezoidal{10}()
    D       = domain(q)
    @test D == ReferenceLine()
    x,w = q()
    @test sum(w) ≈ 1
    # integrate a periodic function. Should be very accurate.
    @test isapprox(integrate(x->cos(2π*x),q),0,atol=1e-10)
    @test integrate(x->sin(2π*x[1])^2,q) ≈ 0.5
end

@testset "TrapezoidalOpen quadrature" begin
    q = TrapezoidalOpen{10}()
    D = domain(q)
    @test D == ReferenceLine()
    x,w = q()
    @test sum(w) ≈ 1
    # integrate a periodic function. Should be very accurate.
    @test isapprox(integrate(x->cos(2π*x[1]),q),0,atol=1e-10)
    @test integrate(x->sin(2π*x[1])^2,q) ≈ 0.5
end

@testset "Fejer quadrature" begin
    N = 5
    q = Fejer{N}()
    x,w = q()
    @test sum(w) ≈ 1
    # integrate all polynomial of degree N-1 exactly
    for n in 1:(N-1)
        @test integrate(x->x[1]^n,q) ≈ 1/(n+1)
    end
end

@testset "Gauss quad on triangle" begin
    d = ReferenceTriangle()
    # exact value for x^a*y^b integrate over reference triangle
    exa = (a,b) -> factorial(a)*factorial(b)/factorial(a+b+2)
    # check all quadrature implemented
    orders = keys(Integration.TRIANGLE_GAUSS_ORDER_TO_NPTS)
    for p in orders
        q     = Gauss(;domain=d,order=p)
        x,w   = q()
        for i in 0:p
            for j in 0:p-i
                @test integrate(x->x[1]^i*x[2]^j,q) ≈ exa(i,j)
            end
        end
    end
end

@testset "Gauss quad on tetrahedron" begin
    d = ReferenceTetrahedron()
    for p in (1,2)
        q = Gauss(;domain=d,order=p)
        x,w = q()
        @test sum(w) ≈ 1/6
    end
    # FIXME: check that we integrate all monomials up to `order` like in the
    # reference triangle
end

@testset "Tensor product quad on square" begin
    N,M = 10,12
    qx  = Fejer(N)
    qy  = Fejer(M)
    q   = TensorProductQuadrature(qx,qy)
    a,b = N-1,M-1 # maximum integration order of monomials
    @test integrate( x -> 1,q) ≈ 1
    f = x -> x[1]^a*x[2]^b
    @test integrate(f,q) ≈ 1/(a+1)*1/(b+1)
end

max_stars_repo_name: WaveProp/WavePropBase

max_issues_repo_name: WaveProp/WavePropBase

max_forks_repo_name: WaveProp/WavePropBase


text:
module ConjugateGradientSolver

export ConjugateGradient

using ..LinearSolvers
const LS = LinearSolvers
using ClimateMachine.MPIStateArrays
using LinearAlgebra
using LazyArrays
using StaticArrays

struct ConjugateGradient{AT1, AT2, FT, RD, RT, IT} <:
       LS.AbstractIterativeLinearSolver
    # tolerances (2)
    rtol::FT
    atol::FT
    # arrays of size reshape_tuple (7)
    r0::AT1
    z0::AT1
    p0::AT1
    r1::AT1
    z1::AT1
    p1::AT1
    Lp::AT1
    # arrays of size(MPIStateArray) which are aliased to two of the previous dimensions (2)
    alias_p0::AT2
    alias_Lp::AT2
    # reduction dimension (1)
    dims::RD
    # reshape dimension (1)
    reshape_tuple::RT
    # maximum number of iterations (1)
    max_iter::IT

end

# Define the outer constructor for the ConjugateGradient struct
"""
# ConjugateGradient
function ConjugateGradient(Q::AT; rtol = eps(eltype(Q)), atol = eps(eltype(Q)), dims = :) where AT

# Description
- Outer constructor for the ConjugateGradient struct

# Arguments
- `Q`:(array). The kind of object that linearoperator! acts on.

# Keyword Arguments
- `rtol`: (float). relative tolerance
- `atol`: (float). absolute tolerance
- `dims`: (tuple or : ). the dimensions to reduce over
- `reshape_tuple`: (tuple). the dimensions that the conjugate gradient solver operators over

# Comment
- The reshape tuple is necessary in case the linearoperator! is defined over vectors of a different size as compared to what plays nicely with the dimension reduction in the ConjugateGradient. It also allows the user to define preconditioners over arrays that are more convenienently shaped.

# Return
- ConjugateGradient struct
"""
function ConjugateGradient(
    Q::AT;
    rtol = eps(eltype(Q)),
    atol = eps(eltype(Q)),
    max_iter = length(Q),
    dims = :,
    reshape_tuple = size(Q),
) where {AT}

    # allocate arrays (5)
    r0 = reshape(similar(Q), reshape_tuple)
    z0 = reshape(similar(Q), reshape_tuple)
    r1 = reshape(similar(Q), reshape_tuple)
    z1 = reshape(similar(Q), reshape_tuple)
    p1 = reshape(similar(Q), reshape_tuple)
    # allocate array of different shape (2)
    alias_p0 = similar(Q)
    alias_Lp = similar(Q)
    # allocate create aliased arrays (2)
    p0 = reshape(alias_p0, reshape_tuple)
    Lp = reshape(alias_Lp, reshape_tuple)

    container = [
        rtol,
        atol,
        r0,
        z0,
        p0,
        r1,
        z1,
        p1,
        Lp,
        alias_p0,
        alias_Lp,
        dims,
        reshape_tuple,
        max_iter,
    ]
    # create struct instance by splatting the container into the default constructor
    return ConjugateGradient{
        typeof(z0),
        typeof(Q),
        eltype(Q),
        typeof(dims),
        typeof(reshape_tuple),
        typeof(max_iter),
    }(container...)
end

# Define the outer constructor for the ConjugateGradient struct
"""
function ConjugateGradient(Q::MPIStateArray; rtol = eps(eltype(Q)), atol = eps(eltype(Q)), dims = :)

# Description
Outer constructor for the ConjugateGradient struct with MPIStateArrays. THIS IS A HACK DUE TO RESHAPE FUNCTIONALITY ON MPISTATEARRAYS.

# Arguments
- `Q`:(array). The kind of object that linearoperator! acts on.

# Keyword Arguments
- `rtol`: (float). relative tolerance
- `atol`: (float). absolute tolerance
- `dims`: (tuple or : ). the dimensions to reduce over
- `reshape_tuple`: (tuple). the dimensions that the conjugate gradient solver operators over

# Comment
- The reshape tuple is necessary in case the linearoperator! is defined over vectors of a different size as compared to what plays nicely with the dimension reduction in the ConjugateGradient. It also allows the user to define preconditioners over arrays that are more convenienently shaped.

# Return
- ConjugateGradient struct
"""
function ConjugateGradient(
    Q::MPIStateArray;
    rtol = eps(eltype(Q)),
    atol = eps(eltype(Q)),
    max_iter = length(Q),
    dims = :,
    reshape_tuple = size(Q),
)

    # create empty container for pushing struct objects
    # allocate arrays (5)
    r0 = reshape(similar(Q.data), reshape_tuple)
    z0 = reshape(similar(Q.data), reshape_tuple)
    r1 = reshape(similar(Q.data), reshape_tuple)
    z1 = reshape(similar(Q.data), reshape_tuple)
    p1 = reshape(similar(Q.data), reshape_tuple)
    # allocate array of different shape (2)
    alias_p0 = similar(Q)
    alias_Lp = similar(Q)
    # allocate create aliased arrays (2)
    p0 = reshape(alias_p0.data, reshape_tuple)
    Lp = reshape(alias_Lp.data, reshape_tuple)
    container = [
        rtol,
        atol,
        r0,
        z0,
        p0,
        r1,
        z1,
        p1,
        Lp,
        alias_p0,
        alias_Lp,
        dims,
        reshape_tuple,
        max_iter,
    ]
    # create struct instance by splatting the container into the default constructor
    return ConjugateGradient{
        typeof(z0),
        typeof(Q),
        eltype(Q),
        typeof(dims),
        typeof(reshape_tuple),
        typeof(max_iter),
    }(container...)
end

"""
LS.initialize!(linearoperator!, Q, Qrhs, solver::ColumnwisePreconditionedConjugateGradient, args...)

# Description

- This function initializes the iterative solver. It is called as part of the AbstractIterativeLinearSolver routine. SEE CODEREF for documentation on AbstractIterativeLinearSolver

# Arguments

- `linearoperator!`: (function). This applies the predefined linear operator on an array. Applies a linear operator to object "y" and overwrites object "z". The function argument i s linearoperator!(z,y, args...) and it returns nothing.
- `Q`: (array). This is an object that linearoperator! outputs
- `Qrhs`: (array). This is an object that linearoperator! acts on
- `solver`: (struct). This is a scruct for dispatch, in this case for ColumnwisePreconditionedConjugateGradient
- `args...`: (arbitrary). This is optional arguments that can be passed into linearoperator! function.

# Keyword Arguments

- There are no keyword arguments

# Return
- `converged`: (bool). A boolean to say whether or not the iterative solver has converged.
- `threshold`: (float). The value of the residual for the first timestep

# Comment
- This function does nothing for conjugate gradient

"""
function LS.initialize!(
    linearoperator!,
    Q,
    Qrhs,
    solver::ConjugateGradient,
    args...,
)

    return false, Inf
end

"""
LS.doiteration!(linearoperator!, Q, Qrhs, solver::ColumnwisePreconditionedConjugateGradient, threshold, args...; applyPC!)

# Description

- This function enacts the iterative solver. It is called as part of the AbstractIterativeLinearSolver routine. SEE CODEREF for documentation on AbstractIterativeLinearSolver

# Arguments

- `linearoperator!`: (function). This applies the predefined linear operator on an array. Applies a linear operator to object "y" and overwrites object "z". It is a function with arguments linearoperator!(z,y, args...), where "z" gets overwritten by "y" and "args..." are additional arguments passed to the linear operator. The linear operator is assumed to return nothing.
- `Q`: (array). This is an object that linearoperator! overwrites
- `Qrhs`: (array). This is an object that linearoperator! acts on. This is the rhs to the linear system
- `solver`: (struct). This is a scruct for dispatch, in this case for ConjugateGradient
- `threshold`: (float). Either an absolute or relative tolerance
- `applyPC!`: (function). Applies a preconditioner to objecy "y" and overwrites object "z". applyPC!(z,y)
- `args...`: (arbitrary). This is necessary for the linearoperator! function which has a signature linearoperator!(b, x, args....)

# Keyword Arguments

- There are no keyword arguments

# Return
- `converged`: (bool). A boolean to say whether or not the iterative solver has converged.
- `iteration`: (int). Iteration number for the iterative solver
- `threshold`: (float). The value of the residual for the first timestep

# Comment
- This function does conjugate gradient

"""
function LS.doiteration!(
    linearoperator!,
    Q,
    Qrhs,
    solver::ConjugateGradient,
    threshold,
    args...;
    applyPC! = (x, y) -> x .= y,
)

    # unroll names for convenience

    rtol = solver.rtol
    atol = solver.atol
    residual_norm = typemax(eltype(Q))
    dims = solver.dims
    converged = false

    max_iter = solver.max_iter
    r0 = solver.r0
    z0 = solver.z0
    p0 = solver.p0
    r1 = solver.r1
    z1 = solver.z1
    p1 = solver.p1
    Lp = solver.Lp
    alias_p0 = solver.alias_p0
    alias_Lp = solver.alias_Lp
    alias_Q = reshape(Q, solver.reshape_tuple)

    # Smack residual by linear operator
    linearoperator!(alias_Lp, Q, args...)
    # make sure that arrays are of the appropriate size
    alias_p0 .= Qrhs
    r0 .= p0 - Lp
    # apply the preconditioner
    applyPC!(z0, r0)
    # update p0
    p0 .= z0

    # TODO: FIX THIS
    absolute_residual = maximum(sqrt.(sum(r0 .* r0, dims = dims)))
    relative_residual =
        absolute_residual / maximum(sqrt.(sum(Qrhs .* Qrhs, dims = :)))
    # TODO: FIX THIS
    if (absolute_residual <= atol) || (relative_residual <= rtol)
        # wow! what a great guess
        converged = true
        return converged, 1, absolute_residual
    end

    for j in 1:max_iter
        linearoperator!(alias_Lp, alias_p0, args...)

        α = sum(r0 .* z0, dims = dims) ./ sum(p0 .* Lp, dims = dims)

        # Update along preconditioned direction, (note that broadcast will indeed work as expected)
        @. alias_Q += α * p0

        @. r1 = r0 - α * Lp

        # TODO: FIX THIS
        absolute_residual = maximum(sqrt.(sum(r1 .* r1, dims = dims)))
        relative_residual =
            absolute_residual / maximum(sqrt.(sum(Qrhs .* Qrhs, dims = :)))
        # TODO: FIX THIS
        converged = false
        if (absolute_residual <= atol) || (relative_residual <= rtol)
            converged = true
            return converged, j, absolute_residual
        end

        applyPC!(z1, r1)

        β = sum(z1 .* r1, dims = dims) ./ sum(z0 .* r0, dims = dims)

        # Update
        @. p0 = z1 + β * p0
        @. z0 = z1
        @. r0 = r1

    end

    # TODO: FIX THIS
    converged = true
    return converged, max_iter, absolute_residual
end

"""
LS.doiteration!(linearoperator!, Q::MPIStateArray, Qrhs::MPIStateArray, solver::ColumnwisePreconditionedConjugateGradient, threshold, args...; applyPC!)

# Description

This function enacts the iterative solver. It is called as part of the AbstractIterativeLinearSolver routine. SEE CODEREF for documentation on AbstractIterativeLinearSolver. THIS IS A HACK TO WORK WITH MPISTATEARRAYS. THE ISSUE IS WITH RESHAPE.

# Arguments

- `linearoperator!`: (function). This applies the predefined linear operator on an array. Applies a linear operator to object "y" and overwrites object "z". It is a function with arguments linearoperator!(z,y, args...), where "z" gets overwritten by "y" and "args..." are additional arguments passed to the linear operator. The linear operator is assumed to return nothing.
- `Q`: (array). This is an object that linearoperator! overwrites
- `Qrhs`: (array). This is an object that linearoperator! acts on. This is the rhs to the linear system
- `solver`: (struct). This is a scruct for dispatch, in this case for ConjugateGradient
- `threshold`: (float). Either an absolute or relative tolerance
- `applyPC!`: (function). Applies a preconditioner to objecy "y" and overwrites object "z". applyPC!(z,y)
- `args...`: (arbitrary). This is necessary for the linearoperator! function which has a signature linearoperator!(b, x, args....)

# Keyword Arguments

- There are no keyword arguments

# Return
- `converged`: (bool). A boolean to say whether or not the iterative solver has converged.
- `iteration`: (int). Iteration number for the iterative solver
- `threshold`: (float). The value of the residual for the first timestep

# Comment
- This function does conjugate gradient

"""
function LS.doiteration!(
    linearoperator!,
    Q::MPIStateArray,
    Qrhs::MPIStateArray,
    solver::ConjugateGradient,
    threshold,
    args...;
    applyPC! = (x, y) -> x .= y,
)

    # unroll names for convenience

    rtol = solver.rtol
    atol = solver.atol
    residual_norm = typemax(eltype(Q))
    dims = solver.dims
    converged = false

    max_iter = solver.max_iter
    r0 = solver.r0
    z0 = solver.z0
    p0 = solver.p0
    r1 = solver.r1
    z1 = solver.z1
    p1 = solver.p1
    Lp = solver.Lp
    alias_p0 = solver.alias_p0
    alias_Lp = solver.alias_Lp
    alias_Q = reshape(Q.data, solver.reshape_tuple)

    # Smack residual by linear operator
    linearoperator!(alias_Lp, Q, args...)
    # make sure that arrays are of the appropriate size
    alias_p0 .= Qrhs.data
    r0 .= p0 - Lp
    # apply the preconditioner
    applyPC!(z0, r0)
    # update p0
    p0 .= z0

    # TODO: FIX THIS
    absolute_residual = maximum(sqrt.(sum(r0 .* r0, dims = dims)))
    relative_residual =
        absolute_residual / maximum(sqrt.(sum(Qrhs .* Qrhs, dims = :)))
    # TODO: FIX THIS
    if (absolute_residual <= atol) || (relative_residual <= rtol)
        # wow! what a great guess
        converged = true
        return converged, 1, absolute_residual
    end

    for j in 1:max_iter
        linearoperator!(alias_Lp, alias_p0, args...)

        α = sum(r0 .* z0, dims = dims) ./ sum(p0 .* Lp, dims = dims)

        # Update along preconditioned direction, (note that broadcast will indeed work as expected)
        @. alias_Q += α * p0

        @. r1 = r0 - α * Lp

        # TODO: FIX THIS
        absolute_residual = maximum(sqrt.(sum(r1 .* r1, dims = dims)))
        relative_residual =
            absolute_residual / maximum(sqrt.(sum(Qrhs .* Qrhs, dims = :)))
        # TODO: FIX THIS
        converged = false
        if (absolute_residual <= atol) || (relative_residual <= rtol)
            converged = true
            return converged, j, absolute_residual
        end

        applyPC!(z1, r1)

        β = sum(z1 .* r1, dims = dims) ./ sum(z0 .* r0, dims = dims)

        # Update
        @. p0 = z1 + β * p0
        @. z0 = z1
        @. r0 = r1

    end

    # TODO: FIX THIS
    converged = true
    return converged, max_iter, absolute_residual
end

end #module


text:
module ComorbidityIndexes

using ICD10Utilities

export comorbiditylists

export flagcomorbidities
export scorecomorbidity

# comorbidity codes
include("quan.jl")
include("macss.jl")
include("comorbiditylists.jl")

# functions

"""
    flagcomorbidities(icdcodes, comorbidities)

Flag comorbidities in an array of ICD-10 codes.

`comorbidities` can be `:elixhauser` or `:charlson` for Quan ICD-10 versions of
Elixhauser and Charlson comorbidities, respectively, or `:macss` for the
Multipurpose Australian Comorbidity Scoring System.

Returns Dict giving true or false for each comorbidity.
"""
function flagcomorbidities(icdcodes, comorbidities)
  if !haskey(comorbiditylists, comorbidities)
    throw(DomainError(comorbidities, "Valid comorbidities are $(keys(comorbiditylists))"))
  elseif comorbidities in [:elixhauser, :charlson]
    return _quan(icdcodes, comorbiditylists[comorbidities])
  elseif comorbidities in [:macss]
    return _macss(icdcodes, comorbiditylists[comorbidities])
  end
end

function _quan(icdcodes, comorbidities)
  RET = Dict{Symbol,Bool}()
  for cmrb in eachindex(comorbidities)
    RET[cmrb] = any(icd -> icd3(icd) in icd3.(comorbidities[cmrb][1]), icdcodes) ||
                any(icd -> icd4(icd) in icd4.(comorbidities[cmrb][2]), icdcodes)
  end

  if RET[:diabetes_comp]
    RET[:diabetes] = false
  end
  if haskey(RET, :severeliverdis) && RET[:severeliverdis]
    RET[:liverdisease] = false
  end
  if RET[:metastaticcancer]
    if haskey(RET, :malignancy)
      RET[:malignancy] = false
    elseif haskey(RET, :solidtumor)
      RET[:solidtumor] = false
    end
  end
  return RET
end

function _macss(icdcodes, comorbidities)
  RET = Dict{Symbol,Bool}()
  for cmrb in eachindex(comorbidities)
    RET[cmrb] = any(icd -> icd in comorbidities[cmrb], icdcodes)
  end
  return RET
end

"""
    scorecomorbidity(comorbidities, score)

Compute summary comorbidity score.

`comorbidities` should be the output from `flagcomorbidities`(@ref).
`score` should be `:charlson` for original Charlson weights, `:charlson_quan2011` for
Quan's updated Charlson weights, or `:vanwalraven` for the van Walraven weights
for the Elixhauser comorbidities.

It is your responsibility to request a score compatible with the comorbidity
flags you provide.
"""
function scorecomorbidity(comorbidities, score)
  if score == :charlson
    return comorbidities[:myocardialinfarction] +
           comorbidities[:congestiveheartfailure] +
           comorbidities[:periphvascdis] +
           comorbidities[:cerebrovascdis] +
           comorbidities[:dementia] +
           comorbidities[:chronicpulmonarydis] +
           comorbidities[:rheumaticdis] +
           comorbidities[:pepticulcerdis] +
           comorbidities[:mildliverdis] +
           comorbidities[:diabetes] +
           2 * (comorbidities[:diabetes_comp] +
            comorbidities[:paralysis] +
            comorbidities[:renaldis] +
            comorbidities[:malignancy]) +
           3 * comorbidities[:severeliverdis] +
           6 * (comorbidities[:metastaticcancer] + comorbidities[:hivaids])
  elseif score == :charlson_quan2011
    return comorbidities[:chronicpulmonarydis] +
           comorbidities[:rheumaticdis] +
           comorbidities[:diabetes_comp] +
           comorbidities[:renaldis] +
           2 * (comorbidities[:congestiveheartfailure] +
            comorbidities[:dementia] +
            comorbidities[:mildliverdis] +
            comorbidities[:paralysis] +
            comorbidities[:malignancy]) +
           4 * (comorbidities[:severeliverdis] + comorbidities[:hivaids]) +
           6 * comorbidities[:metastaticcancer]
  elseif score == :vanwalraven
    return -7 * (comorbidities[:drugabuse]) +
           -4 * (comorbidities[:obesity]) +
           -3 * (comorbidities[:depression]) +
           -2 * (comorbidities[:bloodlossanaemia] + comorbidities[:deficiencyanaemia]) +
           -1 * (comorbidities[:valvulardisease]) +
           2 * (comorbidities[:periphvascdis]) +
           3 * (comorbidities[:chronicpulmonarydis] + comorbidities[:coagulopathy]) +
           4 * (comorbidities[:pulmonarycircdisorders] + comorbidities[:solidtumor]) +
           5 * (comorbidities[:cardiacarrythmias] +
            comorbidities[:renalfailure] +
            comorbidities[:fluidelectrolytedisorders]) +
           6 * (comorbidities[:otherneurodisorders] + comorbidities[:weightloss]) +
           7 * (comorbidities[:congestiveheartfailure] + comorbidities[:paralysis]) +
           9 * (comorbidities[:lymphoma]) +
           11 * (comorbidities[:liverdisease]) +
           12 * (comorbidities[:metastaticcancer])
  end
end

end # module


text:
"""
    bands_full(sys::TBHamiltonian{D, T}, kpath::KPath{D}) where {D, T}

Evaluate the band structure of `sys` along `kpath`, using full diagonalization.
"""
function bands_full(sys::TBHamiltonian{D, T}, kpath::KPath{D}) where {D, T}

    res = pmap(k->eigen(Hermitian(hamiltoniank_full(sys, k))), kpath.kpts)
    empty = Vector{Array{1, ComplexF64}}(undef, 0)

    return KBandAmplitudes(
                            kpath,
                            Vector{Array{1, ComplexF64}}(undef, 0),
                            [res[k].values for k=1:kpath.nkpts]
                        )

end

"""
    bands_sparse(sys::TBHamiltonian{D, T}, kpath::KPath{D}; nbands, arround, ishift=0.1, kw...) where {D, T}

Evaluate the `nbands` of the band structure of `sys` closest to the energy `arround` for the `kpath`, using ARPACK.jl function `eigs`, using the shift and invert method. `ishift` is a shift in the imaginary direction to avoid poles in the real line.
"""
function bands_sparse(sys::TBHamiltonian{D, T}, kpath::KPath{D}; nbands, arround, ishift=0.1, kw...) where {D, T}

    res = pmap(k->eigs(hamiltoniank_sparse(sys, k); nev=nbands, which=:LM, sigma = arround + ishift*im, kw...), kpath.kpts)
    empty = Vector{Array{1, ComplexF64}}(undef, 0)

    return KBandAmplitudes(
                            kpath,
                            Vector{Array{1, ComplexF64}}(undef, 0),
                            [real.(res[k][1]) for k=1:kpath.nkpts]
                        )

end
# add dos_full and dos_sparse
# add arpes_full and arpes_sparse
# add pdos_full and pdos_sparse


text:
"""
    rekons(A,B,Xawal)
adalah fungsi yang digunakan untuk mencari solusi SPL `Ax=B` secara iteratif dengan
nilai tebakan awal `Xawal` menggunakan metode rekonstruksi aljabar.

# Example
```jl
julia> A = [4 -1 1
            4 -8 1
           -2 1 5];

julia> B = [7;-21;15];

julia> Xa = [1,2,2];

julia> X,flag,err,M = rekons(A,B,Xa);

julia> X
3-element Array{Float64,1}:
 1.9999987931831584
 3.9999993707528194
 2.9999996431226994

julia> flag
0

julia> err
7.477270454885512e-8

julia> M
18×5 Array{Float64,2}:
  0.0  1.0       2.0      2.0      NaN
  1.0  0.677229  3.45151  2.58059    0.658544
  2.0  1.51207   3.75397  2.85403    0.0644652
  3.0  1.79762   3.89498  2.94005    0.0145717
  ⋮
 14.0  1.99998   3.99999  3.0        9.83924e-7
 15.0  1.99999   4.0      3.0        4.16762e-7
 16.0  2.0       4.0      3.0        1.76529e-7
 17.0  2.0       4.0      3.0        7.47727e-8
```
return solusi `X` dengan `flag` bernilai 0 jika metode rekonstruksi aljabar
berhasil menemukan solusi dan gagal jika tidak 0.
Serta, matriks `M` yang berisi catatan proses tiap iterasi `[k, X', error]`
"""
function rekons(A,B,Xawal)
    delta = 10^-7;
    maxi = 100;
    flag = 1;
    X = Xawal[:];
    Xlama = X;
    M = [0 X' NaN];
    for k = 1:maxi
        for i = 1:length(B)
            Xlama = X;
            X = Xlama + A[i,:]*(B[i]-A[i,:]'Xlama)/(A[i,:]'A[i,:])
        end
        err = norm(X-Xlama)
        M = [M; [k X' err] ];
        if err<delta
            flag = 0;
            break
        end
    end
    err = M[end,end]
    return X, flag, err, M
end


text:
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
Pkg.instantiate()
using PackageCompiler
PackageCompiler.create_sysimage([
                                 :Pluto,
                                 :MLJ,
                                 :MLJLinearModels,
#                                  :DataFrames,
#                                  :Plots,
#                                  :PlutoUI, :Distributions,
#                                  :Flux, :Zygote, :Random,
#                                  :NearestNeighborModels, :Statistics, :LinearAlgebra,
#                                  :DataFrames, :CSV, :StatsPlots
                                ];
                                 replace_default = true,
#                                 sysimage_path=joinpath(@__DIR__, "mlcourse.so"),
                                precompile_execution_file=joinpath(@__DIR__, "warmup.jl"))

max_stars_repo_name: baltisberger/MLCourse

max_issues_repo_name: baltisberger/MLCourse

max_forks_repo_name: baltisberger/MLCourse


text:
include("Density.jl");
include("BiasPotential.jl");

function calc_delta_potential_of_mean_force(umbrella_center, data, sigma, kbt, input)
    delta_pmf = -kbt * log(calc_density(input, data, sigma))
    delta_pmf += kbt * log(calc_density(umbrella_center, data, sigma))
    delta_pmf -= calc_bias_potential(umbrella_center, input) - calc_bias_potential(umbrella_center, umbrella_center)
    return delta_pmf
end

function calc_potential_of_mean_force(umbrella_centers, data, weight, sigma)
    pmf = []
    for i = 1:length(data)
        tmp = 0.0
        for j = 1:length(umbrella_centers)
              tmp += weight[j] * get_gaussian(data[i], umbrella_centers[j], sigma)
        end
        push!(pmf, tmp)
    end

    return pmf
end


text:
"""
    segmkf(n, K; rep = 1)
Build segments for K-fold cross-validation.
* `n` : Total nb. observations in the dataset. The sampling
    is implemented with 1:n.
* `K` : Nb. folds (segments) splitting the data.
* `rep` : Nb. replications of the sampling.

Build K segments splitting the data (eventually with replications) that can
be used to validate a model.

The function returns a list (vector) of rep elements.
Each element of the list contains K vectors (= K segments).
Each segment contains the indexes (position within 1:n) of the sampled observations.

## Examples
```julia
n = 10 ; m = 3 ; rep = 4
segm = segmts(n, m; rep)
i = 1
segm[i]
segm[i][1]
```
"""
function segmkf(n, K; rep = 1)
    m = K - n % K ;
    s = list(rep)
    @inbounds for i = 1:rep
        s[i] = list(K, Vector{Int64})
        v = [randperm(n) ; repeat([0], outer = m)]
        v = reshape(v, K, :)
        @inbounds for j = 1:K
            s[i][j] = sort(filter(x -> x > 0, v[j, :]))
        end
    end
    s
end

"""
    segmkf(n, K, group; rep = 1)
Build segments (with block-sampling) for K-fold cross-validation.
* `n` : Total nb. observations in the dataset. The sampling
    is implemented with 1:n.
* `group` : A vector (n,) defining the blocks.
* `K` : Nb. folds (segments) splitting the data.
* `rep` : Nb. replications of the sampling.

Build K segments splitting the data (eventually with replications) that can
be used to validate a model.

The function samples entire blocks of observations instead of observations.
Vector `group` (defining the blocks) must be of length n.

Such a block-sampling is required when data is structured by blocks and
when the response to predict is correlated within blocks.
It prevents high underestimation of the generalization error.

The function returns a list (vector) of rep elements.
Each element of the list contains K vectors (= K segments).
Each segment contains the indexes (position within 1:n) of the sampled observations.

## Examples
```julia
n = 10
group = ["A", "B", "C", "D", "E", "A", "B", "C", "D", "E"]    # The blocks of the observations
unique(group)    # Print of the blocks
K = 3 ; rep = 4
segm = segmkf(n, K, group; rep = rep)
i = 1
segm[i]
segm[i][1]
group[segm[i][1]]
group[segm[i][2]]
group[segm[i][3]]
```
"""
function segmkf(n, K, group; rep = 1)
    # n is not used but is kept for multiple dispatch
    group = vec(group) # must be of size (n,)
    s = list(rep)
    yagg = unique(group)
    zn = length(yagg)
    K = min(K, zn)
    m = K - zn % K ;
    s = list(rep)
    @inbounds for i = 1:rep
        s[i] = list(K, Vector{Int64})
        v = [randperm(zn) ; repeat([0], outer = m)]
        v = reshape(v, K, :)
        @inbounds for j = 1:K
            s[i][j] = sort(filter(x -> x > 0, v[j, :]))
        end
    end
    zs = copy(s)
    @inbounds for i = 1:rep
        @inbounds for j = 1:K
            u = s[i][j]
            v = findall(in(yagg[u]).(group))
            zs[i][j] = v
        end
    end
    zs
end

"""
    segmts(n, m; rep = 1)
Build segments for "test-set" validation.
* `n` : Total nb. observations in the dataset. The sampling
    is implemented within 1:n.
* `m` : Nb. observations in each segment.
* `rep` : Nb. replications of the sampling.

This builds a test set (eventually replicated) that can be used to validate a model.

The function returns a list (vector) of rep elements.
Each element of the list contains a vector (= segment)
of the indexes (position within 1:n) of the m sampled observations.

## Examples
```julia
n = 10 ; m = 3 ; rep = 4
segm = segmts(n, m; rep)
i = 1
segm[i]
segm[i][1]
```
"""
function segmts(n, m; rep = 1)
    s = list(rep)
    for i = 1:rep
        s[i] = list(1, Vector{Int64})
        s[i][1] = sample(1:n, m, replace = false, ordered = true)
    end
    s
end

"""
    segmts(n, m, group; rep = 1)
Build segments with block-sampling for "test-set" validation.
* `n` : Total nb. observations in the dataset. The sampling
    is implemented with 1:n.
* `group` : A vector (n,) defining the blocks.
* `m` : Nb. blocks in the segment.
* `rep` : Nb. replications of the sampling.

Build a test set (eventually replicated) that can be used to validate a model.

The function samples entire blocks of observations instead of observations.
Vector `group` (defining the blocks) must be of length n.

Such a block-sampling is required when data is structured by blocks and
when the response to predict is correlated within blocks.
It prevents high underestimation of the generalization error.

The function returns a list (Vector::Any) of rep elements.
Each element of the list contains a vector (= segment)
of the indexes (position within 1:n) of the sampled observations.
Each segment contains m blocks.

## Examples
```julia
n = 10
group = ["A", "B", "C", "D", "E", "A", "B", "C", "D", "E"]    # The blocks of the observations
unique(group)    # Print of the blocks
m = 2 ; rep = 4
segm = segmts(n, m, group; rep)
i = 1
segm[i]
segm[i][1]
group[segm[i][1]]
```
"""
function segmts(n, m, group; rep = 1)
    # n is not used but is kept for multiple dispatch
    group = vec(group) # must be of size (n,)
    s = list(rep)
    yagg = unique(group)
    zn = length(yagg)
    m = min(m, zn)
    @inbounds for i = 1:rep
        s[i] = list(1, Vector{Int64})
        s[i][1] = sample(1:zn, m, replace = false, ordered = true)
    end
    zs = copy(s)
    @inbounds for i = 1:rep
        u = s[i][1]
        v = findall(in(yagg[u]).(group))        # which(group %in% yagg[u])
        zs[i][1] = v
    end
    zs
end


text:
include("static.jl")
export Static
include("backtracking.jl")
include("hagerzhangline.jl")
export HZAW
export FFQuadInterp


text:
AP = CSV.read("../examples/AirPassengers.csv")
logAP = log.(Vector{Float64}(AP.Passengers))

@testset "Air passengers with Kalman filter" begin
    model1 = structural(logAP, 12)

    @test isa(model1, StateSpaceModel)
    @test model1.mode == "time-invariant"

    ss1 = statespace(model1)

    @test ss1.filter_type <: KalmanFilter
    @test isa(ss1, StateSpace)
    compare_forecast_simulation(ss1, 20, 1000, 1e-2)

    model2 = structural(logAP, 12)

    ss2 = statespace(model2; filter_type = SquareRootFilter{Float64})

    @test ss2.filter_type <: SquareRootFilter
    @test isa(ss2, StateSpace)
    compare_forecast_simulation(ss2, 20, 1000, 1e-2)

    model3 = structural(logAP, 12)

    ss3 = statespace(model3; filter_type = UnivariateKalmanFilter{Float64})

    @test ss3.filter_type <: UnivariateKalmanFilter
    @test isa(ss3, StateSpace)
    compare_forecast_simulation(ss3, 20, 1000, 1e-2)

    @test ss2.smoother.alpha ≈ ss1.smoother.alpha rtol = 1e-2
    @test ss3.smoother.alpha ≈ ss1.smoother.alpha rtol = 1e-3
    @test ss2.filter.a ≈ ss1.filter.a rtol = 1e-3
    @test ss3.filter.a ≈ ss1.filter.a rtol = 1e-3

    diags1 = diagnostics(ss1)
    diags2 = diagnostics(ss2)
    diags3 = diagnostics(ss3)

    @test diags1.p_jarquebera ≈ diags2.p_jarquebera atol = 1e-2
    @test diags1.p_jarquebera ≈ diags3.p_jarquebera atol = 1e-2
    @test diags1.p_ljungbox ≈ diags2.p_ljungbox atol = 1e-2
    @test diags1.p_ljungbox ≈ diags3.p_ljungbox atol = 1e-2
    @test diags1.p_homo ≈ diags2.p_homo atol = 1e-2
    @test diags1.p_homo ≈ diags3.p_homo atol = 1e-2
end


text:
export perceptron

import ..PerceptronModel

using mlpack._Internal.io

import mlpack_jll
const perceptronLibrary = mlpack_jll.libmlpack_julia_perceptron

# Call the C binding of the mlpack perceptron binding.
function perceptron_mlpackMain()
  success = ccall((:perceptron, perceptronLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module perceptron_internal
  import ..perceptronLibrary

import ...PerceptronModel

# Get the value of a model pointer parameter of type PerceptronModel.
function IOGetParamPerceptronModel(paramName::String)::PerceptronModel
  PerceptronModel(ccall((:IO_GetParamPerceptronModelPtr, perceptronLibrary), Ptr{Nothing}, (Cstring,), paramName))
end

# Set the value of a model pointer parameter of type PerceptronModel.
function IOSetParamPerceptronModel(paramName::String, model::PerceptronModel)
  ccall((:IO_SetParamPerceptronModelPtr, perceptronLibrary), Nothing, (Cstring, Ptr{Nothing}), paramName, model.ptr)
end

# Serialize a model to the given stream.
function serializePerceptronModel(stream::IO, model::PerceptronModel)
  buf_len = UInt[0]
  buf_ptr = ccall((:SerializePerceptronModelPtr, perceptronLibrary), Ptr{UInt8}, (Ptr{Nothing}, Ptr{UInt}), model.ptr, Base.pointer(buf_len))
  buf = Base.unsafe_wrap(Vector{UInt8}, buf_ptr, buf_len[1]; own=true)
  write(stream, buf)
end
# Deserialize a model from the given stream.
function deserializePerceptronModel(stream::IO)::PerceptronModel
  buffer = read(stream)
  PerceptronModel(ccall((:DeserializePerceptronModelPtr, perceptronLibrary), Ptr{Nothing}, (Ptr{UInt8}, UInt), Base.pointer(buffer), length(buffer)))
end
end # module

"""
    perceptron(; [input_model, labels, max_iterations, test, training, verbose])

This program implements a perceptron, which is a single level neural network.
The perceptron makes its predictions based on a linear predictor function
combining a set of weights with the feature vector.  The perceptron learning
rule is able to converge, given enough iterations (specified using the
`max_iterations` parameter), if the data supplied is linearly separable.  The
perceptron is parameterized by a matrix of weight vectors that denote the
numerical weights of the neural network.

This program allows loading a perceptron from a model (via the `input_model`
parameter) or training a perceptron given training data (via the `training`
parameter), or both those things at once.  In addition, this program allows
classification on a test dataset (via the `test` parameter) and the
classification results on the test set may be saved with the `predictions`
output parameter.  The perceptron model may be saved with the `output_model`
output parameter.

Note: the following parameter is deprecated and will be removed in mlpack 4.0.0:
`output`.
Use `predictions` instead of `output`.

The training data given with the `training` option may have class labels as its
last dimension (so, if the training data is in CSV format, labels should be the
last column).  Alternately, the `labels` parameter may be used to specify a
separate matrix of labels.

All these options make it easy to train a perceptron, and then re-use that
perceptron for later classification.  The invocation below trains a perceptron
on `training_data` with labels `training_labels`, and saves the model to
`perceptron_model`.

```julia
julia> using CSV
julia> training_data = CSV.read("training_data.csv")
julia> training_labels = CSV.read("training_labels.csv"; type=Int)
julia> _, perceptron_model, _ = perceptron(labels=training_labels,
            training=training_data)
```

Then, this model can be re-used for classification on the test data `test_data`.
 The example below does precisely that, saving the predicted classes to
`predictions`.

```julia
julia> using CSV
julia> test_data = CSV.read("test_data.csv")
julia> _, _, predictions = perceptron(input_model=perceptron_model,
            test=test_data)
```

Note that all of the options may be specified at once: predictions may be
calculated right after training a model, and model training can occur even if an
existing perceptron model is passed with the `input_model` parameter.  However,
note that the number of classes and the dimensionality of all data must match.
So you cannot pass a perceptron model trained on 2 classes and then re-train
with a 4-class dataset.  Similarly, attempting classification on a 3-dimensional
dataset with a perceptron that has been trained on 8 dimensions will cause an
error.

# Arguments

 - `input_model::PerceptronModel`: Input perceptron model.
 - `labels::Array{Int, 1}`: A matrix containing labels for the training
      set.
 - `max_iterations::Int`: The maximum number of iterations the perceptron
      is to be run  Default value `1000`.

 - `test::Array{Float64, 2}`: A matrix containing the test set.
 - `training::Array{Float64, 2}`: A matrix containing the training set.
 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

# Return values

 - `output::Array{Int, 1}`: The matrix in which the predicted labels for
      the test set will be written.
 - `output_model::PerceptronModel`: Output for trained perceptron model.
 - `predictions::Array{Int, 1}`: The matrix in which the predicted labels
      for the test set will be written.

"""
function perceptron(;
                    input_model::Union{PerceptronModel, Missing} = missing,
                    labels = missing,
                    max_iterations::Union{Int, Missing} = missing,
                    test = missing,
                    training = missing,
                    verbose::Union{Bool, Missing} = missing,
                    points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, perceptronLibrary), Nothing, ());

  IORestoreSettings("Perceptron")

  # Process each input argument before calling mlpackMain().
  if !ismissing(input_model)
    perceptron_internal.IOSetParamPerceptronModel("input_model", convert(PerceptronModel, input_model))
  end
  if !ismissing(labels)
    IOSetParamURow("labels", labels)
  end
  if !ismissing(max_iterations)
    IOSetParam("max_iterations", convert(Int, max_iterations))
  end
  if !ismissing(test)
    IOSetParamMat("test", test, points_are_rows)
  end
  if !ismissing(training)
    IOSetParamMat("training", training, points_are_rows)
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("output")
  IOSetPassed("output_model")
  IOSetPassed("predictions")
  # Call the program.
  perceptron_mlpackMain()

  return IOGetParamURow("output"),
         perceptron_internal.IOGetParamPerceptronModel("output_model"),
         IOGetParamURow("predictions")
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Compute the satellite eclipse time summary.
#
# References
# ==============================================================================
#
# [1] Longo, C. R. O., Rickman, S. L (1995). Method for the Calculation of
#     Spacecraft Umbra and Penumbra Shadow Terminator Points. NASA Technical
#     Paper 3547.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export eclipse_time_summary

"""
    eclipse_time_summary([io::IO,] orbp::OrbitPropagator; kwargs...)

Compute the eclipse time summary for the orbit propagator `orbp`. The summary is
computed as the total time the object stays in the sunlight, penumbra, and umbra
regions per day.

# Returns

If the `io` parameter **is not** present, then the function returns three
vectors:

- The total sunlight time per orbit at each day [s].
- The total penumbra time per orbit at each day [s].
- The total umbra time per orbit at each day [s].

If the `io` parameter **is** present, then the function prints to `io` the
results in table format.

# Keywords

- `num_days::Number`: Number of days in which the analysis will be performed.
    (**Default** = 365)
- `step::Number`: The step in which the propagation will occur. Notice that this
    function has a crossing estimation to accurately estimate the transition
    between the regions. However, if this step is very large, we may miss some
    small regions. If it is negative, then it will be selected as the time in
    which the mean anomaly advances 0.5°.

If the parameter `io` is passed, then the following additional keywords are
available:

- `unit::Symbol`: Select the unit in which the results will be printed. The
    possible values are:
    - `:s` for seconds (**Default**);
    - `:m` for minutes; or
    - `:h` for hours.
- `use_pager::Bool`: If `true`, then the result will be printed using
    **TerminalPager.jl**.
"""
function eclipse_time_summary(
    orbp::OrbitPropagator;
    num_days::Number = 365,
    step::Number = -1
)
    jd₀ = get_epoch(orbp)

    # TODO: Improve how the orbit period is computed.
    # We must obtain the mean elements to compute the orbit period. Maybe there
    # is a better way to do this.
    mean_elements = get_mean_elements(orbp)

    # We need the orbit period because we will propagate one orbit per day.
    orbit_period = period(mean_elements)

    # Check the propagation step we need to use.
    Δt₀ = step < 0 ? orbit_period * 0.5 / 360 : Float64(step)

    # Vector of the days in which the eclipse time will be computed.
    days = 0:1:num_days-1

    # Preallocate the output variables.
    sunlight_time = zeros(num_days)
    penumbra_time = zeros(num_days)
    umbra_time    = zeros(num_days)

    # Loop
    # ==========================================================================

    @inbounds for d in days
        # TODO: Should we transform between MOD => TOD/TEME?
        # Get the Sun position represented in the inertial reference frame.
        s_i = sun_position_i(jd₀ + d)

        # Initial state.
        old_state = _get_lighting_condition(orbp, 0, d, s_i)

        # Compute the eclipse time during one orbit.
        Δt  = Δt₀
        t_k = Δt

        while true
            new_state = _get_lighting_condition(orbp, t_k, d, s_i)

            # Check if the state has changed.
            if new_state != old_state
                # Refine to find the edge.
                t_k₀ = t_k - Δt
                t_k₁ = t_k
                t_kc = find_crossing(
                    _fc,
                    t_k₀,
                    t_k₁,
                    true,
                    false,
                    orbp,
                    d,
                    s_i,
                    old_state
                )

                # Times to be added in the old and new states.
                Δts₀ = t_kc - t_k₀
                Δts₁ = t_k₁ - t_kc

                _accumulate(
                    Δts₀,
                    old_state,
                    d + 1,
                    sunlight_time,
                    penumbra_time,
                    umbra_time
                )

                _accumulate(
                    Δts₁,
                    new_state,
                    d + 1,
                    sunlight_time,
                    penumbra_time,
                    umbra_time
                )

            # If not, just add the time step to the current state.
            else
                _accumulate(
                    Δt,
                    new_state,
                    d + 1,
                    sunlight_time,
                    penumbra_time,
                    umbra_time
                )
            end

            old_state = new_state

            abs(t_k - orbit_period) < 1e-3 && break

            # Make sure that the last interval will have the exact size so that
            # the end of the analysis is the end of the orbit.
            if (t_k + Δt > orbit_period)
                Δt = orbit_period - t_k
            end

            t_k += Δt
        end
    end

    return sunlight_time, penumbra_time, umbra_time
end

function eclipse_time_summary(
    io::IO,
    orbp::OrbitPropagator;
    num_days::Number = 365,
    step::Number = -1,
    unit::Symbol = :s,
    use_pager::Bool = false
)
    # Compute the summary.
    sunlight_time, penumbra_time, umbra_time = eclipse_time_summary(
        orbp;
        num_days,
        step
    )

    # Get the date vector.
    vjd = get_epoch(orbp) .+ collect(0:1:num_days-1)
    vdate = string.(jd_to_date.(DateTime, vjd))

    # Obtain the conversion given the unit selection.
    if unit == :h
        unit_label = "h"
        unit_factor = 3600
    elseif unit == :m
        unit_label = "min"
        unit_factor = 60
    else
        unit_label = "s"
        unit_factor = 1
    end

    sunlight_time ./= unit_factor
    penumbra_time ./= unit_factor
    umbra_time ./= unit_factor

    str = pretty_table(
        String,
        vcat(
            hcat(vdate, sunlight_time, penumbra_time, umbra_time),
            hcat(
                AnsiTextCell(string(crayon"red bold") * "Maximum"),
                maximum(sunlight_time),
                maximum(penumbra_time),
                maximum(umbra_time)
            ),
            hcat(
                AnsiTextCell(string(crayon"yellow bold") * "Mean"),
                mean(sunlight_time),
                mean(penumbra_time),
                mean(umbra_time)
            ),
            hcat(
                AnsiTextCell(string(crayon"blue bold") * "Minimum"),
                minimum(sunlight_time),
                minimum(penumbra_time),
                minimum(umbra_time)
            ),
        );
        color = get(stdout, :color, false),
        crop = use_pager ? :none : :horizontal,
        display_size = use_pager ? (-1, -1) : displaysize(stdout),
        header = [
            "Day",
            "Sunlight time [$unit_label]",
            "Penumbra time [$unit_label]",
            "Umbra time [$unit_label]"
        ],
        hlines = [:header, num_days + 1],
        alignment_anchor_regex = Dict(
            2 => [r"\."],
            3 => [r"\."],
            4 => [r"\."]
        ),
        vlines = [1]
    )

    # Check if we need to use the pager.
    if use_pager
        pager(
            str;
            freeze_columns = max(length.(vdate)..., 6) + 3,
            freeze_rows = 2
        )
    else
        println(str)
    end

    return nothing
end

#                              Private functions
# ==============================================================================

function _accumulate(
    Δts::Number,
    state::Symbol,
    ind::Integer,
    sunlight_time::AbstractArray,
    penumbra_time::AbstractArray,
    umbra_time::AbstractArray
)
    @inbounds if state == :sunlight
        sunlight_time[ind] += Δts
    elseif state == :penumbra
        penumbra_time[ind] += Δts
    elseif state == :umbra
        umbra_time[ind] += Δts
    end

    return nothing
end

# Without the `@noinline` we get **a lot** of allocations when calling the
# function `find_crossing`.
@noinline function _get_lighting_condition(
    orbp::OrbitPropagator,
    t::Number,
    d::Number,
    s_i::AbstractVector
)
    r_i, ~ = propagate!(orbp, 86400d + t)
    return lighting_condition(r_i, s_i)
end

@noinline function _fc(
    t::Number,
    orbp::OrbitPropagator,
    d::Number,
    s_i::AbstractVector,
    old_state::Symbol
)
    return _get_lighting_condition(orbp, t, d, s_i) == old_state
end


text:
using DifferentialEquations, Flux, DiffEqFlux
using Plots, NPZ
using IterTools: ncycle
using BSON: @save, @load

true_data = npzread("data.npy")[1:20:4000,:]
data_len = size(true_data)[1]
x_data = true_data[:,1]
state_data = true_data[:,2:3]

u0 = [true_data[1,2],true_data[1,3]]
tspan = (0.0,true_data[data_len,1])
tsteps = 0.0:true_data[data_len,1]-true_data[data_len-1,1]:true_data[data_len,1]

my_nn = Chain(
  Dense(2, 32, relu),
  Dense(32, 2)) |> f32

p, re = Flux.destructure(my_nn)

function right_hand_side(du,u,p,t)
  m = re(p)
  nn_output = m(u)
  du[1] = nn_output[1]
  du[2] = nn_output[2]
end

prob = ODEProblem(right_hand_side,u0,tspan,p)
sol = solve(prob, Tsit5())

# Plot the solution
plot(sol)
savefig("Initial_NN_ODE.png")

function predict_ode(state_batch,time_batch)
  tmp_prob = remake(prob; p = p,u=state_batch[1])
  pred = Array(solve(tmp_prob, Tsit5(), saveat = time_batch))
  return pred
end

function loss_ode(state_batch,time_batch)
    pred = predict_ode(state_batch,time_batch)
    loss = sum(abs2, state_batch-pred)
    return loss
end

k = 20
train_loader = Flux.Data.DataLoader((transpose(state_data), tsteps), batchsize = k)
numEpochs = 200
losses=[]

# test_batch = iterate(train_loader)
# print(check[1][2])
# print(size(test_batch[1][1]))
# loss_ode(test_batch[1][1],test_batch[1][2])

cb() = begin

  l=loss_ode(transpose(state_data), tsteps)
  push!(losses, l)
  @show l

  pred = predict_ode(transpose(state_data),tsteps)

  plt = plot(x_data,hcat(transpose(pred),state_data),
    line=(4, [:solid :solid :dash :dash]), ylim = (-6, 6),
    label=["y1-pred" "y2-pred" "y1-true" "y2-true"])

  display(plt)
  # Tell sciml_train to not halt the optimization. If return true, then
  # optimization stops.
  return false
end

opt=ADAM(0.01)
Flux.train!(loss_ode, Flux.params(p), ncycle(train_loader,numEpochs), opt, cb=Flux.throttle(cb, 10))

my_nn_final = re(p)
@save "mymodel.bson" my_nn_final

# Inference
@load "mymodel.bson" my_nn_final
p, re = Flux.destructure(my_nn_final)

prob_final = ODEProblem(right_hand_side, u0, tspan, p)
sol_final = solve(prob_final,Tsit5())
plot(sol_final)
savefig("NN_ode_optimized.png") # Both state components close to 1

max_stars_repo_name: Romit-Maulik/Tutorials-Demos-Practice

max_issues_repo_name: Romit-Maulik/Tutorials-Demos-Practice

max_forks_repo_name: Romit-Maulik/Tutorials-Demos-Practice


text:
"""
    isinERA5Region(
        point :: Point2{<:Real},
        e5geo :: ERA5Region;
        tlon  :: Real = 0,
        tlat  :: Real = 0,
        throw :: Bool = true
    ) -> Bool

Check if a geographical point `Point` is within an ERA5Region defined by `e5geo`.

Arguments
=========

- `point` : A geographical point of Type `Point2`.  Pass `Point2(plon,plat)`, where `plon` and `plat` are the longitude and latitudes of the point.
- `e5geo` : The ERA5Region struct container

Keyword Arguments
=================

- `tlon`  : Threshold for the longitude bound
- `tlat`  : Threshold for the latitude bound
- `throw` : If `true`, then if `Point` is not within `geo`, an error is thrown and the program stops running.
"""
function isinERA5Region(
    point :: Point2{<:Real},
    e5geo :: ERA5Region;
    tlon  :: Real = 0,
    tlat  :: Real = 0,
    throw :: Bool = true
)

    return isinGeoRegion(point,e5geo.geo,tlon=tlon,tlat=tlat,throw=throw)

end

"""
    isinERA5Region(
        geo    :: GeoRegion,
        e5geo  :: ERA5Region;
        domask :: Bool = false,
        throw  :: Bool = true
    ) -> Bool

Check if a child GeoRegion defined by `geo` is within a ERA5Region `e5geo`.

Arguments
=========

- `geo`   : A GeoRegion that we postulate to be a "child", or a subset of the ERA5Region defined by `e5geo`
- `e5geo` : An ERA5Region that we postulate to be a "parent", or containing the GeoRegion defined by `geo`

Keyword Arguments
=================

- `throw`  : If `true`, then if `geo` is not within `e5geo`, an error is thrown and the program stops running
- `domask` : If `throw` is `false` and `domask` is `true`, return a mask (with bounds defined by the `geo` GeoRegion) showing the region where `geo` and `e5geo` do not overlap
"""
function isinERA5Region(
    geo    :: GeoRegion,
    e5geo  :: ERA5Region;
    domask :: Bool = false,
    throw  :: Bool = true
)

    return isinGeoRegion(geo,e5geo.geo,domask=domask,throw=throw)

end

"""
    ERA5RegionGrid(
        e5geo :: ERA5Region{ST,FT},
        lon   :: Vector{<:Real},
        lat   :: Vector{<:Real}
    ) -> GeoRegion.RegionGrid

Creates a RegionGrid containing information and mask information required to extract regional data for the ERA5Region from the raw data.

Arguments
=========

- `e5geo` : A ERA5Region struct type
- `lon`   : A vector containing the longitude points
- `lat`   : A vector containing the latitude points
"""
function ERA5RegionGrid(
    e5geo :: ERA5Region,
    lon   :: Vector{<:Real},
    lat   :: Vector{<:Real}
)

    return RegionGrid(e5geo.geo,lon,lat)

end


text:
module SuiteSparseGraphBLAS
__precompile__(true)
using Libdl: dlsym, dlopen, dlclose
using Preferences
include("find_binary.jl")
const libgraphblas_handle = Ref{Ptr{Nothing}}()
@static if artifact_or_path == "default"
    using SSGraphBLAS_jll
    const libgraphblas = SSGraphBLAS_jll.libgraphblas
else
    const libgraphblas = artifact_or_path
end
using SparseArrays
using SparseArrays: nonzeroinds
using MacroTools
using LinearAlgebra
using Random: randsubseq, default_rng, AbstractRNG, GLOBAL_RNG
using CEnum
using SpecialFunctions: lgamma, gamma, erf, erfc
using Base.Broadcast
include("abstracts.jl")
include("libutils.jl")
include("lib/LibGraphBLAS.jl")
using .libgb

include("operators/libgbops.jl")
include("types.jl")
include("gbtypes.jl")

include("operators/operatorutils.jl")
include("operators/unaryops.jl")
include("operators/binaryops.jl")
include("operators/monoids.jl")
include("operators/semirings.jl")
include("operators/selectops.jl")
include("descriptors.jl")
using .UnaryOps
using .BinaryOps
using .Monoids
using .Semirings
_createunaryops()
_createbinaryops()
_createmonoids()
_createsemirings()
include("operators/oplist.jl")
include("indexutils.jl")

export lgamma, gamma, erf, erfc #reexport of SpecialFunctions.
export frexpe, frexpx, positioni, positionj #UnaryOps not found in Julia/stdlibs.
#BinaryOps not found in Julia/stdlibs.
export second, rminus, pair, iseq, isne, isgt, islt, isge, isle, ∨, ∧, lxor, fmod, firsti,
    firstj, secondi, secondj
#SelectOps not found in Julia/stdlibs
export offdiag

const GBVecOrMat{T} = Union{GBVector{T}, GBMatrix{T}}
const GBMatOrTranspose{T} = Union{GBMatrix{T}, Transpose{<:Any, GBMatrix{T}}}
const GBArray{T} = Union{GBVector{T}, GBMatOrTranspose{T}}
const ptrtogbtype = Dict{Ptr, AbstractGBType}()

const GrBOp = Union{
    libgb.GrB_Monoid,
    libgb.GrB_UnaryOp,
    libgb.GrB_Semiring,
    libgb.GrB_BinaryOp,
    libgb.GxB_SelectOp
}

const TypedOp = Union{
    TypedUnaryOperator,
    TypedBinaryOperator,
    TypedMonoid,
    TypedSemiring
}

const MonoidBinaryOrRig = Union{
    TypedMonoid,
    TypedSemiring,
    TypedBinaryOperator,
    AbstractSemiring,
    AbstractBinaryOp,
    AbstractMonoid
}

const OperatorUnion = Union{
    AbstractOp,
    GrBOp
}

export T1, T0, T0T1, C, CT1, CT0, CT0T1, S, ST1, ST0, ST0T1, SC, SCT1, SCT0, SCT0T1, R, RT1,
    RT0, RT0T1, RC, RCT1, RCT0, RCT0T1, RS, RST1, RST0, RST0T1, RSC, RSCT1, RSCT0, RSCT0T1
include("scalar.jl")
include("vector.jl")
include("matrix.jl")
include("random.jl")

include("operations/operationutils.jl")
include("operations/transpose.jl")
include("operations/mul.jl")
include("operations/ewise.jl")
include("operations/map.jl")
include("operations/select.jl")
include("operations/reduce.jl")
include("operations/kronecker.jl")
include("operations/concat.jl")
include("operations/resize.jl")

include("print.jl")
include("import.jl")
include("export.jl")
include("options.jl")
#EXPERIMENTAL
include("operations/argminmax.jl")
include("operations/broadcasts.jl")
include("chainrules/chainruleutils.jl")
include("chainrules/mulrules.jl")
include("chainrules/ewiserules.jl")
include("chainrules/maprules.jl")
include("chainrules/reducerules.jl")
include("chainrules/selectrules.jl")
include("chainrules/constructorrules.jl")
#include("random.jl")
include("misc.jl")
export libgb
export UnaryOps, BinaryOps, Monoids, Semirings #Submodules
export UnaryOp, BinaryOp, Monoid, Semiring #UDFs
export Descriptor #Types
export xtype, ytype, ztype, validtypes #Determine input/output types of operators
export GBScalar, GBVector, GBMatrix #arrays
export clear!, extract, extract!, subassign!, assign!, hvcat! #array functions

#operations
export mul, select, select!, eadd, eadd!, emul, emul!, map, map!, gbtranspose, gbtranspose!,
gbrand
# Reexports.
export diag, diagm, mul!, kron, kron!, transpose, reduce, tril, triu
export nnz, sprand, findnz, nonzeros, nonzeroinds
function __init__()
    @static if artifact_or_path != "default"
        libgraphblas_handle[] = dlopen(libgraphblas)
    else
        #The artifact does dlopen for us.
        libgraphblas_handle[] = SSGraphBLAS_jll.libgraphblas_handle
    end
    _load_globaltypes()
    # We initialize GraphBLAS by giving it Julia's GC wrapped memory management functions.
    # In the future this should hopefully allow us to do no-copy passing of arrays between Julia and SS:GrB.
    # In the meantime it helps Julia respond to memory pressure from SS:GrB and finalize things in a timely fashion.
    libgb.GxB_init(libgb.GrB_NONBLOCKING, cglobal(:jl_malloc), cglobal(:jl_calloc), cglobal(:jl_realloc), cglobal(:jl_free), true)
    _loaddescriptors()
    _loadselectops()
    # Set printing to base-1 rather than base-0.
    gbset(BASE1, true)
    atexit() do
        # Finalize the lib. Frees a small internal memory pool.
        libgb.GrB_finalize()
        @static if artifact_or_path != "default"
            dlclose(libgraphblas_handle[])
        end
    end
end

include("operators/ztypes.jl")
end #end of module


text:
include(joinpath(@__DIR__, "..", "src/dynamics", "flamingo", "visuals.jl"))
T = Float64
vis = Visualizer()
open(vis)

s = get_simulation("flamingo", "flat_2D_lc", "flat")
model = s.model
env = s.env
const ContactImplicitMPC = Main
ref_traj_ = deepcopy(ContactImplicitMPC.get_trajectory(s.model, s.env,
    joinpath(module_dir(), "src/dynamics/flamingo/gaits/gait_forward_36_4.jld2"),
    load_type = :split_traj_alt))
ref_traj = deepcopy(ref_traj_)

H = ref_traj.H
h = ref_traj.h
N_sample = 5
H_mpc = 15
h_sim = h / N_sample
H_sim = 1000 #35000

p = pratt_policy(model, h_sim)

# q1_ref = copy(ref_traj.q[2])
# q0_ref = copy(ref_traj.q[1])
q1_ref = [0.0, 0.849, -0.00, 0.1, 0.295, -0.3, 0.1, π/2, π/2]
q0_ref = [0.0, 0.849, -0.00, 0.1, 0.295, -0.3, 0.1, π/2, π/2]
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

sim = simulator(s, q0_sim, q1_sim, h_sim, H_sim,
    p = p,
	ip_opts = InteriorPointOptions(
		undercut = Inf,
		γ_reg = 0.0,
		# verbose = true,
        r_tol = 1.0e-8,
        κ_tol = 1.0e-8),
    sim_opts = SimulatorOptions(warmstart = true),
    )

@time status = simulate!(sim, verbose = true)
l = 9
lu = 1
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.([q[l:l] for q in sim.traj.q])...)', color=:blue, linewidth=1.0)
plot!(plt[2,1], hcat(Vector.(vcat([fill(ref_traj.u[i][lu:lu], N_sample) for i=1:H]...))...)',
    color=:red, linewidth=3.0)
plot!(plt[3,1], hcat(Vector.(vcat([fill(ref_traj.γ[i][1:nc], N_sample) for i=1:H]...))...)',
    color=:red, linewidth=3.0)
plot!(plt[2,1], hcat(Vector.([u[lu:lu] 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!(plt[3,1], hcat(Vector.([b[1:nb] for b in sim.traj.b]*N_sample)...)', color=:red, linewidth=1.0)

plot_surface!(vis, env, xlims=[-0.5, 1.5], ylims = [-0.5, 0.5])
plot_lines!(vis, model, sim.traj.q)
# anim = visualize_robot!(vis, model, sim.traj, sample=10)
anim = visualize_meshrobot!(vis, model, sim.traj, sample=2)
anim = visualize_force!(vis, model, env, sim.traj, anim=anim, h=h_sim, sample=10)

# filename = "flamingo_100_steps"
# MeshCat.convert_frames_to_video(
#     "/home/simon/Downloads/$filename.tar",
#     "/home/simon/Documents/$filename.mp4", overwrite=true)
#
# convert_video_to_gif(
#     "/home/simon/Documents/$filename.mp4",
#     "/home/simon/Documents/$filename.gif", overwrite=true)
#
#
# settransform!(vis["/Cameras/default"],
# 		compose(Translation(0.0, 0.5, -1.0),LinearMap(RotZ(-pi / 2.0))))
#
# # Ghost
# flamingo_ghost!(vis, sim, x -> 0.0)
#
# # Animation
# flamingo_animation!(vis, sim, x -> 0.0)

qref = [0.0, 0.849, -0.00, 0.1, 0.295, -0.3, 0.1, π/2, π/2]
set_robot!(vis, model, qref)
qref = [0.0, 0.849, -0.00, -0.3, 0.1, 0.1, 0.295, π/2, π/2]
set_robot!(vis, model, qref)


text:
function evolve_boundary!(
    KS::SolverSet,
    ctr::AM{TC},
    a1face::AM{TF},
    a2face::AM{TF},
    dt,
    mode,
    bc,
) where {TC<:Union{ControlVolume2F,ControlVolume2D2F},TF<:Union{Interface2F,Interface2D2F}}

    nx, ny, dx, dy = begin
        if KS.ps isa CSpace2D
            KS.ps.nr, KS.ps.nθ, KS.ps.dr, KS.ps.darc
        else
            KS.ps.nx, KS.ps.ny, KS.ps.dx, KS.ps.dy
        end
    end

    bcs = ifelse(bc isa Symbol, [bc, bc, bc, bc], bc)

    if bcs[1] == :maxwell
        @inbounds Threads.@threads for j = 1:ny
            n = -KS.ps.n[1, j, 4]
            vn, vt = local_velocity(KS.vs.u, KS.vs.v, n)
            xc = (KS.ps.vertices[1, j, 1, 1] + KS.ps.vertices[1, j, 4, 1]) / 2
            yc = (KS.ps.vertices[1, j, 1, 2] + KS.ps.vertices[1, j, 4, 2]) / 2
            bcL = local_frame(KS.ib.bc(xc, yc, KS.ib.p), n)
            flux_boundary_maxwell!(
                a1face[1, j].fw,
                a1face[1, j].fh,
                a1face[1, j].fb,
                bcL, # left
                ctr[1, j].h,
                ctr[1, j].b,
                vn,
                vt,
                KS.vSpace.weights,
                KS.gas.K,
                dt,
                dy[1, j],
                1,
            )
            a1face[1, j].fw .= global_frame(a1face[1, j].fw, n)
        end
    end

    if bcs[2] == :maxwell
        @inbounds Threads.@threads for j = 1:ny
            n = KS.ps.n[nx, j, 2]
            vn, vt = local_velocity(KS.vs.u, KS.vs.v, n)
            xc = (KS.ps.vertices[nx, j, 2, 1] + KS.ps.vertices[nx, j, 3, 1]) / 2
            yc = (KS.ps.vertices[nx, j, 2, 2] + KS.ps.vertices[nx, j, 3, 2]) / 2
            bcR = local_frame(KS.ib.bc(xc, yc, KS.ib.p), n)
            flux_boundary_maxwell!(
                a1face[nx+1, j].fw,
                a1face[nx+1, j].fh,
                a1face[nx+1, j].fb,
                bcR, # right
                ctr[nx, j].h,
                ctr[nx, j].b,
                vn,
                vt,
                KS.vSpace.weights,
                KS.gas.K,
                dt,
                dy[nx, j],
                -1,
            )
            a1face[nx+1, j].fw .= global_frame(a1face[nx+1, j].fw, n)
        end
    end

    if bcs[3] == :maxwell
        @inbounds Threads.@threads for i = 1:nx
            n = -KS.ps.n[i, 1, 1]
            vn, vt = local_velocity(KS.vs.u, KS.vs.v, n)
            xc = (KS.ps.vertices[i, 1, 1, 1] + KS.ps.vertices[i, 1, 2, 1]) / 2
            yc = (KS.ps.vertices[i, 1, 1, 2] + KS.ps.vertices[i, 1, 2, 2]) / 2
            bcD = local_frame(KS.ib.bc(xc, yc, KS.ib.p), n)

            flux_boundary_maxwell!(
                a2face[i, 1].fw,
                a2face[i, 1].fh,
                a2face[i, 1].fb,
                bcD, # left
                ctr[i, 1].h,
                ctr[i, 1].b,
                vn,
                vt,
                KS.vSpace.weights,
                KS.gas.K,
                dt,
                dx[i, 1],
                1,
            )
            a2face[i, 1].fw .= global_frame(a2face[i, 1].fw, n)
        end
    end

    if bcs[4] == :maxwell
        @inbounds Threads.@threads for i = 1:nx
            n = KS.ps.n[i, ny, 3]
            vn, vt = local_velocity(KS.vs.u, KS.vs.v, n)
            xc = (KS.ps.vertices[i, ny, 3, 1] + KS.ps.vertices[i, ny, 4, 1]) / 2
            yc = (KS.ps.vertices[i, ny, 3, 2] + KS.ps.vertices[i, ny, 4, 2]) / 2
            bcU = local_frame(KS.ib.bc(xc, yc, KS.ib.p), n)

            flux_boundary_maxwell!(
                a2face[i, ny+1].fw,
                a2face[i, ny+1].fh,
                a2face[i, ny+1].fb,
                bcU, # right
                ctr[i, ny].h,
                ctr[i, ny].b,
                vn,
                vt,
                KS.vSpace.weights,
                KS.gas.K,
                dt,
                dx[i, ny],
                -1,
            )
            a2face[i, ny+1].fw .= global_frame(a2face[i, ny+1].fw, n)
        end
    end

end


text:
@testset "Constructs_test" begin

    mt_a = Mod("a")
    mt_b = Mod("b")
    mt_c = Mod("c")
    mt_d = Mod("d")
    mt_e = Mod("e")
    # some test constructs
    oc1 = OrderedConstruct([mt_a,mt_d,mt_b])
    oc2 = OrderedConstruct([mt_a,mt_b,mt_e,mt_d])
    unoc1 = UnorderedConstruct([mt_e,mt_d,mt_c])
    unoc2 = UnorderedConstruct([mt_e,mt_a,mt_d,mt_e])
     @testset "Base" begin
        @test length(oc1) == 3
        @test length(oc2) == 4
        @test length(unoc1) == 3
        @test length(unoc2) == 4

        @test oc1[2] == mt_d
        @test oc2[3] == mt_e

        @test isequal(oc1[1:2],OrderedConstruct([mt_a,mt_d])) == true
        @test isequal(oc1[1:2],OrderedConstruct([mt_a,mt_d])) == true
        @test length(collect(oc1)) == length(oc1)
        @test length(collect(unoc1)) == length(unoc1)
        @test eltype(collect(oc1)) == typeof(mt_a)
    end

    @testset "*_ordered" begin
        @test isequal(mt_a*mt_b,OrderedConstruct([mt_a,mt_b])) == true
        @test isequal(OrderedConstruct([mt_a,mt_b])*mt_c,OrderedConstruct([mt_a,mt_b,mt_c])) == true
        @test isequal(mt_c*OrderedConstruct([mt_a,mt_b]),OrderedConstruct([mt_c,mt_a,mt_b])) == true
        @test isequal((OrderedConstruct([mt_a,mt_b])*OrderedConstruct([mt_d,mt_e])),OrderedConstruct([mt_a,mt_b,mt_d,mt_e])) == true
    end

    @testset "+_unordered" begin
        @test isequal((mt_a+mt_b),UnorderedConstruct([mt_a,mt_b])) == true
        @test isequal(UnorderedConstruct([mt_a,mt_b])+mt_c,UnorderedConstruct([mt_a,mt_b,mt_c])) == true
        @test isequal(mt_c+UnorderedConstruct([mt_a,mt_b]),UnorderedConstruct([mt_c,mt_a,mt_b])) == true
        @test isequal((UnorderedConstruct([mt_a,mt_b])+UnorderedConstruct([mt_d,mt_e])),UnorderedConstruct([mt_a,mt_b,mt_d,mt_e])) == true
    end

end

max_stars_repo_name: SamDier/BOMoD

max_issues_repo_name: SamDier/BOMoD

max_forks_repo_name: SamDier/BOMoD


text:
module LeagueV4 # LOLTools

using ..LOLTools: AbstractDTO, lol_api_server, http_action, call_api, nothing_in_event
using HTTP

struct LeagueListDTO <: AbstractDTO
    leagueId
    tier
    entries::Vector{NamedTuple{(:summonerId, :summonerName, :leaguePoints, :rank, :wins, :losses, :veteran, :inactive, :freshBlood, :hotStreak)}}
    queue
    name
    LeagueListDTO(leagueId, tier, entries, queue, name) = new(leagueId, tier, entries, queue, name)
end

struct LeagueEntryDTO <: AbstractDTO
    queueType
    summonerName
    hotStreak
    miniSeries
    wins
    veteran
    losses
    rank
    leagueId
    inactive
    freshBlood
    tier
    summonerId
    leaguePoints
    LeagueEntryDTO(queueType, summonerName, hotStreak, miniSeries, wins, veteran, losses, rank, leagueId, inactive, freshBlood, tier, summonerId, leaguePoints) = new(queueType, summonerName, hotStreak, miniSeries, wins, veteran, losses, rank, leagueId, inactive, freshBlood, tier, summonerId, leaguePoints)
end

"""
         challengerleagues(api_key::String,
                           platform::String,
                           queue::String ;
                           endpoint::HTTP.URI = lol_api_server(platform),
                           action::Function = http_action,
                           event::Function = nothing_in_event)::LeagueListDTO
"""
function challengerleagues(api_key::String,
                           platform::String,
                           queue::String ;
                           endpoint::HTTP.URI = lol_api_server(platform),
                           action::Function = http_action,
                           event::Function = nothing_in_event)::LeagueListDTO
    call_api(LeagueListDTO, api_key, action, endpoint, "/lol/league/v4/challengerleagues/by-queue/$queue", event, challengerleagues)
end

"""
        grandmasterleagues(api_key::String,
                           platform::String,
                           queue::String ;
                           endpoint::HTTP.URI = lol_api_server(platform),
                           action::Function = http_action,
                           event::Function = nothing_in_event)::LeagueListDTO
"""
function grandmasterleagues(api_key::String,
                           platform::String,
                           queue::String ;
                           endpoint::HTTP.URI = lol_api_server(platform),
                           action::Function = http_action,
                           event::Function = nothing_in_event)::LeagueListDTO
    call_api(LeagueListDTO, api_key, action, endpoint, "/lol/league/v4/grandmasterleagues/by-queue/$queue", event, grandmasterleagues)
end

"""
         masterleagues(api_key::String,
                       platform::String,
                       queue::String ;
                       endpoint::HTTP.URI = lol_api_server(platform),
                       action::Function = http_action,
                       event::Function = nothing_in_event)::LeagueListDTO
"""
function masterleagues(api_key::String,
                       platform::String,
                       queue::String ;
                       endpoint::HTTP.URI = lol_api_server(platform),
                       action::Function = http_action,
                       event::Function = nothing_in_event)::LeagueListDTO
    call_api(LeagueListDTO, api_key, action, endpoint, "/lol/league/v4/masterleagues/by-queue/$queue", event, masterleagues)
end

"""
         entries_by_summoner_id(api_key::String,
                                platform::String,
                                encryptedSummonerId::String ;
                                endpoint::HTTP.URI = lol_api_server(platform),
                                action::Function = http_action,
                                event::Function = nothing_in_event)::Set{LeagueEntryDTO}
"""
function entries_by_summoner_id(api_key::String,
                                platform::String,
                                encryptedSummonerId::String ;
                                endpoint::HTTP.URI = lol_api_server(platform),
                                action::Function = http_action,
                                event::Function = nothing_in_event)::Set{LeagueEntryDTO}
    call_api(Set{LeagueEntryDTO}, api_key, action, endpoint, "/lol/league/v4/entries/by-summoner/$encryptedSummonerId", event, entries_by_summoner_id)
end

"""
         entries_by_queue_tier_division(api_key::String,
                                        platform::String,
                                        queue::String,
                                        tier::String,
                                        division::String ;
                                        page::Int = 1,
                                        endpoint::HTTP.URI = lol_api_server(platform),
                                        action::Function = http_action,
                                        event::Function = nothing_in_event)::Set{LeagueEntryDTO}
"""
function entries_by_queue_tier_division(api_key::String,
                                        platform::String,
                                        queue::String,
                                        tier::String,
                                        division::String ;
                                        page::Int = 1,
                                        endpoint::HTTP.URI = lol_api_server(platform),
                                        action::Function = http_action,
                                        event::Function = nothing_in_event)::Set{LeagueEntryDTO}
    path = string(merge(HTTP.URI(path="/lol/league/v4/entries/$queue/$tier/$division"), query=HTTP.escapeuri("page", page)))
    call_api(Set{LeagueEntryDTO}, api_key, action, endpoint, path, event, entries_by_queue_tier_division)
end

end # module LOLTools.LeagueV4