text:
import Base: rand

export Zonotope,
       scale,
       scale!,
       reduce_order,
       remove_zero_generators,
       quadratic_map

using LazySets.Arrays: _vector_type, _matrix_type

"""
    Zonotope{N<:Real, VN<:AbstractVector{N}, MN<:AbstractMatrix{N}} <: AbstractZonotope{N}

Type that represents a zonotope.

### Fields

- `center`     -- center of the zonotope
- `generators` -- matrix; each column is a generator of the zonotope

### Notes

Mathematically, a zonotope is defined as the set

```math
Z = \\left\\{ x ∈ \\mathbb{R}^n : x = c + ∑_{i=1}^p ξ_i g_i,~~ ξ_i \\in [-1, 1]~~ ∀ i = 1,…, p \\right\\},
```
where ``c \\in \\mathbb{R}^n`` is its *center* and ``\\{g_i\\}_{i=1}^p``,
``g_i \\in \\mathbb{R}^n``, is the set of *generators*.
This characterization defines a zonotope as the finite Minkowski sum of line
segments.
Zonotopes can be equivalently described as the image of a unit infinity-norm
ball in ``\\mathbb{R}^n`` by an affine transformation.

Zonotopes can be constructed in two different ways: either passing the generators
as a matrix, where each column represents a generator, or passing a list of vectors
where each vector represents a generator. Below we illustrate both ways.

### Examples

A two-dimensional zonotope with given center and set of generators:

```jldoctest zonotope_label
julia> Z = Zonotope([1.0, 0.0], [0.1 0.0; 0.0 0.1])
Zonotope{Float64,Array{Float64,1},Array{Float64,2}}([1.0, 0.0], [0.1 0.0; 0.0 0.1])

julia> dim(Z)
2

julia> center(Z)
2-element Array{Float64,1}:
 1.0
 0.0

julia> genmat(Z)
2×2 Array{Float64,2}:
 0.1  0.0
 0.0  0.1
```
Here, the first vector in the `Zonotope` constructor corresponds to the zonotope's
center, and each column of the second argument corresponds to a generator. The
functions `center` and `genmat` return the center and the generator matrix of this
zonotope respectively.

We can collect its vertices using `vertices_list`:

```jldoctest zonotope_label
julia> vertices_list(Z)
4-element Array{Array{Float64,1},1}:
 [1.1, 0.1]
 [0.9, 0.1]
 [0.9, -0.1]
 [1.1, -0.1]
```

The support vector along a given direction can be computed using `σ`
(resp. the support function can be computed using `ρ`):

```jldoctest zonotope_label
julia> σ([1., 1.], Z)
2-element Array{Float64,1}:
 1.1
 0.1
```

Zonotopes admit an alternative constructor that receives a list of
vectors, each vector representing a generator:

```jldoctest
julia> Z = Zonotope(ones(2), [[1., 0.], [0., 1.], [1., 1.]])
Zonotope{Float64,Array{Float64,1},Array{Float64,2}}([1.0, 1.0], [1.0 0.0 1.0; 0.0 1.0 1.0])

julia> genmat(Z)
2×3 Array{Float64,2}:
 1.0  0.0  1.0
 0.0  1.0  1.0
```
"""
struct Zonotope{N<:Real, VN<:AbstractVector{N}, MN<:AbstractMatrix{N}} <: AbstractZonotope{N}
    center::VN
    generators::MN

    function Zonotope(center::VN, generators::MN) where {N<:Real,
                                                         VN<:AbstractVector{N},
                                                         MN<:AbstractMatrix{N}}
        @assert length(center) == size(generators, 1) "the dimension of the " *
            "center ($(length(center))) and the generators " *
            "($(size(generators, 1))) need to match"
        return new{N, VN, MN}(center, generators)
    end
end

isoperationtype(::Type{<:Zonotope}) = false
isconvextype(::Type{<:Zonotope}) = true

# constructor from center and list of generators
function Zonotope(center::VN, generators_list::AbstractVector{VN}) where {N<:Real, VN<:AbstractVector{N}}
    MT = _matrix_type(VN)
    G = MT(undef, length(center), length(generators_list))
    for (j, gj) in enumerate(generators_list)
        @inbounds G[:, j] = gj
    end
    return Zonotope(center, G)
end

"""
    remove_zero_generators(Z::Zonotope{N, VN, MN}) where {N<:Real,
                                                          VN<:AbstractVector{N},
                                                          MN<:AbstractMatrix{N}}

Return a new zonotope removing the generators which are zero of the given zonotope.

### Input

- `Z` -- zonotope

### Output

If there are no zero generators, the result is the original zonotope `Z`.
Otherwise the result is a new zonotope that has the center and generators as `Z`
except for those generators that are zero.
"""
function remove_zero_generators(Z::Zonotope{N, VN, MN}) where {N<:Real,
                                                               VN<:AbstractVector{N},
                                                               MN<:AbstractMatrix{N}}
    G = Z.generators
    G2 = remove_zero_columns(G)
    if G === G2
        return Z
    end
    return Zonotope(Z.center, G2)
end

# --- AbstractCentrallySymmetric interface functions ---

"""
    center(Z::Zonotope{N}) where {N<:Real}

Return the center of a zonotope.

### Input

- `Z` -- zonotope

### Output

The center of the zonotope.
"""
function center(Z::Zonotope{N}) where {N<:Real}
    return Z.center
end

# --- AbstractZonotope interface functions ---

"""
    togrep(Z::Zonotope)

Return a generator representation of a zonotope.

### Input

- `Z` -- zonotopic set

### Output

The same set `Z`.
"""
function togrep(Z::Zonotope)
    return Z
end

"""
   genmat(Z::Zonotope)

Return the generator matrix of a zonotope.

### Input

- `Z` -- zonotope

### Output

A matrix where each column represents one generator of the zonotope `Z`.
"""
function genmat(Z::Zonotope)
    return Z.generators
end

"""
    generators(Z::Zonotope)

Return an iterator over the generators of a zonotope.

### Input

- `Z` -- zonotope

### Output

An iterator over the generators of `Z`.
"""
function generators(Z::Zonotope)
    return generators_fallback(Z)
end

"""
    ngens(Z::Zonotope)

Return the number of generators of a zonotope.

### Input

- `Z` -- zonotope

### Output

Integer representing the number of generators.
"""
ngens(Z::Zonotope) = size(Z.generators, 2)

# --- LazySet interface functions ---

"""
    rand(::Type{Zonotope}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
         [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing)

Create a random zonotope.

### Input

- `Zonotope`       -- type for dispatch
- `N`              -- (optional, default: `Float64`) numeric type
- `dim`            -- (optional, default: 2) dimension
- `rng`            -- (optional, default: `GLOBAL_RNG`) random number generator
- `seed`           -- (optional, default: `nothing`) seed for reseeding
- `num_generators` -- (optional, default: `-1`) number of generators of the
                      zonotope (see comment below)

### Output

A random zonotope.

### Algorithm

All numbers are normally distributed with mean 0 and standard deviation 1.

The number of generators can be controlled with the argument `num_generators`.
For a negative value we choose a random number in the range `dim:2*dim` (except
if `dim == 1`, in which case we only create a single generator).
"""
function rand(::Type{Zonotope};
              N::Type{<:Real}=Float64,
              dim::Int=2,
              rng::AbstractRNG=GLOBAL_RNG,
              seed::Union{Int, Nothing}=nothing,
              num_generators::Int=-1)
    rng = reseed(rng, seed)
    center = randn(rng, N, dim)
    if num_generators < 0
        num_generators = (dim == 1) ? 1 : rand(dim:2*dim)
    end
    generators = randn(rng, N, dim, num_generators)
    return Zonotope(center, generators)
end

# --- Zonotope functions ---

"""
    scale(α::Real, Z::Zonotope)

Concrete scaling of a zonotope.

### Input

- `α` -- scalar
- `Z` -- zonotope

### Output

The zonotope obtained by applying the numerical scale to the center and
generators of ``Z``.
"""
function scale(α::Real, Z::Zonotope)
    c = α .* Z.center
    gi = α .* Z.generators
    return Zonotope(c, gi)
end

"""
    scale!(α::Real, Z::Zonotope)

Concrete scaling of a zonotope modifing `Z` in-place

### Input

- `α` -- scalar
- `Z` -- zonotope

### Output

The zonotope `Z` after applying the numerical scale `α` to its center and generators.
"""
function scale!(α::Real, Z::Zonotope)
    c = Z.center
    G = Z.generators
    c .= α .* c
    G .= α .* G
    return Z
end

"""
    reduce_order(Z::Zonotope, r::Union{Integer, Rational})

Reduce the order of a zonotope by overapproximating with a zonotope with less
generators.

### Input

- `Z` -- zonotope
- `r` -- desired order

### Output

A new zonotope with less generators, if possible.

### Algorithm

See `overapproximate(Z::Zonotope{N}, ::Type{<:Zonotope}, r::Union{Integer, Rational}) where {N<:Real}` for details.
"""
function reduce_order(Z::Zonotope{N}, r::Union{Integer, Rational}) where {N<:Real}
    return overapproximate(Z, Zonotope, r)
end

# ============================
# Zonotope splitting methods
# ============================

function _split(Z::Zonotope, j::Int)
    c, G = Z.center, Z.generators

    c₁ = similar(c)
    G₁ = similar(G)
    Z₁ = Zonotope(c₁, G₁)

    c₂ = similar(c)
    G₂ = similar(G)
    Z₂ = Zonotope(c₂, G₂)

    split!(Z₁, Z₂, Z, j)
end

function split!(Z₁::Zonotope, Z₂::Zonotope, Z::Zonotope, j::Int)
    c, G = Z.center, Z.generators
    n, p = size(G)
    @assert 1 <= j <= p "cannot split a zonotope with $p generators along index $j"

    c₁, G₁ = Z₁.center, Z₁.generators
    c₂, G₂ = Z₂.center, Z₂.generators
    copyto!(G₁, G)

    @inbounds for i in 1:n
        α = G[i, j] / 2
        c₁[i] = c[i] - α
        c₂[i] = c[i] + α
        G₁[i, j] = α
    end
    copyto!(G₂, G₁)

    return _split_ret(Z₁, Z₂)
end

_split_ret(Z₁::Zonotope, Z₂::Zonotope) = (Z₁, Z₂)

function load_split_static()
return quote

function _split_ret(Z₁::Zonotope{N, SV, SM}, Z₂::Zonotope{N, SV, SM}) where {N, n, p, SV<:MVector{n, N}, SM<:MMatrix{n, p, N}}
    Z₁ = Zonotope(SVector{n}(Z₁.center), SMatrix{n, p}(Z₁.generators))
    Z₂ = Zonotope(SVector{n}(Z₂.center), SMatrix{n, p}(Z₂.generators))
    return Z₁, Z₂
end

end end  # quote / load_split_static

function _split(Z::Zonotope, gens::AbstractVector, n::AbstractVector)
    p = length(gens)
    @assert p == length(n) "the number of generators doesn't match the " *
                            "number of indicated partitions ($p and $(length(n)))"

    @assert p <= ngens(Z) "the number of generators to split is greater " *
                          "than the number of generators of the zonotope ($p and $(ngens(Z)))"

    Zs = [Z]
    for (i, g) in enumerate(gens)
        for j = 1:n[i]
            km = length(Zs)
            for k = 1:km
                append!(Zs, split(Zs[k], g))
            end
            deleteat!(Zs, 1:km)
        end
    end
    return Zs
end

"""
    linear_map!(Zout::Zonotope, M::AbstractMatrix, Z::Zonotope)

Compute the concrete linear map of a zonotope storing the result in `Zout`.

### Input

- `Zout` -- zonotope (output)
- `M`    -- matrix
- `Z`    -- zonotope

### Output

The zonotope `Zout`, which is modified in-place.
"""
function linear_map!(Zout::Zonotope, M::AbstractMatrix, Z::Zonotope)
    mul!(Zout.center, M, Z.center)
    mul!(Zout.generators, M, Z.generators)
    return Zout
end

"""
    quadratic_map(Q::Vector{MT}, Z::Zonotope{N}) where {N, MT<:AbstractMatrix{N}}

Return an overapproximation of the quadratic map of the given zonotope.

### Input

- `Z` -- zonotope
- `Q` -- array of square matrices

### Output

An overapproximation of the quadratic map of the given zonotope.

### Notes

Mathematically, a quadratic map of a zonotope is defined as:

```math
Z_Q = \\right\\{ \\lambda | \\lambda_i = x^T Q\\^{(i)} x,~i = 1, \\ldots, n,~x \\in Z \\left\\}
```
such that each coordinate ``i`` of the resulting zonotope is influenced by ``Q\\^{(i)}``

### Algorithm

This function implements [Lemma 1, 1].

[1] *Matthias Althoff and Bruce H. Krogh. 2012. Avoiding geometric intersection
operations in reachability analysis of hybrid systems. In Proceedings of the
15th ACM international conference on Hybrid Systems: Computation and Control
(HSCC ’12). Association for Computing Machinery, New York, NY, USA, 45–54.*
"""
function quadratic_map(Q::Vector{MT}, Z::Zonotope{N}) where {N, MT<:AbstractMatrix{N}}
    @assert length(Q) == dim(Z) "the number of matrices needs to match the dimension of the zonotope"
    G = genmat(Z)
    c = center(Z)
    n, p = size(G)
    h = Matrix{N}(undef, n, binomial(p+2, 2)-1)
    d = Vector{N}(undef, n)
    g(x) = view(G, :, x)
    cᵀ = c'
    for (i, Qᵢ) in enumerate(Q)
        cᵀQᵢ = cᵀ * Qᵢ
        Qᵢc = Qᵢ * c
        aux = zero(N)
        for j=1:p
            aux += g(j)' * Qᵢ * g(j)
            h[i, j] = cᵀQᵢ * g(j) + g(j)' * Qᵢc
            h[i, p+j] = 0.5 * g(j)' * Qᵢ * g(j)
        end
        d[i] = cᵀQᵢ * c + 0.5 * aux
        l = 0
        for j=1:p-1
            gjᵀQᵢ = g(j)' * Qᵢ
            Qᵢgj = Qᵢ * g(j)
            for k=j+1:p
                l += 1
                h[i, 2p+l] = gjᵀQᵢ * g(k) + g(k)' * Qᵢgj
            end
        end
    end
    return Zonotope(d, remove_zero_columns(h))
end

"""
    _bound_intersect_2D(Z::Zonotope, L::Line2D)

Return the support function in the direction [0, 1] of the intersection between
the given zonotope and line.

### Input

- `Z` -- zonotope
- `L` -- vertical line 2D

### Output

The support function in the direction [0, 1] of the intersection between the
given zonotope and line.

### Notes

The algorithm assumes that the given line is vertical and that the intersection
between the given sets is not empty.

### Algorithm

This function implements [Algorithm 8.2, 1].

[1] *Colas Le Guernic. Reachability Analysis of Hybrid Systems with Linear
Continuous Dynamics. Computer Science [cs]. Université Joseph-Fourier - Grenoble
I, 2009. English. fftel-00422569v2f*
"""
function _bound_intersect_2D(Z::Zonotope, L::Line2D)
    c = center(Z)
    P = copy(c)
    G = genmat(Z)
    r = ngens(Z)
    g(x) = view(G, :, x)
    for i = 1:r
        gi = g(i)
        if !_isupwards(gi)
            gi .= -gi
        end
        P .= P - gi
    end
    G = sortslices(G, dims=2, by=x->atan(x[2], x[1])) # sort gens
    if P[1] < L.b
        G .= G[:,end:-1:1]
    end
    j = 1
    while isdisjoint(LineSegment(P, P+2g(j)), L)
        P .= P + 2g(j)
        j += 1
        if j > size(G, 2)
            error("Got unexpected error, check that the sets intersect")
        end
    end
    singleton = intersection(LineSegment(P, P+2g(j)), L)
    return element(singleton)[2]
end
# ====================================
# Zonotope vertex enumeration methods
# ====================================

function _vertices_list_2D(c::AbstractVector{N}, G::AbstractMatrix{N}; apply_convex_hull::Bool) where {N}
    if same_sign(G)
        return _vertices_list_2D_positive(c, G)
    else
        # FIXME generalized 2D vertices list function is not implemented yet
        # See LazySets#2209
        return _vertices_list_iterative(c, G, apply_convex_hull=apply_convex_hull)
    end
end

function _vertices_list_2D_positive(c::AbstractVector{N}, G::AbstractMatrix{N}) where {N}
    n, p = size(G)

    # TODO special case p = 1 or p = 2 ?

    sorted_G = sortslices(G, dims=2, by=x->atan(x[2], x[1]))
    index = ones(N, p, 2*p)
    @inbounds for i in 1:p
        index[i, i+1:i+p-1] .= -one(N)
    end
    index[:, 1] .= -one(N)
    V = sorted_G * index .+ c
    return [V[:, i] for i in 1:2*p]
end

function _vertices_list_iterative(c::VN, G::MN; apply_convex_hull::Bool) where {N, VN<:AbstractVector{N}, MN<:AbstractMatrix{N}}
    p = size(G, 2)
    vlist = Vector{VN}()
    sizehint!(vlist, 2^p)

    for ξi in Iterators.product([(1, -1) for i = 1:p]...)
        push!(vlist, c .+ G * collect(ξi))
    end

    return apply_convex_hull ? convex_hull!(vlist) : vlist
end

# special case 2D zonotope of order 1/2
function _vertices_list_2D_order_one_half(c::VN, G::MN; apply_convex_hull::Bool) where {N, VN<:AbstractVector{N}, MN}
    vlist = Vector{VN}(undef, 2)
    g = view(G, :, 1)
    @inbounds begin
        vlist[1] = c .+ g
        vlist[2] = c .- g
    end
    return apply_convex_hull ? _two_points_2d!(vlist) : vlist
end

# special case 2D zonotope of order 1
function _vertices_list_2D_order_one(c::VN, G::MN; apply_convex_hull::Bool) where {N, VN<:AbstractVector{N}, MN}
    vlist = Vector{VN}(undef, 4)
    a = [one(N), one(N)]
    b = [one(N), -one(N)]
    @inbounds begin
        vlist[1] = c .+ G * a
        vlist[2] = c .- G * a
        vlist[3] = c .+ G * b
        vlist[4] = c .- G * b
    end
    return apply_convex_hull ? _four_points_2d!(vlist) : vlist
end


text:
module SOAPlite

using ASE
using JuLIP
using PyCall

py_soaplite = pyimport("soaplite")

export soaplite

function soaplite(at, cutoff; n_max=1, l_max=1, sigma=0.5, cross_over=true)

all_unique_nums = unique(collect(at.Z))

#Setting SOAPlite Descriptors
alphas, betas = py_soaplite.genBasis.getBasisFunc(cutoff, n_max)
soaplite_desc = py_soaplite.get_soap_structure

# Calculate descriptor
atom_struct = ASEAtoms(at)
soaplite_rtn_desc = soaplite_desc(atom_struct.po, alphas, betas,
                                  rCut=cutoff, nMax=n_max,
                                  Lmax=l_max, crossOver=cross_over,
                                  all_atomtypes=all_unique_nums)

return soaplite_rtn_desc
end

end # module

ext: jl


text:
using GeometricIntegrators
using ChargedParticleDynamics.GuidingCenter4d.TokamakFastSurface
using GeometricExamples

const Δt     = 10.
const ntime  = 2500
const nx     = 200
const ny     = 200
const nsave  = 5
const nplot  = 5
const run_id = "poincare_2nd_fast_vprk_dt10"

include("guiding_center_4d_settings_vprk.jl")

tableau_list = get_tableau_list_vprk_projection()

pinv = guiding_center_4d_iode_poincare_invariant_2nd(Δt, nx, ny, ntime, nsave)

include("guiding_center_4d_fast_poincare_invariant_2nd.jl")


text:
"""
    AbstractInverseRetractionMethod

Abstract type for methods for inverting a retraction (see [`inverse_retract`](@ref)).
"""
abstract type AbstractInverseRetractionMethod end

"""
    AbstractRetractionMethod

Abstract type for methods for [`retract`](@ref)ing a tangent vector to a manifold.
"""
abstract type AbstractRetractionMethod end

"""
    ApproximateInverseRetraction <: AbstractInverseRetractionMethod

An abstract type for representing approximate inverse retraction methods.

"""
abstract type ApproximateInverseRetraction <: AbstractInverseRetractionMethod end

"""
    ApproximateRetraction <: AbstractInverseRetractionMethod

An abstract type for representing approximate retraction methods.
"""
abstract type ApproximateRetraction <: AbstractRetractionMethod end

@doc raw"""
    EmbeddedRetraction{T<:AbstractRetractionMethod} <: AbstractRetractionMethod

Compute a retraction by using the retraction of type `T` in the embedding and projecting the result.

# Constructor

    EmbeddedRetraction(r::AbstractRetractionMethod)

Generate the retraction with retraction `r` to use in the embedding.
"""
struct EmbeddedRetraction{T<:AbstractRetractionMethod} <: AbstractRetractionMethod
    retraction::T
end

"""
    ExponentialRetraction <: AbstractRetractionMethod

Retraction using the exponential map.
"""
struct ExponentialRetraction <: AbstractRetractionMethod end

@doc raw"""
    ODEExponentialRetraction{T<:AbstractRetractionMethod, B<:AbstractBasis} <: AbstractRetractionMethod

Approximate the exponential map on the manifold by evaluating the ODE descripting the geodesic at 1,
assuming the default connection of the given manifold by solving the ordinary differential
equation

```math
\frac{d^2}{dt^2} p^k + Γ^k_{ij} \frac{d}{dt} p_i \frac{d}{dt} p_j = 0,
```

where ``Γ^k_{ij}`` are the Christoffel symbols of the second kind, and
the Einstein summation convention is assumed.

# Constructor

    ODEExponentialRetraction(
        r::AbstractRetractionMethod,
        b::AbstractBasis=DefaultOrthogonalBasis(),
    )

Generate the retraction with a retraction to use internally (for some approaches)
and a basis for the tangent space(s).
"""
struct ODEExponentialRetraction{T<:AbstractRetractionMethod,B<:AbstractBasis} <:
       AbstractRetractionMethod
    retraction::T
    basis::B
end
function ODEExponentialRetraction(r::T) where {T<:AbstractRetractionMethod}
    return ODEExponentialRetraction(r, DefaultOrthonormalBasis())
end
function ODEExponentialRetraction(::T, b::CachedBasis) where {T<:AbstractRetractionMethod}
    return throw(
        DomainError(
            b,
            "Cached Bases are currently not supported, since the basis has to be implemented in a surrounding of the start point as well.",
        ),
    )
end
function ODEExponentialRetraction(r::ExponentialRetraction, ::AbstractBasis)
    return throw(
        DomainError(
            r,
            "You can not use the exponential map as an inner method to solve the ode for the exponential map.",
        ),
    )
end
function ODEExponentialRetraction(r::ExponentialRetraction, ::CachedBasis)
    return throw(
        DomainError(
            r,
            "Neither the exponential map nor a Cached Basis can be used with this retraction type.",
        ),
    )
end

"""
    PolarRetraction <: AbstractRetractionMethod

Retractions that are based on singular value decompositions of the matrix / matrices
for point and tangent vector on a [`AbstractManifold`](@ref)
"""
struct PolarRetraction <: AbstractRetractionMethod end

"""
    ProjectionRetraction <: AbstractRetractionMethod

Retractions that are based on projection and usually addition in the embedding.
"""
struct ProjectionRetraction <: AbstractRetractionMethod end

"""
    QRRetraction <: AbstractRetractionMethod

Retractions that are based on a QR decomposition of the
matrix / matrices for point and tangent vector on a [`AbstractManifold`](@ref)
"""
struct QRRetraction <: AbstractRetractionMethod end

"""
    SoftmaxRetraction <: AbstractRetractionMethod

Describes a retraction that is based on the softmax function.
"""
struct SoftmaxRetraction <: AbstractRetractionMethod end

@doc raw"""
    PadeRetraction{m} <: AbstractRetractionMethod

A retraction based on the Padé approximation of order $m$
"""
struct PadeRetraction{m} <: AbstractRetractionMethod end

function PadeRetraction(m::Int)
    (m < 1) && error(
        "The Padé based retraction is only available for positive orders, not for order $m.",
    )
    return PadeRetraction{m}()
end
@doc raw"""
    CayleyRetraction <: AbstractRetractionMethod

A retraction based on the Cayley transform, which is realized by using the
[`PadeRetraction`](@ref)`{1}`.
"""
const CayleyRetraction = PadeRetraction{1}

@doc raw"""
   EmbeddedInverseRetraction{T<:AbstractInverseRetractionMethod} <: AbstractInverseRetractionMethod

Compute an inverse retraction by using the inverse retraction of type `T` in the embedding and projecting the result

# Constructor

    EmbeddedInverseRetraction(r::AbstractInverseRetractionMethod)

Generate the inverse retraction with inverse retraction `r` to use in the embedding.
"""
struct EmbeddedInverseRetraction{T<:AbstractInverseRetractionMethod} <:
       AbstractInverseRetractionMethod
    inverse_retraction::T
end

"""
    LogarithmicInverseRetraction <: AbstractInverseRetractionMethod

Inverse retraction using the [`log`](@ref)arithmic map.
"""
struct LogarithmicInverseRetraction <: AbstractInverseRetractionMethod end

@doc raw"""
    PadeInverseRetraction{m} <: AbstractRetractionMethod

An inverse retraction based on the Padé approximation of order $m$ for the retraction.
"""
struct PadeInverseRetraction{m} <: AbstractInverseRetractionMethod end

function PadeInverseRetraction(m::Int)
    (m < 1) && error(
        "The Padé based inverse retraction is only available for positive orders, not for order $m.",
    )
    return PadeInverseRetraction{m}()
end
@doc raw"""
    CayleyInverseRetraction <: AbstractInverseRetractionMethod

A retraction based on the Cayley transform, which is realized by using the
[`PadeRetraction`](@ref)`{1}`.
"""
const CayleyInverseRetraction = PadeInverseRetraction{1}

"""
    PolarInverseRetraction <: AbstractInverseRetractionMethod

Inverse retractions that are based on a singular value decomposition of the
matrix / matrices for point and tangent vector on a [`AbstractManifold`](@ref)
"""
struct PolarInverseRetraction <: AbstractInverseRetractionMethod end

"""
    ProjectionInverseRetraction <: AbstractInverseRetractionMethod

Inverse retractions that are based on a projection (or its inversion).
"""
struct ProjectionInverseRetraction <: AbstractInverseRetractionMethod end

"""
    QRInverseRetraction <: AbstractInverseRetractionMethod

Inverse retractions that are based on a QR decomposition of the
matrix / matrices for point and tangent vector on a [`AbstractManifold`](@ref)
"""
struct QRInverseRetraction <: AbstractInverseRetractionMethod end

"""
    NLSolveInverseRetraction{T<:AbstractRetractionMethod,TV,TK} <:
        ApproximateInverseRetraction

An inverse retraction method for approximating the inverse of a retraction using `NLsolve`.

# Constructor

    NLSolveInverseRetraction(
        method::AbstractRetractionMethod[, X0];
        project_tangent=false,
        project_point=false,
        nlsolve_kwargs...,
    )

Constructs an approximate inverse retraction for the retraction `method` with initial guess
`X0`, defaulting to the zero vector. If `project_tangent` is `true`, then the tangent
vector is projected before the retraction using `project`. If `project_point` is `true`,
then the resulting point is projected after the retraction. `nlsolve_kwargs` are keyword
arguments passed to `NLsolve.nlsolve`.
"""
struct NLSolveInverseRetraction{TR<:AbstractRetractionMethod,TV,TK} <:
       ApproximateInverseRetraction
    retraction::TR
    X0::TV
    project_tangent::Bool
    project_point::Bool
    nlsolve_kwargs::TK
    function NLSolveInverseRetraction(m, X0, project_point, project_tangent, nlsolve_kwargs)
        return new{typeof(m),typeof(X0),typeof(nlsolve_kwargs)}(
            m,
            X0,
            project_point,
            project_tangent,
            nlsolve_kwargs,
        )
    end
end
function NLSolveInverseRetraction(
    m,
    X0 = nothing;
    project_tangent::Bool = false,
    project_point::Bool = false,
    nlsolve_kwargs...,
)
    return NLSolveInverseRetraction(m, X0, project_point, project_tangent, nlsolve_kwargs)
end

"""
    SoftmaxInverseRetraction <: AbstractInverseRetractionMethod

Describes an inverse retraction that is based on the softmax function.
"""
struct SoftmaxInverseRetraction <: AbstractInverseRetractionMethod end

"""
    default_inverse_retraction_method(M::AbstractManifold)

The [`AbstractInverseRetractionMethod`](@ref) that is used when calling
[`inverse_retract`](@ref) without specifying the inverse retraction method.
By default, this is the [`LogarithmicInverseRetraction`](@ref).
"""
default_inverse_retraction_method(::AbstractManifold) = LogarithmicInverseRetraction()

"""
    default_retraction_method(M::AbstractManifold)

The [`AbstractRetractionMethod`](@ref) that is used when calling [`retract`](@ref) without
specifying the retraction method. By default, this is the [`ExponentialRetraction`](@ref).
"""
default_retraction_method(::AbstractManifold) = ExponentialRetraction()

"""
    inverse_retract!(M::AbstractManifold, X, p, q[, method::AbstractInverseRetractionMethod])

Compute the inverse retraction, a cheaper, approximate version of the
[`log`](@ref)arithmic map), of points `p` and `q` on the [`AbstractManifold`](@ref) `M`.
Result is saved to `X`.

Inverse retraction method can be specified by the last argument, defaulting to
[`default_inverse_retraction_method`](@ref)`(M)`. See the documentation of respective manifolds for
available methods.

See also [`retract!`](@ref).
"""
function inverse_retract!(
    M::AbstractManifold,
    X,
    p,
    q,
    m::AbstractInverseRetractionMethod = default_inverse_retraction_method(M),
)
    return _inverse_retract!(M, X, p, q, m)
end
function _inverse_retract!(M::AbstractManifold, X, p, q, ::LogarithmicInverseRetraction)
    return log!(M, X, p, q)
end

#
# dispatch to lower level
function _inverse_retract!(M::AbstractManifold, X, p, q, ::CayleyInverseRetraction)
    return inverse_retract_caley!(M, X, p, q)
end
function _inverse_retract!(M::AbstractManifold, X, p, q, m::EmbeddedInverseRetraction)
    return inverse_retract_embedded!(M, X, p, q, m.inverse_retraction)
end
function _inverse_retract!(
    M::AbstractManifold,
    X,
    p,
    q,
    ::PadeInverseRetraction{n},
) where {n}
    return inverse_retract_pade!(M, X, p, q, n)
end
function _inverse_retract!(M::AbstractManifold, X, p, q, ::PolarInverseRetraction)
    return inverse_retract_polar!(M, X, p, q)
end
function _inverse_retract!(M::AbstractManifold, X, p, q, ::ProjectionInverseRetraction)
    return inverse_retract_project!(M, X, p, q)
end
function _inverse_retract!(M::AbstractManifold, X, p, q, ::QRInverseRetraction)
    return inverse_retract_qr!(M, X, p, q)
end
function _inverse_retract!(M::AbstractManifold, X, p, q, ::SoftmaxInverseRetraction)
    return inverse_retract_softmax!(M, X, p, q)
end
function _inverse_retract!(M::AbstractManifold, X, p, q, m::NLSolveInverseRetraction)
    return inverse_retract_nlsolve!(M, X, p, q, m)
end
"""
    inverse_retract_embedded!(M::AbstractManifold, X, p, q, m::AbstractInverseRetractionMethod)

computes the mutating variant of the [`EmbeddedInverseRetraction`](@ref) using
the [`AbstractInverseRetractionMethod`](@ref) `m` in the embedding (see [`get_embedding`](@ref))
and projecting the result back.
"""
function inverse_retract_embedded!(
    M::AbstractManifold,
    X,
    p,
    q,
    m::AbstractInverseRetractionMethod,
)
    return project!(
        M,
        X,
        p,
        inverse_retract(
            get_embedding(M),
            embed(get_embedding(M), p),
            embed(get_embedding(M), q),
            m,
        ),
    )
end

"""
    inverse_retract_caley!(M::AbstractManifold, X, p, q)

computes the mutating variant of the [`CayleyInverseRetraction`](@ref),
which by default calls the first order [`PadeInverseRetraction`§(@ref).
"""
function inverse_retract_caley!(M::AbstractManifold, X, p, q)
    return inverse_retract_pade!(M, X, p, q, 1)
end

"""
    inverse_retract_pade!(M::AbstractManifold, p, q, n)

computes the mutating variant of the [`PadeInverseRetraction`](@ref)`(n)`,
"""
inverse_retract_pade!(M::AbstractManifold, p, q, n)

function inverse_retract_pade! end

"""
    inverse_retract_qr!(M::AbstractManifold, X, p, q)

computes the mutating variant of the [`QRInverseRetraction`](@ref).
"""
inverse_retract_qr!(M::AbstractManifold, X, p, q)

function inverse_retract_qr! end

"""
    inverse_retract_project!(M::AbstractManifold, X, p, q)

computes the mutating variant of the [`ProjectionInverseRetraction`](@ref).
"""
inverse_retract_project!(M::AbstractManifold, X, p, q)

function inverse_retract_project! end

"""
    inverse_retract_polar!(M::AbstractManifold, X, p, q)

computes the mutating variant of the [`PolarInverseRetraction`](@ref).
"""
inverse_retract_polar!(M::AbstractManifold, X, p, q)

function inverse_retract_polar! end

"""
    inverse_retract_nlsolve!(M::AbstractManifold, X, p, q, m::NLSolveInverseRetraction)

computes the mutating variant of the [`NLSolveInverseRetraction`](@ref) `m`.
"""
inverse_retract_nlsolve!(M::AbstractManifold, X, p, q, m::NLSolveInverseRetraction)

function inverse_retract_nlsolve! end

"""
    inverse_retract_softmax!(M::AbstractManifold, X, p, q)

computes the mutating variant of the [`SoftmaxInverseRetraction`](@ref).
"""
inverse_retract_softmax!(M::AbstractManifold, X, p, q)

function inverse_retract_softmax! end

"""
    inverse_retract(M::AbstractManifold, p, q)
    inverse_retract(M::AbstractManifold, p, q, method::AbstractInverseRetractionMethod

Compute the inverse retraction, a cheaper, approximate version of the
[`log`](@ref)arithmic map), of points `p` and `q` on the [`AbstractManifold`](@ref) `M`.

Inverse retraction method can be specified by the last argument, defaulting to
[`default_inverse_retraction_method`](@ref)`(M)`.
For available inverse retractions on certain manifolds see the documentation on the
corresponding manifold.

See also [`retract`](@ref).
"""
function inverse_retract(
    M::AbstractManifold,
    p,
    q,
    m::AbstractInverseRetractionMethod = default_inverse_retraction_method(M),
)
    return _inverse_retract(M, p, q, m)
end
function _inverse_retract(M::AbstractManifold, p, q, ::CayleyInverseRetraction)
    return inverse_retract_caley(M, p, q)
end
function _inverse_retract(M::AbstractManifold, p, q, m::EmbeddedInverseRetraction)
    return inverse_retract_embedded(M, p, q, m.inverse_retraction)
end
_inverse_retract(M::AbstractManifold, p, q, ::LogarithmicInverseRetraction) = log(M, p, q)
function _inverse_retract(M::AbstractManifold, p, q, m::NLSolveInverseRetraction)
    return inverse_retract_nlsolve(M, p, q, m)
end
function _inverse_retract(M::AbstractManifold, p, q, ::PadeInverseRetraction{n}) where {n}
    return inverse_retract_pade(M, p, q, n)
end
function _inverse_retract(M::AbstractManifold, p, q, ::PolarInverseRetraction)
    return inverse_retract_polar(M, p, q)
end
function _inverse_retract(M::AbstractManifold, p, q, ::ProjectionInverseRetraction)
    return inverse_retract_project(M, p, q)
end
function _inverse_retract(M::AbstractManifold, p, q, ::QRInverseRetraction)
    return inverse_retract_qr(M, p, q)
end
function _inverse_retract(M::AbstractManifold, p, q, ::SoftmaxInverseRetraction)
    return inverse_retract_softmax(M, p, q)
end
"""
    inverse_retract_embedded(M::AbstractManifold, p, q, m::AbstractInverseRetractionMethod)

computes the allocating variant of the [`EmbeddedInverseRetraction`](@ref) using
the [`AbstractInverseRetractionMethod`](@ref) `m` in the embedding (see [`get_embedding`](@ref))
and projecting the result back.
"""
function inverse_retract_embedded(
    M::AbstractManifold,
    p,
    q,
    m::AbstractInverseRetractionMethod,
)
    return project(
        M,
        p,
        inverse_retract(
            get_embedding(M),
            embed(get_embedding(M), p),
            embed(get_embedding(M), q),
            m,
        ),
    )
end

"""
    inverse_retract_caley(M::AbstractManifold, p, q)

computes the allocating variant of the [`CayleyInverseRetraction`](@ref),
which by default allocates and calls [`inverse_retract_caley!`](@ref).
"""
function inverse_retract_caley(M::AbstractManifold, p, q)
    X = allocate_result(M, inverse_retract, p, q)
    return inverse_retract_caley!(M, X, p, q)
end
"""
    inverse_retract_pade(M::AbstractManifold, p, q)

computes the allocating variant of the [`PadeInverseRetraction`](@ref)`(n)`,
which by default allocates and calls [`inverse_retract_pade!`](@ref ManifoldsBase.inverse_retract_pade!).
"""
function inverse_retract_pade(M::AbstractManifold, p, q, n)
    X = allocate_result(M, inverse_retract, p, q)
    return inverse_retract_pade!(M, X, p, q, n)
end

"""
    inverse_retract_polar(M::AbstractManifold, p, q)

computes the allocating variant of the [`PolarInverseRetraction`](@ref),
which by default allocates and calls [`inverse_retract_polar!`](@ref ManifoldsBase.inverse_retract_polar!).
"""
function inverse_retract_polar(M::AbstractManifold, p, q)
    X = allocate_result(M, inverse_retract, p, q)
    return inverse_retract_polar!(M, X, p, q)
end
"""
    inverse_retract_project(M::AbstractManifold, p, q)

computes the allocating variant of the [`ProjectionInverseRetraction`](@ref),
which by default allocates and calls [`inverse_retract_project!`](@ref ManifoldsBase.inverse_retract_project!).
"""
function inverse_retract_project(M::AbstractManifold, p, q)
    X = allocate_result(M, inverse_retract, p, q)
    return inverse_retract_project!(M, X, p, q)
end
"""
    inverse_retract_qr(M::AbstractManifold, p, q)

computes the allocating variant of the [`QRInverseRetraction`](@ref),
which by default allocates and calls [`inverse_retract_qr!`](@ref ManifoldsBase.inverse_retract_qr!).
"""
function inverse_retract_qr(M::AbstractManifold, p, q)
    X = allocate_result(M, inverse_retract, p, q)
    return inverse_retract_qr!(M, X, p, q)
end
"""
    inverse_retract_nlsolve(M::AbstractManifold, p, q, m::NLSolveInverseRetraction)

computes the allocating variant of the [`NLSolveInverseRetraction`](@ref) `m`,
which by default allocates and calls [`inverse_retract_nlsolve!`](@ref).
"""
function inverse_retract_nlsolve(M::AbstractManifold, p, q, m::NLSolveInverseRetraction)
    X = allocate_result(M, inverse_retract, p, q)
    return inverse_retract_nlsolve!(M, X, p, q, m)
end
"""
    inverse_retract_softmax(M::AbstractManifold, p, q)

computes the allocating variant of the [`SoftmaxInverseRetraction`](@ref),
which by default allocates and calls [`inverse_retract_softmax!`](@ref ManifoldsBase.inverse_retract_softmax!).
"""
function inverse_retract_softmax(M::AbstractManifold, p, q)
    X = allocate_result(M, inverse_retract, p, q)
    return inverse_retract_softmax!(M, X, p, q)
end

"""
    retract!(M::AbstractManifold, q, p, X)
    retract!(M::AbstractManifold, q, p, X, t::Real=1)
    retract!(M::AbstractManifold, q, p, X, method::AbstractRetractionMethod)
    retract!(M::AbstractManifold, q, p, X, t::Real=1, method::AbstractRetractionMethod)

Compute a retraction, a cheaper, approximate version of the [`exp`](@ref)onential map,
from `p` into direction `X`, scaled by `t`, on the [`AbstractManifold`](@ref) manifold `M`.
Result is saved to `q`.

Retraction method can be specified by the last argument, defaulting to
[`default_retraction_method`](@ref)`(M)`. See the documentation of respective manifolds for available
methods.

See [`retract`](@ref) for more details.
"""
function retract!(
    M::AbstractManifold,
    q,
    p,
    X,
    m::AbstractRetractionMethod = default_retraction_method(M),
)
    return _retract!(M, q, p, X, m)
end
function retract!(
    M::AbstractManifold,
    q,
    p,
    X,
    t::Real,
    method::AbstractRetractionMethod = default_retraction_method(M),
)
    return retract!(M, q, p, t * X, method)
end
# dispatch to lower level
_retract!(M::AbstractManifold, q, p, X, ::CayleyRetraction) = retract_caley!(M, q, p, X)
function _retract!(M::AbstractManifold, q, p, X, m::EmbeddedRetraction)
    return retract_embedded!(M, q, p, X, m.retraction)
end
_retract!(M::AbstractManifold, q, p, X, ::ExponentialRetraction) = exp!(M, q, p, X)
function _retract!(M::AbstractManifold, q, p, X, m::ODEExponentialRetraction)
    return retract_exp_ode!(M, q, p, X, m.retraction, m.basis)
end
_retract!(M::AbstractManifold, q, p, X, ::PolarRetraction) = retract_polar!(M, q, p, X)
function _retract!(M::AbstractManifold, q, p, X, ::ProjectionRetraction)
    return retract_project!(M, q, p, X)
end
_retract!(M::AbstractManifold, q, p, X, ::QRRetraction) = retract_qr!(M, q, p, X)
_retract!(M::AbstractManifold, q, p, X, ::SoftmaxRetraction) = retract_softmax!(M, q, p, X)
function _retract!(M::AbstractManifold, q, p, X, ::PadeRetraction{n}) where {n}
    return retract_pade!(M, q, p, X, n)
end

"""
    retract_embedded!(M::AbstractManifold, X, p, q, m::AbstractRetractionMethod)

computes the mutating variant of the [`EmbeddedRetraction`](@ref) using
the [`AbstractRetractionMethod`](@ref) `m` in the embedding (see [`get_embedding`](@ref))
and projecting the result back.
"""
function retract_embedded!(M::AbstractManifold, q, p, X, m::AbstractRetractionMethod)
    return project!(
        M,
        q,
        retract(
            get_embedding(M),
            embed(get_embedding(M), p),
            embed(get_embedding(M), p, X),
            m,
        ),
    )
end

"""
    retract_caley!(M::AbstractManifold, X, p, q)

computes the mutating variant of the [`CayleyRetraction`](@ref),
which by default falls back to calling the first order [`PadeRetraction`](@ref).
"""
function retract_caley!(M::AbstractManifold, q, p, X)
    return retract_pade!(M, q, p, X, 1)
end

"""
    retract_exp_ode!(M::AbstractManifold, q, p, X, m::AbstractRetractionMethod, B::AbstractBasis)

computes the mutating variant of the [`ODEExponentialRetraction`](@ref)`(m, B)`.
"""
retract_exp_ode!(
    M::AbstractManifold,
    q,
    p,
    X,
    m::AbstractRetractionMethod,
    B::AbstractBasis,
)

function retract_exp_ode! end

"""
    retract_pade!(M::AbstractManifold, q, p, n)

computes the mutating variant of the [`PadeRetraction`](@ref)`(n)`.
"""
retract_pade!(M::AbstractManifold, q, p, n)

function retract_pade! end

"""
    retract_project!(M::AbstractManifold, q, p, X)

computes the mutating variant of the [`ProjectionRetraction`](@ref).
"""
retract_project!(M::AbstractManifold, q, p, X)

function retract_project! end

"""
    retract_polar!(M::AbstractManifold, q, p, X)

computes the mutating variant of the [`PolarRetraction`](@ref).
"""
retract_polar!(M::AbstractManifold, q, p, X)

function retract_polar! end

"""
    retract_qr!(M::AbstractManifold, q, p, X)

computes the mutating variant of the [`QRRetraction`](@ref).
"""
retract_qr!(M::AbstractManifold, q, p, X)

function retract_qr! end

"""
    retract_softmax!(M::AbstractManifold, q, p, X)

computes the mutating variant of the [`SoftmaxRetraction`](@ref).
"""
retract_softmax!(M::AbstractManifold, q, p, X)

function retract_softmax! end

@doc raw"""
    retract(M::AbstractManifold, p, X, method::AbstractRetractionMethod=default_retraction_method(M))
    retract(M::AbstractManifold, p, X, t::Real=1, method::AbstractRetractionMethod=default_retraction_method(M))

Compute a retraction, a cheaper, approximate version of the [`exp`](@ref)onential map,
from `p` into direction `X`, scaled by `t`, on the [`AbstractManifold`](@ref) `M`.

A retraction ``\operatorname{retr}_p: T_p\mathcal M → \mathcal M`` is a smooth map that fulfils

1. ``\operatorname{retr}_p(0) = p``
2. ``D\operatorname{retr}_p(0): T_p\mathcal M \to T_p\mathcal M`` is the identity map,
i.e. ``D\operatorname{retr}_p(0)[X]=X`` holds for all ``X\in T_p\mathcal M``,

where ``D\operatorname{retr}_p`` denotes the differential of the retraction

The retraction is called of second order if for all ``X`` the curves ``c(t) = R_p(tX)``
have a zero acceleration at ``t=0``, i.e. ``c''(0) = 0``.

Retraction method can be specified by the last argument, defaulting to
[`default_retraction_method`](@ref)`(M)`. For further available retractions see the documentation of respective manifolds.

Locally, the retraction is invertible. For the inverse operation, see [`inverse_retract`](@ref).
"""
function retract(
    M::AbstractManifold,
    p,
    X,
    m::AbstractRetractionMethod = default_retraction_method(M),
)
    return _retract(M, p, X, m)
end
function retract(
    M::AbstractManifold,
    p,
    X,
    t::Real,
    m::AbstractRetractionMethod = default_retraction_method(M),
)
    return retract(M, p, t * X, m)
end
function _retract(M::AbstractManifold, p, X, m::EmbeddedRetraction)
    return retract_embedded(M, p, X, m.retraction)
end
_retract(M::AbstractManifold, p, X, ::ExponentialRetraction) = exp(M, p, X)
function _retract(M::AbstractManifold, p, X, m::ODEExponentialRetraction)
    return retract_exp_ode(M, p, X, m.retraction, m.basis)
end
_retract(M::AbstractManifold, p, X, ::PolarRetraction) = retract_polar(M, p, X)
_retract(M::AbstractManifold, p, X, ::ProjectionRetraction) = retract_project(M, p, X)
_retract(M::AbstractManifold, p, X, ::QRRetraction) = retract_qr(M, p, X)
_retract(M::AbstractManifold, p, X, ::SoftmaxRetraction) = retract_softmax(M, p, X)
_retract(M::AbstractManifold, p, X, ::CayleyRetraction) = retract_caley(M, p, X)
function _retract(M::AbstractManifold, p, X, ::PadeRetraction{n}) where {n}
    return retract_pade(M, p, X, n)
end
"""
    retract_embedded(M::AbstractManifold, p, X, m::AbstractRetractionMethod)

computes the allocating variant of the [`EmbeddedRetraction`](@ref) using
the [`AbstractRetractionMethod`](@ref) `m` in the embedding (see [`get_embedding`](@ref))
and projecting the result back.
"""
function retract_embedded(M::AbstractManifold, p, X, m::AbstractRetractionMethod)
    return project(
        M,
        retract(
            get_embedding(M),
            embed(get_embedding(M), p),
            embed(get_embedding(M), p, X),
            m,
        ),
    )
end
"""
    retract_polar(M::AbstractManifold, p, q)

computes the allocating variant of the [`PolarRetraction`](@ref),
which by default allocates and calls [`retract_polar!`](@ref ManifoldsBase.retract_polar!).
"""
function retract_polar(M::AbstractManifold, p, X)
    q = allocate_result(M, retract, p, X)
    return retract_polar!(M, q, p, X)
end
"""
    retract_project(M::AbstractManifold, p, q)

computes the allocating variant of the [`ProjectionRetraction`](@ref),
which by default allocates and calls [`retract_project!`](@ref ManifoldsBase.retract_project!).
"""
function retract_project(M::AbstractManifold, p, X)
    q = allocate_result(M, retract, p, X)
    return retract_project!(M, q, p, X)
end
"""
    retract_qr(M::AbstractManifold, p, q)

computes the allocating variant of the [`QRRetraction`](@ref),
which by default allocates and calls [`retract_qr!`](@ref ManifoldsBase.retract_qr!).
"""
function retract_qr(M::AbstractManifold, p, X)
    q = allocate_result(M, retract, p, X)
    return retract_qr!(M, q, p, X)
end
"""
    retract_exp_ode(M::AbstractManifold, p, q, m::AbstractRetractionMethod, B::AbstractBasis)

computes the allocating variant of the [`ODEExponentialRetraction`](@ref)`(m,B)`,
which by default allocates and calls [`retract_exp_ode!`](@ref ManifoldsBase.retract_exp_ode!).
"""
function retract_exp_ode(
    M::AbstractManifold,
    p,
    X,
    m::AbstractRetractionMethod,
    B::AbstractBasis,
)
    q = allocate_result(M, retract, p, X)
    return retract_exp_ode!(M, q, p, X, m, B)
end
"""
    retract_softmax(M::AbstractManifold, p, q)

computes the allocating variant of the [`SoftmaxRetraction`](@ref),
which by default allocates and calls [`retract_softmax!`](@ref ManifoldsBase.retract_softmax!).
"""
function retract_softmax(M::AbstractManifold, p, X)
    q = allocate_result(M, retract, p, X)
    return retract_softmax!(M, q, p, X)
end

"""
    retract_caley(M::AbstractManifold, p, q)

computes the allocating variant of the [`CayleyRetraction`](@ref),
which by default allocates and calls [`retract_caley!`](@ref ManifoldsBase.retract_caley!).
"""
function retract_caley(M::AbstractManifold, p, X)
    q = allocate_result(M, retract, p, X)
    return retract_caley!(M, q, p, X)
end
"""
    retract_pade(M::AbstractManifold, p, q)

computes the allocating variant of the [`PadeRetraction`](@ref)`(n)`,
which by default allocates and calls [`retract_pade!`](@ref ManifoldsBase.retract_pade!).
"""
function retract_pade(M::AbstractManifold, p, X, n)
    q = allocate_result(M, retract, p, X)
    return retract_pade!(M, q, p, X, n)
end

Base.show(io::IO, ::CayleyRetraction) = print(io, "CayleyRetraction()")
Base.show(io::IO, ::PadeRetraction{m}) where {m} = print(io, "PadeRetraction($(m))")


text:
"""
    prediction( sys::StateSpace, N::Integer )

Form the prediction matrix ``Γ`` that computes the states of the linear system `sys` over the horizon of length `N`
by using all the inputs across the entire horizon.

# Extended help

This matrix consists of the blocks ``A^i B`` where ``i`` is the diagonal number, forming the matrix
```math
⎡ B   •   •  •  ⋯⎤
  ⎢ AB  B   •  •  ⋯⎥
  ⎢ A²B AB  B  •  ⋯⎥
  ⎢ A³B A²B AB B  ⋯⎥
  ⎣ ⋮   ⋮   ⋮  ⋮ ⋱ ⎦
```

It is used to compute the component of the state vector for a horizon that is caused by the system inputs
during the horizon using the equation ``x = Γu``.

The full state vector for a horizon can be computed as ``x = Γu + Φx₀``, where ``Φ`` is the initial condition
propagation matrix from [`initialpropagation`](@ref).
"""
function prediction( sys::StateSpace, N::Integer )
    if( N < 1 )
        throw( DomainError( N, "Horizon length must be greater than 1" ) )
    end

    nu = ninputs( sys )
    nx = nstates( sys )
    mt = eltype( sys.A )

    mainDiag  = kron( I(N), sys.B )

    otherDiag = kron( I(N), to_matrix( mt, I, nx ) )

    D = zeros( mt, N, N )
    D[diagind( D, -1 )] .= 1

    otherDiag = otherDiag + kron( D, -sys.A )

    Γ = otherDiag\mainDiag

    return BlockArray{ mt }( Γ, [nx for i = 1:N], [nu for i = 1:N] )
end

"""
    initialpropagation( sys::StateSpace, N::Integer; ive::Bool = false )

Form the matrix ``Φ`` that propagates the initial condition ``x₀`` in the linear system `sys` across the horizon `N`.
When `ive` is true, the initial condition ``x₀`` is embedded in the resulting state vector

# Extended help

When `ive` is false, this matrix consists of a column of the blocks ``A^i`` where ``i`` is the row number, forming the matrix
```math
⎡ A  ⎤
  ⎢ A² ⎥
  ⎢ A³ ⎥
  ⎢ A⁴ ⎥
  ⎣ ⋮  ⎦
```

When `ive` is true, this matrix consists of a column of the blocks ``A^i-1`` where ``i`` is the row number, forming the matrix
```math
⎡ I  ⎤
  ⎢ A² ⎥
  ⎢ A² ⎥
  ⎢ A³ ⎥
  ⎣ ⋮  ⎦
```

It is used to compute the component of the state vector for a horizon that is caused by the intial state of the system using
the equation ``x = Φx₀``.

The full state vector for a horizon can be computed as ``x = Γu + Φx₀``, where ``Γ`` is the prediction matrix
from [`prediction`](@ref).
"""
function initialpropagation( sys::StateSpace, N::Integer; ive::Bool = false )
    if( N < 1 )
        throw( DomainError( N, "Horizon length must be greater than 1" ) )
    end

    sub = ive ? 1 : 0

    blockGen = (sys.A^(i-sub) for i=1:N)

    # The mortar function requires the passed array to have ndims=2, so reshape to
    # have an array with 2 dimensions even though there is only one column
    blocks = reshape( collect( blockGen ),  (N, 1) )

    return mortar( blocks )
end

"""
    hessian( clqr::ConstrainedTimeInvariantLQR, ::Type{PrimalFullCondensing} )

Form the complete Hessian for the constrained time-invariant LQR problem `clqr` using the fully
condensed form.
"""
function hessian( clqr::ConstrainedTimeInvariantLQR, ::Type{PrimalFullCondensing} )
    sysₖ = getsystem( clqr, prestabilized = true )
    N    = clqr.N
    nₓ   = nstates( sysₖ )
    nᵤ   = ninputs( sysₖ )

    Iₙ = 1.0*I(N)
    Γₖ = prediction( sysₖ, N )

    Q = clqr.Q
    S = clqr.S
    K = clqr.K
    R = clqr.R

    Qₖ = Q - K'*S' - S*K + K'*R*K
    Sₖ = S - K'*R

    Q̄ₖ = blockkron( Iₙ, Qₖ )
    S̄ₖ = blockkron( Iₙ, Sₖ )
    R̄  = blockkron( Iₙ, R  )

    Q̄ₖ[Block( N, N )] = clqr.P
    S̄ₖ[Block( N, N )] = zeros( eltype( clqr.S ), size( clqr.S ) ) - K'*R

    # Pass the result through the Hermitian type to cleanup numerical errors that can occur
    # making it non-symmetric (they are on the order of 10e-17)
    return Hermitian( Γₖ'*Q̄ₖ*Γₖ + Γₖ'*S̄ₖ + S̄ₖ'*Γₖ + R̄  )
end

function hessian_shifted( clqr::ConstrainedTimeInvariantLQR, ::Type{PrimalFullCondensing} )
    sysₖ = getsystem( clqr, prestabilized = true )
    mt   = eltype( sysₖ.A )
    N    = clqr.N
    N₁   = N+1
    nₓ   = nstates( sysₖ )
    nᵤ   = ninputs( sysₖ )

    Iₙ  = 1.0*I(N)
    Iₙ₁ = 1.0*I(N₁)

    # Form the prediction matrix with the initial value embedded
    Γₖ = prediction( sysₖ, N₁ )
    Γ̃ₖ = BlockArray{ mt }( zeros( mt, nₓ*N₁, nᵤ*N ), [nₓ for i = 1:N₁], [nᵤ for i = 1:N] )
    [Γ̃ₖ[Block(i, j)] = Γₖ[Block(i-1, j)] for i=2:N₁, j=1:N]

    Q̄ₖ = blockkron( Iₙ₁,  clqr.Q )
    S̄  = blockkron(  Iₙ,  clqr.S )
    R̄  = blockkron(  Iₙ,  clqr.R )
    K̄  = blockkron(  Iₙ, -clqr.K )

    Q̄ₖ[Block( N₁, N₁ )] = clqr.P
    S̄ = vcat( S̄, zeros( eltype( clqr.S ), (nₓ, N*nᵤ ) ) )
    K̄ = hcat( K̄, zeros( eltype( clqr.K ), (N*nᵤ, nₓ ) ) )

    # Pass the result through the Hermitian type to cleanup numerical errors that can occur
    # making it non-symmetric (they are on the order of 10e-17)
    return Hermitian( Γ̃ₖ'*(Q̄ₖ + K̄'*R̄*K̄ + S̄ *K̄ + K̄'*S̄' )*Γ̃ₖ + Γ̃ₖ'*(K̄' *R̄ + S̄ ) + (R̄ *K̄ + S̄')*Γ̃ₖ + R̄  )
end

"""
    linearcoefficients( clqr::ConstrainedTimeInvariantLQR, ::Type{PrimalFullCondensing} )

Form the vector of coefficients for the linear term of the fully condensed quadratic QP for the
CLQR problem defined by `clqr`.
"""
function linearcoefficients( clqr::ConstrainedTimeInvariantLQR, ::Type{PrimalFullCondensing} )
    sysₖ = getsystem( clqr, prestabilized = true )
    N    = clqr.N
    nx   = nstates( sysₖ )
    nu   = ninputs( sysₖ )

    Iₙ = 1.0*I(N)
    Γₖ = prediction( sysₖ, N )
    Φₖ = initialpropagation( sysₖ, N )

    Q̄ₖ = blockkron( Iₙ,  clqr.Qₖ )
    S̄  = blockkron( Iₙ,  clqr.S  )
    R̄  = blockkron( Iₙ,  clqr.R  )
    K̄  = blockkron( Iₙ, -clqr.K  )

    Q̄ₖ[Block( N, N )] = clqr.P
    S̄[Block( N, N )]  = zeros( eltype( clqr.S ), size( clqr.S ) )

    # Temporary product used in two places
    T = 1.0*I(N*nu) + Γₖ'*K̄'

    return ((Γₖ'*S̄ + T*R̄')*K̄ + T*S̄' + Γₖ'*Q̄ₖ')*Φₖ
end

function inequalityconstraints( clqr::ConstrainedTimeInvariantLQR, ::Type{PrimalFullCondensing} )
    sysₖ = getsystem( clqr, prestabilized = true )

    mt   = eltype( sysₖ.A )
    N    = clqr.N
    nx   = nstates( sysₖ )
    nu   = ninputs( sysₖ )
    Iₙ   = 1.0*I(N)

    # Get the condensed initial propagation matrix
    Φ̃ₖ = initialpropagation( sysₖ, N, ive = true )

    # Form the prediction matrix with the initial value embedded
    Γₖ = prediction( sysₖ, N )
    Γ̃ₖ = BlockArray{ mt }( zeros( mt, nx*N, nu*N ), [nx for i = 1:N], [nu for i = 1:N] )
    [Γ̃ₖ[Block(i, j)] = Γₖ[Block(i-1, j)] for i=2:N, j=1:N]

    # Create the component matrices
    Ē  = blockkron( Iₙ, clqr.E )
    F̄  = blockkron( Iₙ, clqr.F )
    K̄  = blockkron( Iₙ, clqr.K  )

    # Common matrix
    # The minus sign here comes from the fact we have a negative feedback controller (u = -Kx)
    T = Ē - F̄ *K̄

    # Compute the coefficient matrix of the linear inequality constraints
    G = T*Γ̃ₖ + F̄

    # Compute the initial state matrix for the RHS of the constraints
    L = -T*Φ̃ₖ;

    # Compute the constant vector for the RHS
    g = blockkron( ones( Float64, (N) ), clqr.g )

    return G, L, g
end

function equalityconstraints( ::ConstrainedTimeInvariantLQR, ::Type{PrimalFullCondensing} )
    # The fully condensed form does not have any equality constraints from the problem
    return nothing
end


text:
#
# This file is part of the Simulate.jl Julia package, MIT license
#
# Paul Bayer, 2019
#
# This is a Julia package for discrete event simulation
#

using Simulate, Random, Unitful, Test, .Threads
import Unitful: Time, ms, s, minute, hr

println(".... testing Simulate.jl .....")
@testset "clock.jl" begin
    include("test_clock.jl")
end

@testset "threads.jl" begin
    include("test_threads.jl")
end

@testset "channel 1 example" begin
    include("test_channels1.jl")
end

@testset "process.jl" begin
    include("test_process.jl")
end


text:
@testset "CSAMaxOrd" begin

  function sum_of_two_squares(a::fmpz)
    for i=1:Int(root(a,2))
      for j=1:Int(root(a,2))
        if a==i^2+j^2
          return true
        end
      end
    end
    return false
  end

  Qx,x=PolynomialRing(FlintQQ, "x")
  for b in Hecke.squarefree_up_to(100)[2:end]
    K,a=NumberField(x^2-b)
    O=maximal_order(K);
    cocval=Array{nf_elem,2}(2, 2)
    G=[Hecke.NfToNfMor(K,K,a),Hecke.NfToNfMor(K,K,-a)]
    cocval[1,1]=K(1)
    cocval[1,2]=K(1)
    cocval[2,1]=K(1)
    cocval[2,2]=K(-1)
    A=Hecke.CrossedProductAlgebra(K,G,cocval)
    if Hecke.issplit(A)
      @test sum_of_two_squares(fmpz(b))
    else
      @test !sum_of_two_squares(fmpz(b))
    end
  end

  A=Hecke.quaternion_algebra(4,36)
  @test Hecke.issplit(A)
  A=Hecke.quaternion_algebra(-1,-1)
  O=Hecke.AlgAssOrd(A, [A[i] for i=1:4])
  @test Hecke.schur_index_at_real_plc(O)==2

end


text:
# This file is part of JuliaFMI.
# Licensed under MIT: https://github.com/AnHeuermann/JuliaFMI/blob/master/LICENSE.txt

# This file contains functions to be exported for JuliaFMI package.

module JuliaFMI

# Dependecies
using Libdl         # For using dlopen, dlclose and so on
using LightXML      # For parsing XML files

export simulateFMU

function simulateFMU(pathToFMU::String)
    main(pathToFMU)
end

end # module

max_stars_repo_name: JKRT/JuliaFMI

max_issues_repo_name: JKRT/JuliaFMI

max_forks_repo_name: JKRT/JuliaFMI


text:
abstract type AbstractProjected{O,Sp,Sa} <: AbstractSampled{O,Sp,Sa} end

# For now we just remove CRS on GPU - it often contains strings
Adapt.adapt_structure(to, m::AbstractProjected) = Sampled(order(m), span(m), sampling(m))

"""
    Projected(order::Order, span, sampling, crs, mappedcrs)
    Projected(; order=Ordered(), span=AutoSpan(), sampling=Points(), crs, mappedcrs=nothing)

An [`AbstractSampled`]($DDabssampleddocs) `IndexMode` with projections attached.

Fields and behaviours are identical to [`Sampled`]($DDsampleddocs)
with the addition of `crs` and `mappedcrs` fields.

If both `crs` and `mappedcrs` fields contain CRS data (in a `GeoFormat` wrapper
from GeoFormatTypes.jl) the selector inputs and plot axes will be converted
from and to the specified `mappedcrs` projection automatically. A common use case
would be to pass `mappedcrs=EPSG(4326)` to the constructor when loading eg. a GDALarray:

```julia
GDALarray(filename; mappedcrs=EPSG(4326))
```

The underlying `crs` will be detected by GDAL.

If `mappedcrs` is not supplied (ie. `mappedcrs=nothing`), the base index will be
shown on plots, and selectors will need to use whatever format it is in.
"""
struct Projected{O<:Order,Sp<:Regular,Sa<:Sampling,PC,MC} <: AbstractProjected{O,Sp,Sa}
    order::O
    span::Sp
    sampling::Sa
    crs::PC
    mappedcrs::MC
end
Projected(; order=Ordered(), span=Regular(),
          sampling=Points(), crs, mappedcrs=nothing) =
    Projected(order, span, sampling, crs, mappedcrs)

crs(mode::Projected) = mode.crs
crs(mode::IndexMode) = nothing

mappedcrs(mode::Projected) = mode.mappedcrs
mappedcrs(mode::IndexMode) = nothing

rebuild(g::Projected, order=order(g), span=span(g),
        sampling=sampling(g), crs=crs(g), mappedcrs=mappedcrs(g)) =
    Projected(order, span, sampling, crs, mappedcrs)

"""
    Mapped(order::Order, span, sampling, crs, mappedcrs)
    Mapped(; order=Ordered(), span=AutoSpan(), sampling=Points(), crs=nothing, mappedcrs)

An [`AbstractSampled`]($DDabssampleddocs) `IndexMode`, where the dimension index has
been mapped to another projection, usually lat/lon or `EPSG(4326)`.

Fields and behaviours are identical to [`Sampled`]($DDsampleddocs) with the addition of
`crs` and `mappedcrs` fields.

The mapped dimension index will be used as for [`Sampled`]($DDsampleddocs),
but to save in another format the underlying `projectioncrs` may be used.
"""
struct Mapped{O<:Order,Sp<:Span,Sa<:Sampling,PC,MC} <: AbstractProjected{O,Sp,Sa}
    order::O
    span::Sp
    sampling::Sa
    crs::PC
    mappedcrs::MC
end
Mapped(; order=Ordered(), span=AutoSpan(), sampling=Points(),
       crs=nothing, mappedcrs) =
    Mapped(order, span, sampling, crs, mappedcrs)

crs(mode::Mapped, dim) = crs(mode)
crs(mode::Mapped) = mode.crs

mappedcrs(mode::Mapped, dim) = mappedcrs(mode)
mappedcrs(mode::Mapped) = mode.mappedcrs

rebuild(g::Mapped, order=order(g), span=span(g),
        sampling=sampling(g), crs=crs(g), mappedcrs=mappedcrs(g)) =
    Mapped(order, span, sampling, crs, mappedcrs)

"""
    convertmode(dstmode::Type{<:IndexMode}, x)

Convert the dimension mode between `Projected` and `Mapped`.
Other dimension modes pass through unchanged.

This is used to e.g. save a netcdf file to GeoTiff.
"""
convertmode(dstmode::Type{<:IndexMode}, A::AbstractArray) =
    rebuild(A, data(A), convertmode(dstmode, dims(A)))
convertmode(dstmode::Type{<:IndexMode}, dims::Tuple) =
    map(d -> convertmode(dstmode, d), dims)
convertmode(dstmode::Type{<:IndexMode}, dim::Dimension) =
    convertmode(dstmode, basetypeof(mode(dim)), dim)
# Non-projected IndexMode modess pass through
convertmode(dstmode::Type, srcmode::Type{<:IndexMode}, dim::Dimension) = dim
# AbstractProjected passes through if it's the same as dstmode
convertmode(dstmode::Type{M}, srcmode::Type{M}, dim::Dimension) where M<:AbstractProjected = dim
# Otherwise AbstractProjected needs ArchGDAL
convertmode(dstmode::Type, srcmode::Type{<:AbstractProjected}, dim::Dimension) =
    error("Load ArchGDAL.jl to convert projected dimensions")
# The rest of these methods are in reprojected.jl as they need ArchGDAL


text:
# Input and output routines

"""
    read_mineos(file) -> ::SteppedLayeredModel

Read a `SteppedLayeredModel` in Mineos tabular format, which contains:

Line 1: Title
Line 2: anisotropy flag (0=false, 1=true), referencefrequency (Hz), tabular format flag (1=true)
Line 3: number of layers, inner core boundary index, core mantle boundary index
Lines 4+: radius (m), density (kg/m³), Vpv (m/s), Vsv (m/s), Qκ, Qμ, Vph (m/s),
          Vsh (m/s), η

### References:

1. https://geodynamics.org/cig/software/mineos/mineos-manual.pdf
"""
function read_mineos(file)
    open(file, "r") do f
        name = readline(f)
        ifanis, tref, ifdeck = parse.((Int, Float64, Int), split(readline(f)))
        aniso = ifanis == 1
        ifdeck != 1 && error("File '$file' is not in tabular format (`ifdeck` is `$ifdeck`)")
        n, nic, noc = parse.(Int, split(readline(f)))
        d = readdlm(f)
        r = d[:,1]./1e3
        a = maximum(r)
        rho = d[:,2]./1e3
        if aniso
            vpv = d[:,3]./1e3
            vsv = d[:,4]./1e3
            vph = d[:,7]./1e3
            vsh = d[:,8]./1e3
            eta = d[:,9]
            voigt_vels = voigt_velocities.(vpv, vsv, vph, vsh, eta)
            vp = [v[1] for v in voigt_vels]
            vs = [v[2] for v in voigt_vels]
        else
            vp = d[:,3]./1e3
            vs = d[:,4]./1e3
            vpv = vph = vsv = vsh = eta = []
        end
        Qκ = d[:,5]
        Qμ = d[:,6]
        attenuation = !(all(isequal(0), Qκ) && all(isequal(0), Qμ))
        SteppedLayeredModel(a, n, r, vp, vs, rho, aniso, vph, vpv, vsh, vsv, eta,
                            attenuation, Qμ, Qκ)
    end
end

"""
    save_mineos(m::SteppedLayeredModel, file, freq=1.0, title="Model from EarthModels.jl")

Save an `EarthModel1D` as a tabular Mineos-format file.  Supply the reference
frequency `freq` in Hz and a `title`.
"""
function write_mineos(m::SteppedLayeredModel, file, freq=1.0, title="Model from EarthModels.jl")
    length(title) > 80 &&
        @warn("Mineos model files can have titles only 80 characters long " *
              "('$title' is $(length(title)) characters)")
    ifanis = m.aniso ? 1 : 0
    tref = freq
    ifdeck = 1
    N = m.n
    N > 350 && @warn("Mineos model files are limited to N ≤ 350 (have $N layers)")
    nic, noc = core_interface_layers(m)
    rho, vp, vs = m.density.*1e3, m.vp.*1e3, m.vs.*1e3
    if m.aniso
        vpv, vph, vsv, vsh, eta = m.vpv.*1e3, m.vph.*1e3, m.vsv.*1e3, m.vsh.*1e3, m.eta
    else
        vpv = vph = vp
        vsh = vsv = vs
        eta = ones(m.n)
    end
    if m.attenuation
        Qμ, Qκ = m.Qμ, m.Qκ
    else
        Qμ = Qκ = zeros(m.n)
    end
    open(file, "w") do f
        println(f, title)
        println(f, ifanis, " ", tref, " ", ifdeck)
        println(f, N, " ", nic, " ", noc)
        writedlm(f, [m.r*1e3 rho vpv vsv Qκ Qμ vph vsh eta])
    end
    nothing
end

"""
    core_interface_layers(m::SteppedLayeredModel) -> i_icb, i_cmb

Return the indices of the inner-core boundary, `i_icb`, and the core-mantle boundary,
`i_cmb`, for the model `m` based on the following assumptions:

1. The outer core is a liquid and has Vs = 0.
2. There is only one outer core
3. The inner core is solid and has Vs > 0.
4. There is only one inner core
5. The mantle above the outer core is solid and has Vs > 0.

These numbers refer to the layer which is the **top** of the inner/outer core.
"""
function core_interface_layers(m::SteppedLayeredModel)
    i = 1
    got_icb = false
    while i <= m.n
        if !got_icb && isapprox(m.vs[i], 0, atol=eps(eltype(m.vs)))
            got_icb = true
            i > 1 || error("Unexpectedly found Vs = 0 at centre of model")
            i_icb = i - 1
        end
        if got_icb && m.vs[i] > eps(eltype(m.vs))
            i_cmb = i - 1
            return i_icb, i_cmb
        end
        i += 1
    end
end

"""
    voigt_kappa_mu(A, C, L, N, F) -> κ, μ

Return the Voigt average bulk modulus `κ` and shear modulus `μ` given the Love parameters
of a transversely isotropic medium.  If A, C, etc., contain density, then κ and μ will be
true moduli, whilst if the Love parameters are provided as velocities, then these
are density-normalised 'moduli' with dimensions of L²⋅T².
"""
voigt_kappa_mu(A, C, L, N, F) = (4A + C + 4F - 4N)/9, (A + C - 2F + 5N + 6L)/15

"""
    voigt_velocities(vpv, vsv, vph, vsh, η) -> vp, vs

Return the Voigt average `vp` and `vs`, in the same units as the input, given four
velocities and the radial anisotropy parameter `η`.
"""
function voigt_velocities(vpv, vsv, vph, vsh, η)
    A = vph^2
    C = vpv^2
    L = vsv^2
    N = vsh^2
    F = η*(A - 2L)
    κ, μ = voigt_kappa_mu(A, C, L, N, F) # In km²/s²
    vp = sqrt(κ + 4/3*μ)
    vs = sqrt(μ)
    vp, vs
end


text:
using Documenter, Ramnet, Ramnet.Encoders, Ramnet.Utils

makedocs(
  sitename="Ramnet.jl"
)


text:
# vecnorm(x) = x./norm.(x[:,i] for i in 1:size(x,2))'
vecnorm(x::AbstractMatrix) = norm.(x[:,i] for i in 1:size(x,2))
function normc!(x)
    for i in 1:size(x,2)
        x[:,i]=x[:,i]./norm(x[:,i])
    end
end
rescale(A; dims=1) = (A .- mean(A, dims=dims)) ./ max.(std(A, dims=dims), eps())

max_stars_repo_name: jamesjcai/ScTenifold

max_issues_repo_name: jamesjcai/ScTenifold

max_forks_repo_name: jamesjcai/ScTenifold


text:
using PSCOPF

using Test
using JuMP
using Dates
using DataStructures

@testset verbose=true "test_tso_out_fo_dmo" begin

    #=
    TS: [11h]
    S: [S1]
                        bus 1
                        |
    (pilotable) prod_1_1|load
    Pmin=10, Pmax=100   |    ?
    Csta=0, Cprop=10    |
    DMO => 8h           |
    INIT : OFF          |
    PREV MARKET :       |
        ON?     ?(?)    |
        PROD    ?(?)    |
    PREV TSO :          |
        ON?     ?(?)    |
        PROD    ?(?)    |
                        |
    =#
    function create_context(TS, ech,
                            load_1,
                            market_on_s1, market_isdefinitive_on::Bool,
                            tso_on_s1, tso_isdefinitive_on::Bool,
                            )
        network = PSCOPF.Networks.Network()
        # Buses
        PSCOPF.Networks.add_new_bus!(network, "bus_1")
        # Pilotables
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "prod_1_1", PSCOPF.Networks.PILOTABLE,
                                                10., 100.,
                                                0., 10.,
                                                Dates.Second(3*60*60), Dates.Second(0))
        # Uncertainties
        uncertainties = PSCOPF.Uncertainties()
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", TS[1], "S1", load_1)
        # initial generators state : need to pay starting cost
        generators_init_state = SortedDict(
            "prod_1_1" => PSCOPF.OFF,
        )
        mode = PSCOPF.ManagementMode("mode_5mins", Dates.Minute(5))

        context = PSCOPF.PSCOPFContext(network, TS, mode,
                                        generators_init_state,
                                        uncertainties, nothing)

        definitive_on =  market_isdefinitive_on ? market_on_s1 : missing
        prod = (market_on_s1 == PSCOPF.ON) ? load_1 : 0.
        context.market_schedule = PSCOPF.Schedule(PSCOPF.Market(), ech-Minute(1), SortedDict(
            "prod_1_1" => PSCOPF.GeneratorSchedule("prod_1_1",
                SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.UncertainValue{PSCOPF.GeneratorState}(definitive_on,
                                                                                SortedDict("S1"=>market_on_s1))),
                SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.UncertainValue{Float64}(missing,
                                                                                SortedDict("S1"=>prod)))
                ),
            )
        )

        definitive_on =  tso_isdefinitive_on ? tso_on_s1 : missing
        prod = (tso_on_s1 == PSCOPF.ON) ? load_1 : 0.
        context.tso_schedule = PSCOPF.Schedule(PSCOPF.TSO(), ech-Minute(1), SortedDict(
            "prod_1_1" => PSCOPF.GeneratorSchedule("prod_1_1",
                SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.UncertainValue{PSCOPF.GeneratorState}(definitive_on,
                                                                                SortedDict("S1"=>tso_on_s1))),
                SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.UncertainValue{Float64}(missing,
                                                                                SortedDict("S1"=>prod)))
                ),
            )
        )

        return context
    end

    TS = [DateTime("2015-01-01T11:00:00")]

    #=
    ech=7h : before DMO
    next_ech = 8h : coincide with DMO
    before DMO + still have an ech to decide on => commitment firmness is FREE
    We can start the unit no matter what the previous decision was

    TS: [11h]
    S: [S1]
                        bus 1
                        |
    (pilotable) prod_1_1|load
    Pmin=10, Pmax=100   |    20.
    Csta=0, Cprop=10    |
    DMO => 8h           |
    INIT : OFF          |
    PREV MARKET :       |
        ON?     OFF
        PROD    0.
    PREV TSO :
        ON?     OFF
        PROD    0.

    The previous decision was to turn off the unit
    But now we need to satisfy 20MW of demand => we need to start the unit
    since we are before the DMO, we can do so!
    =#
    @testset "tso_can_start_unit_before_DMO" begin
        ech = DateTime("2015-01-01T07:00:00")
        next_ech = DateTime("2015-01-01T08:00:00") # corresponds to ECH-DMO

        context = create_context(TS, ech, 20.,
                                PSCOPF.OFF, false, #previous market
                                PSCOPF.OFF, false, #previous tso
                                )

        tso = PSCOPF.TSOOutFO()
        result, firmness = PSCOPF.run_step!(context, tso, ech, next_ech)

        # firmness
        expected_firmness = PSCOPF.Firmness(
                SortedDict("prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE), ),
                SortedDict("prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE), )
                )
        @test firmness == expected_firmness

        result = PSCOPF.run(tso, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)
        PSCOPF.update_tso_schedule!(context, ech, result, firmness, tso)

        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL
        # we started prod_1_1
        @test PSCOPF.ON == PSCOPF.get_commitment_value(context.tso_schedule, "prod_1_1", TS[1], "S1")
        @test 20. ≈ PSCOPF.get_prod_value(context.tso_schedule, "prod_1_1", TS[1], "S1")
        @test value(result.lol_model.p_loss_of_load["bus_1", TS[1], "S1"]) < 1e-09
    end

    #=
    ech=10h : after DMO
    after DMO, We can no longer start the unit

    TS: [11h]
    S: [S1]
                        bus 1
                        |
    (pilotable) prod_1_1|load
    Pmin=10, Pmax=100   |    20.
    Csta=0, Cprop=10    |
    DMO => 8h           |
    INIT : OFF          |
    PREV MARKET :       |
        ON?     OFF(definitive)
        PROD    0.
    PREV TSO :
        ON?     OFF(definitive)
        PROD    0.

    The previous decisions were to turn off the unit
    But now we need to satisfy 20MW of demand => we need to start the unit
    now we are past DMO, so we can no longer start the unit
    =#
    @testset "tso_cannot_start_unit_after_DMO" begin
        ech = DateTime("2015-01-01T10:00:00")

        context = create_context(TS, ech, 20.,
                                PSCOPF.OFF, true, #previous market
                                PSCOPF.OFF, true, #previous tso
                                )

        tso = PSCOPF.TSOOutFO()
        next_ech = DateTime("2015-01-01T10:30:00")
        result, firmness = PSCOPF.run_step!(context, tso, ech, next_ech)

        expected_firmness = PSCOPF.Firmness(
            SortedDict("prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.DECIDED), ),
            SortedDict("prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE), )
            )
        @test firmness == expected_firmness

        # Solution uses slack
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_HAS_SLACK
        # we could not start prod_1_1
        @test PSCOPF.OFF == PSCOPF.get_commitment_value(context.tso_schedule, "prod_1_1", TS[1], "S1")
        @test PSCOPF.get_prod_value(context.tso_schedule, "prod_1_1", TS[1], "S1") < 1e-09
        @test 20. ≈ value(result.lol_model.p_loss_of_load["bus_1", TS[1], "S1"])
    end

    #=
    ech=10h : after DMO
    after DMO, We can no longer start the unit

    TS: [11h]
    S: [S1]
                        bus 1
                        |
    (pilotable) prod_1_1|load
    Pmin=10, Pmax=100   |    20.
    Csta=0, Cprop=10    |
    DMO => 8h           |
    INIT : OFF          |
    PREV MARKET :       |
        ON?     ON(definitive)
        PROD    0.
    PREV TSO :
        ON?     OFF(definitive)
        PROD    0.

    The previous TSO decision was to turn off the unit
    The previous Market decision was to turn on the unit
    (The two decisions are not coherent, this is hypothetical and should not happen)
    Currently, we need to satisfy 20MW of demand => we need to start the unit
    now we are past DMO, so we can no longer start the unit if it's OFF
    our reference is the TSO
    => the unit appear shutdown for us
    => we cannot use it (even if market is using it)
    =#
    @testset "tso_cannot_start_unit_after_DMO_even_if_market_does" begin
        ech = DateTime("2015-01-01T10:00:00")
        context = create_context(TS, ech, 20.,
                                PSCOPF.ON, true, #previous market
                                PSCOPF.OFF, true, #previous tso : this is our reference
                                )

        tso = PSCOPF.TSOOutFO()
        @test PSCOPF.is_tso(tso.configs.REF_SCHEDULE_TYPE) #The reference for decided values is the preceding TSO schedule
        next_ech = DateTime("2015-01-01T10:30:00")
        result, firmness = PSCOPF.run_step!(context, tso, ech, next_ech)

        expected_firmness = PSCOPF.Firmness(
            SortedDict("prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.DECIDED), ),
            SortedDict("prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE), )
            )
        @test firmness == expected_firmness

        # Solution uses slack
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_HAS_SLACK
        # we could not start prod_1_1
        @test PSCOPF.OFF == PSCOPF.get_commitment_value(context.tso_schedule, "prod_1_1", TS[1], "S1")
        @test PSCOPF.get_prod_value(context.tso_schedule, "prod_1_1", TS[1], "S1") < 1e-09
        @test 20. ≈ value(result.lol_model.p_loss_of_load["bus_1", TS[1], "S1"])
    end

    #=
    ech=10h : after DMO
    after DMO, We can no longer start the unit

    TS: [11h]
    S: [S1]
                        bus 1
                        |
    (pilotable) prod_1_1|load
    Pmin=10, Pmax=100   |    20.
    Csta=0, Cprop=10    |
    DMO => 8h           |
    INIT : OFF          |
    PREV MARKET :       |
        ON?     ON(definitive)
        PROD    0.
    PREV TSO :
        ON?     OFF(definitive)
        PROD    0.

    The previous TSO decision was to turn off the unit
    The previous Market decision was to turn on the unit
    (The two decisions are not coherent, this is hypothetical and should not happen)
    Currently, we need to satisfy 20MW of demand => we need to start the unit
    now we are past DMO, so we can no longer start the unit if it's OFF
    our reference is the MARKET (not the default behaviour)
    => the unit appears started and we can use it
    =#
    @testset "dmo_constraint_is_linked_to_reference_schedule" begin
        ech = DateTime("2015-01-01T10:00:00")
        context = create_context(TS, ech, 20.,
                                PSCOPF.ON, true, #previous market : this will be set to be our reference (non default)
                                PSCOPF.OFF, true, #previous tso
                                )

        tso = PSCOPF.TSOOutFO(PSCOPF.TSOConfigs(REF_SCHEDULE_TYPE=PSCOPF.Market()))
        @test PSCOPF.is_market(tso.configs.REF_SCHEDULE_TYPE) #The reference for decided values is the preceding TSO schedule
        next_ech = DateTime("2015-01-01T10:30:00")
        result, firmness = PSCOPF.run_step!(context, tso, ech, next_ech)

        expected_firmness = PSCOPF.Firmness(
            SortedDict("prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.DECIDED), ),
            SortedDict("prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE), )
            )
        @test firmness == expected_firmness

        # Solution uses slack
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL
        # we could not start prod_1_1
        @test PSCOPF.ON == PSCOPF.get_commitment_value(context.tso_schedule, "prod_1_1", TS[1], "S1")
        @test 20. ≈ PSCOPF.get_prod_value(context.tso_schedule, "prod_1_1", TS[1], "S1")
        @test value(result.lol_model.p_loss_of_load["bus_1", TS[1], "S1"]) < 1e-09
    end

    #=
    ech=10h : after DMO
    after DMO, We can still shutdown a unit

    TS: [11h]
    S: [S1]
                        bus 1
                        |
    (pilotable) prod_1_1|load
    Pmin=10, Pmax=100   |    20.
    Csta=0, Cprop=10    |
    DMO => 8h           |
    INIT : OFF          |
    PREV MARKET :       |
        ON?     ON(definitive)
        PROD    20.
    PREV TSO :
        ON?     ON(definitive)
        PROD    20.

    The previous decisions were to turn ON the unit
    Currently, we no longer need the unit
    => we can shut it down even if we are past DMO
    =#
    @testset "tso_can_shutdown_unit_after_DMO" begin
        ech = DateTime("2015-01-01T10:00:00")
        context = create_context(TS, ech,
                                0., #No demand
                                PSCOPF.ON, true, #previous market
                                PSCOPF.ON, true, #previous tso
                                )

        tso = PSCOPF.TSOOutFO()
        next_ech = DateTime("2015-01-01T10:30:00")
        result, firmness = PSCOPF.run_step!(context, tso, ech, next_ech)

        expected_firmness = PSCOPF.Firmness(
            SortedDict("prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.DECIDED), ),
            SortedDict("prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE), )
            )
        @test firmness == expected_firmness

        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL
        # 0 demand => we shutdown the unit
        @test PSCOPF.OFF == PSCOPF.get_commitment_value(context.tso_schedule, "prod_1_1", TS[1], "S1")
        @test PSCOPF.get_prod_value(context.tso_schedule, "prod_1_1", TS[1], "S1") < 1e-09
        @test value(result.lol_model.p_loss_of_load["bus_1", TS[1], "S1"]) < 1e-09
    end

end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
function encode_string_literal(buf::IOBuffer, data::Array{UInt8}; huffman=true)
    mask::UInt8 = 0x0
    if huffman
        mask = 0x80
        data = huffman_encode_bytes(data)
    end

    h = mask | UInt8(length(data))
    write(buf, h)
    write(buf, data)
end

function encode_literal(buf::IOBuffer, header::HeaderBinary; options...)
    mask::UInt8 = 0x0
    write(buf, mask)

    encode_string_literal(buf, header[1]; options...)
    encode_string_literal(buf, header[2]; options...)
end

function encode_indexed(buf::IOBuffer, index::UInt8, header::HeaderBinary; options...)
    mask::UInt8 = 0x40
    write(buf, mask | index)

    encode_string_literal(buf, header[2]; options...)
end

function encode_totally_indexed(buf::IOBuffer, index::UInt8; options...)
    mask::UInt8 = 0x80
    write(buf, mask | index)
end

function find_totally_indexed(name::Array{UInt8, 1}, value::Array{UInt8, 1})
    for i = 1:length(STATIC_TABLE)
        if STATIC_TABLE[i][1] == name && STATIC_TABLE[i][2] == value
            return UInt8(i)
        end
    end
    return nothing
end

function encode(table::DynamicTable, headers; options...)
    buf = IOBuffer()
    for header in sort(collect(headers))
        index = find_totally_indexed(bytearr(header[1]), bytearr(header[2]))
        if index === nothing
            encode_literal(buf, (bytearr(header[1]), bytearr(header[2])); options...)
        else
            encode_totally_indexed(buf, index; options...)
        end
    end

    return take!(buf)
end


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL, LinearAlgebra
GL = ViewerGL
using SparseArrays
const Lar.Verbose = false

function arrange3Dfaces(V, copEV, copFE)
	EVs = Lar.FV2EVs(copEV, copFE) # polygonal face fragments

	triangulated_faces = Lar.triangulate2D(V, [copEV, copFE])
	FVs = convert(Array{Lar.Cells}, triangulated_faces)
	V = convert(Lar.Points,V')
	return V,FVs,EVs
end

V = [0.6540617999999998 0.2054992 0.2972308000000001; 0.7142364999999998 0.1455625 0.969203; 0.6872816200343089 0.1724107693563844 0.6681972618528844; 0.5941250999999997 0.8769965 0.36249240000000005; 0.6542998079614951 0.8170598007732168 1.0344646993560203; 0.6929229746417847 0.3843473621810608 0.9924101013070119; 0.6660927570227951 0.3906268097799347 0.66824889023436; 0.6924256885846046 0.3894687221287142 0.9924113129852241; 1.3260341 0.2707609000000001 0.2428771; 1.3862088 0.21082410000000004 0.9148494000000003; 1.2660973 0.9422582000000003 0.3081388; 1.3262719920385087 0.8823213992267827 0.9801110006439793; 0.7946185131077457 0.1533691699791747 0.9627011677876783; 0.795311803386693 0.18207414025178234 0.6679835359595829; 0.8049461769791598 0.3844586322271989 0.9823884300888783; 0.6978893364451173 0.3892254686183435 0.9924004714668809; -0.38740630000000004 0.49022260000000006 0.45363390000000015; 0.3249123 0.707347 0.5231232000000001; -0.17028190000000007 -0.08642419999999995 0.029717700000000125; 0.542036694733776 0.13070013212634024 0.09920695360173326; -0.317917 0.06630639999999999 1.0658723; 0.3944015999999999 0.28343080000000004 1.1353616; -0.10079260000000007 -0.5103404 0.6419561000000001; 0.6115259947337761 -0.2932160678736597 0.7114453536017331; 0.4729938678063347 0.07470259548304928 0.9819171231357612; 0.4673971777131105 0.08956648867512645 0.992844167328483; 0.4666328572520818 0.07498584149882567 0.9828352213039846; 0.19232632825298282 -0.12794307378859587 0.27191838727041745; 0.06793576156374279 -0.27817482626798234 0.4125651841369061; 0.11898101245114942 -0.28171847559315566 0.44353026597502554; -0.22025999018892456 -0.06114183775734239 0.4630346767698157; -0.21520721568246745 -0.07974293126956722 0.4693603419651943; -0.22039804348690156 -0.06683826697471296 0.4822487320160556; 0.5039882186318663 0.07330882745780876 0.6685562398308738; 0.4707836778272794 0.18222015069807973 0.6686268661789595; 0.4673140236772549 0.08983923498806379 0.9928443441979011; 0.7899025999999999 0.060579300000000016 0.6679888999999999; 0.46601 0.0749997 0.6686315999999999; 0.8043229999999999 0.38447249999999994 0.6679745999999999; 0.4804303999999999 0.3988929 0.6686172999999999; 0.7905452000000002 0.060564999999999994 0.9922023000000001; 0.46665270000000003 0.0749854 0.9928450000000001; 0.8049656 0.3844582 0.9921880000000001; 0.4810731000000001 0.3988786 0.9928307000000001; -0.22619069999999997 -0.07204550000000004 0.47156350000000014; -0.04998879999999989 0.08634890000000003 0.7965885000000001; -0.06779629999999998 0.21916400000000003 0.24378000000000005; 0.1084055869015021 0.3775584091796722 0.5688049026274027; 0.09883429999999994 -0.29982910000000007 0.4063673; 0.27503619999999995 -0.14143470000000002 0.7313923; 0.2572285999999999 -0.008619600000000005 0.17858379999999996; 0.4334305216010028 0.14977485494346238 0.503608715725915];
copEV = SparseArrays.sparse([1, 4, 14, 2, 5, 15, 1, 2, 7, 17, 3, 4, 19, 3, 6, 20, 5, 6, 9, 23, 7, 8, 52, 53, 8, 9, 60, 61, 10, 12, 14, 10, 13, 16, 11, 12, 19, 11, 13, 20, 15, 16, 18, 21, 17, 18, 54, 55, 21, 22, 68, 69, 22, 23, 59, 60, 24, 26, 36, 24, 27, 37, 25, 26, 38, 25, 27, 39, 28, 30, 36, 28, 31, 37, 29, 30, 38, 29, 33, 39, 32, 33, 34, 47, 31, 32, 35, 49, 71, 72, 73, 34, 35, 66, 67, 40, 41, 83, 84, 40, 42, 87, 88, 41, 42, 80, 81, 43, 44, 77, 78, 43, 45, 86, 87, 44, 45, 74, 75, 46, 47, 50, 51, 46, 48, 56, 57, 73, 48, 49, 63, 64, 50, 54, 65, 51, 56, 66, 52, 55, 68, 53, 57, 70, 58, 62, 65, 58, 63, 67, 71, 59, 62, 69, 61, 64, 70, 72, 74, 77, 86, 75, 79, 89, 76, 78, 90, 76, 79, 91, 80, 83, 88, 81, 85, 89, 82, 84, 90, 82, 85, 91], [1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 10, 10, 10, 11, 11, 11, 12, 12, 12, 13, 13, 13, 13, 14, 14, 14, 14, 15, 15, 15, 15, 16, 16, 16, 16, 17, 17, 17, 18, 18, 18, 19, 19, 19, 20, 20, 20, 21, 21, 21, 22, 22, 22, 23, 23, 23, 24, 24, 24, 25, 25, 25, 25, 26, 26, 26, 26, 26, 26, 26, 27, 27, 27, 27, 28, 28, 28, 28, 29, 29, 29, 29, 30, 30, 30, 30, 31, 31, 31, 31, 32, 32, 32, 32, 33, 33, 33, 33, 34, 34, 34, 34, 35, 35, 35, 35, 35, 36, 36, 36, 36, 37, 37, 37, 38, 38, 38, 39, 39, 39, 40, 40, 40, 41, 41, 41, 42, 42, 42, 42, 43, 43, 43, 44, 44, 44, 44, 45, 45, 45, 46, 46, 46, 47, 47, 47, 48, 48, 48, 49, 49, 49, 50, 50, 50, 51, 51, 51, 52, 52, 52], Int8[-1, -1, -1, -1, -1, -1, 1, 1, -1, -1, -1, 1, -1, 1, -1, -1, 1, 1, -1, -1, 1, -1, -1, -1, 1, 1, -1, -1, -1, -1, 1, 1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, -1, -1, 1, 1, -1, -1, 1, -1, -1, -1, 1, 1, -1, 1, -1, -1, -1, 1, -1, -1, -1, 1, -1, 1, 1, -1, -1, -1, 1, 1, -1, 1, -1, 1, 1, 1, -1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, -1, -1, -1, -1, -1, -1, 1, -1, -1, -1, 1, 1, -1, -1, -1, -1, -1, -1, 1, -1, -1, 1, 1, 1, -1, -1, -1, 1, -1, -1, 1, -1, -1, -1, 1, 1, 1, -1, -1, 1, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, -1, -1, -1, 1, -1, 1, 1, -1, 1, 1, 1, 1, -1, 1, -1, 1, 1, 1, 1, 1]);
copFE = SparseArrays.sparse([1, 4, 2, 5, 1, 6, 1, 7, 2, 8, 1, 9, 1, 2, 20, 22, 1, 2, 27, 28, 1, 2, 23, 24, 3, 4, 3, 6, 3, 7, 3, 9, 4, 7, 5, 8, 4, 9, 4, 5, 20, 22, 4, 5, 30, 31, 6, 7, 6, 9, 8, 9, 30, 31, 8, 9, 27, 29, 8, 9, 23, 24, 10, 13, 10, 14, 10, 16, 10, 18, 11, 13, 11, 14, 11, 16, 11, 18, 12, 19, 11, 18, 11, 12, 25, 26, 11, 12, 32, 34, 13, 16, 13, 18, 14, 16, 14, 18, 14, 15, 44, 45, 14, 15, 37, 38, 14, 15, 40, 41, 16, 17, 43, 44, 16, 17, 35, 36, 16, 17, 39, 40, 18, 19, 20, 21, 18, 19, 25, 26, 18, 19, 18, 19, 20, 25, 21, 26, 22, 27, 20, 28, 20, 30, 22, 31, 21, 32, 20, 33, 23, 25, 23, 29, 24, 27, 23, 28, 23, 30, 23, 23, 25, 30, 26, 32, 25, 34, 27, 31, 29, 30, 28, 33, 34, 33, 32, 33, 35, 39, 36, 40, 36, 42, 35, 43, 36, 44, 36, 46, 37, 41, 38, 40, 38, 42, 37, 45, 38, 44, 38, 46, 39, 43, 40, 44, 41, 45, 40, 46, 42, 44, 42, 46], [1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 14, 14, 15, 15, 16, 16, 17, 17, 17, 17, 18, 18, 18, 18, 19, 19, 20, 20, 21, 21, 21, 21, 22, 22, 22, 22, 23, 23, 23, 23, 24, 24, 25, 25, 26, 26, 27, 27, 28, 28, 29, 29, 30, 30, 31, 31, 32, 32, 33, 33, 34, 34, 34, 34, 35, 35, 35, 35, 36, 36, 37, 37, 38, 38, 39, 39, 40, 40, 40, 40, 41, 41, 41, 41, 42, 42, 42, 42, 43, 43, 43, 43, 44, 44, 44, 44, 45, 45, 45, 45, 46, 46, 46, 46, 47, 47, 47, 47, 48, 48, 49, 49, 50, 50, 51, 51, 52, 52, 53, 53, 54, 54, 55, 55, 56, 56, 57, 57, 58, 58, 59, 59, 60, 60, 61, 61, 62, 62, 63, 64, 65, 65, 66, 66, 67, 67, 68, 68, 69, 69, 70, 70, 71, 72, 73, 73, 74, 74, 75, 75, 76, 76, 77, 77, 78, 78, 79, 79, 80, 80, 81, 81, 82, 82, 83, 83, 84, 84, 85, 85, 86, 86, 87, 87, 88, 88, 89, 89, 90, 90, 91, 91], Int8[1, 1, -1, -1, -1, 1, -1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, -1, 1, -1, -1, 1, -1, 1, 1, -1, -1, -1, -1, -1, 1, -1, -1, -1, 1, -1, -1, 1, 1, -1, 1, -1, -1, 1, -1, 1, -1, 1, 1, 1, -1, 1, 1, -1, -1, 1, -1, 1, 1, -1, 1, -1, 1, 1, -1, 1, -1, 1, 1, 1, 1, -1, -1, -1, -1, -1, 1, -1, -1, 1, -1, 1, -1, 1, 1, -1, 1, -1, 1, -1, -1, -1, 1, -1, -1, 1, 1, 1, 1, -1, 1, -1, -1, 1, -1, 1, 1, -1, -1, 1, -1, 1, 1, -1, 1, -1, -1, 1, -1, 1, -1, 1, 1, -1, 1, -1, 1, -1, -1, 1, 1, -1, -1, 1, -1, -1, 1, 1, 1, -1, -1, 1, 1, -1, -1, 1, -1, -1, 1, 1, 1, -1, 1, 1, 1, 1, -1, -1, -1, -1, -1, 1, -1, -1, -1, 1, 1, -1, 1, -1, -1, 1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, -1, 1, 1, -1, -1, 1, 1, 1, -1, 1, 1, -1, -1, -1, 1, 1, -1, -1, 1, -1, 1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, 1]);

V,FVs,EVs = arrange3Dfaces(V, copEV, copFE);
GL.VIEW(GL.GLExplode(V,FVs,1.1,1.1,1.1,99,1));
GL.VIEW(GL.GLExplode(V,EVs,1.5,1.5,1.5,99,1));

V = [0.6540617999999998 0.2054992 0.2972308000000001; 0.7142364999999998 0.1455625 0.969203; 0.6872816200343089 0.1724107693563844 0.6681972618528844; 0.5941250999999997 0.8769965 0.36249240000000005; 0.6542998079614951 0.8170598007732168 1.0344646993560203; 0.6929229746417847 0.3843473621810608 0.9924101013070119; 0.6660927570227951 0.3906268097799347 0.66824889023436; 0.6924256885846046 0.3894687221287142 0.9924113129852241; 1.3260341 0.2707609000000001 0.2428771; 1.3862088 0.21082410000000004 0.9148494000000003; 1.2660973 0.9422582000000003 0.3081388; 1.3262719920385087 0.8823213992267827 0.9801110006439793; 0.7946185131077457 0.1533691699791747 0.9627011677876783; 0.795311803386693 0.18207414025178234 0.6679835359595829; 0.8049461769791598 0.3844586322271989 0.9823884300888783; 0.6978893364451173 0.3892254686183435 0.9924004714668809; -0.38740630000000004 0.49022260000000006 0.45363390000000015; 0.3249123 0.707347 0.5231232000000001; -0.17028190000000007 -0.08642419999999995 0.029717700000000125; 0.542036694733776 0.13070013212634024 0.09920695360173326; -0.317917 0.06630639999999999 1.0658723; 0.3944015999999999 0.28343080000000004 1.1353616; -0.10079260000000007 -0.5103404 0.6419561000000001; 0.6115259947337761 -0.2932160678736597 0.7114453536017331; 0.4729938678063347 0.07470259548304928 0.9819171231357612; 0.4673971777131105 0.08956648867512645 0.992844167328483; 0.4666328572520818 0.07498584149882567 0.9828352213039846; 0.5039882186318663 0.07330882745780876 0.6685562398308738; 0.4707836778272794 0.18222015069807973 0.6686268661789595; 0.4673140236772549 0.08983923498806379 0.9928443441979011; 0.7899025999999999 0.060579300000000016 0.6679888999999999; 0.46601 0.0749997 0.6686315999999999; 0.8043229999999999 0.38447249999999994 0.6679745999999999; 0.4804303999999999 0.3988929 0.6686172999999999; 0.7905452000000002 0.060564999999999994 0.9922023000000001; 0.46665270000000003 0.0749854 0.9928450000000001; 0.8049656 0.3844582 0.9921880000000001; 0.4810731000000001 0.3988786 0.9928307000000001];
copEV = SparseArrays.sparse([1, 4, 14, 2, 5, 15, 1, 2, 7, 17, 3, 4, 19, 3, 6, 20, 5, 6, 9, 23, 7, 8, 46, 47, 8, 9, 54, 55, 10, 12, 14, 10, 13, 16, 11, 12, 19, 11, 13, 20, 15, 16, 18, 21, 17, 18, 48, 49, 21, 22, 62, 63, 22, 23, 53, 54, 24, 26, 36, 24, 27, 37, 25, 26, 38, 25, 27, 39, 28, 30, 36, 28, 31, 37, 29, 30, 38, 29, 33, 39, 32, 33, 34, 41, 31, 32, 35, 43, 65, 66, 67, 34, 35, 60, 61, 40, 41, 44, 45, 40, 42, 50, 51, 67, 42, 43, 57, 58, 44, 48, 59, 45, 50, 60, 46, 49, 62, 47, 51, 64, 52, 56, 59, 52, 57, 61, 65, 53, 56, 63, 55, 58, 64, 66], [1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 10, 10, 10, 11, 11, 11, 12, 12, 12, 13, 13, 13, 13, 14, 14, 14, 14, 15, 15, 15, 15, 16, 16, 16, 16, 17, 17, 17, 18, 18, 18, 19, 19, 19, 20, 20, 20, 21, 21, 21, 22, 22, 22, 23, 23, 23, 24, 24, 24, 25, 25, 25, 25, 26, 26, 26, 26, 26, 26, 26, 27, 27, 27, 27, 28, 28, 28, 28, 29, 29, 29, 29, 29, 30, 30, 30, 30, 31, 31, 31, 32, 32, 32, 33, 33, 33, 34, 34, 34, 35, 35, 35, 36, 36, 36, 36, 37, 37, 37, 38, 38, 38, 38], Int8[-1, -1, -1, -1, -1, -1, 1, 1, -1, -1, -1, 1, -1, 1, -1, -1, 1, 1, -1, -1, 1, -1, -1, -1, 1, 1, -1, -1, -1, -1, 1, 1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, -1, -1, 1, 1, -1, -1, 1, -1, -1, -1, 1, 1, -1, 1, -1, -1, -1, 1, -1, -1, -1, 1, -1, 1, 1, -1, -1, -1, 1, 1, -1, 1, -1, 1, 1, 1, -1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, -1, -1, -1, 1, -1, -1, 1, -1, -1, -1, 1, 1, 1, -1, -1, 1, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]);
copFE = SparseArrays.sparse([1, 4, 2, 5, 1, 6, 1, 7, 2, 8, 1, 9, 1, 2, 18, 20, 1, 2, 25, 26, 1, 2, 21, 22, 3, 4, 3, 6, 3, 7, 3, 9, 4, 7, 5, 8, 4, 9, 4, 5, 18, 20, 4, 5, 28, 29, 6, 7, 6, 9, 8, 9, 28, 29, 8, 9, 25, 27, 8, 9, 21, 22, 10, 13, 10, 14, 10, 15, 10, 16, 11, 13, 11, 14, 11, 15, 11, 16, 12, 17, 11, 16, 11, 12, 23, 24, 11, 12, 30, 32, 13, 15, 13, 16, 14, 15, 14, 16, 16, 17, 18, 19, 16, 17, 23, 24, 16, 17, 16, 17, 18, 23, 19, 24, 20, 25, 18, 26, 18, 28, 20, 29, 19, 30, 18, 31, 21, 23, 21, 27, 22, 25, 21, 26, 21, 28, 21, 21, 23, 28, 24, 30, 23, 32, 25, 29, 27, 28, 26, 31, 32, 31, 30, 31], [1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 14, 14, 15, 15, 16, 16, 17, 17, 17, 17, 18, 18, 18, 18, 19, 19, 20, 20, 21, 21, 21, 21, 22, 22, 22, 22, 23, 23, 23, 23, 24, 24, 25, 25, 26, 26, 27, 27, 28, 28, 29, 29, 30, 30, 31, 31, 32, 32, 33, 33, 34, 34, 34, 34, 35, 35, 35, 35, 36, 36, 37, 37, 38, 38, 39, 39, 40, 40, 40, 40, 41, 41, 41, 41, 42, 42, 43, 43, 44, 44, 45, 45, 46, 46, 47, 47, 48, 48, 49, 49, 50, 50, 51, 51, 52, 52, 53, 53, 54, 54, 55, 55, 56, 56, 57, 58, 59, 59, 60, 60, 61, 61, 62, 62, 63, 63, 64, 64, 65, 66, 67, 67], Int8[1, 1, -1, -1, -1, 1, -1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, -1, 1, -1, -1, 1, -1, 1, 1, -1, -1, -1, -1, -1, 1, -1, -1, -1, 1, -1, -1, 1, 1, -1, 1, -1, -1, 1, -1, 1, -1, 1, 1, 1, -1, 1, 1, -1, -1, 1, -1, 1, 1, -1, 1, -1, 1, 1, -1, 1, -1, 1, 1, 1, 1, -1, -1, -1, -1, -1, 1, -1, -1, 1, -1, 1, -1, 1, 1, -1, 1, -1, 1, -1, -1, -1, 1, -1, -1, 1, 1, 1, 1, -1, 1, -1, 1, -1, -1, 1, 1, -1, -1, 1, -1, -1, 1, 1, 1, -1, -1, 1, 1, -1, -1, 1, -1, -1, 1, 1, 1, -1, 1, 1, 1, 1, -1, -1, -1, -1, -1, 1, -1, -1, -1, 1, 1, -1, 1, -1, -1, 1, 1, 1, 1, -1, -1, 1]);

V,FVs,EVs = arrange3Dfaces(V, copEV, copFE);
GL.VIEW(GL.GLExplode(V,FVs,1.1,1.1,1.1,99,1));
GL.VIEW(GL.GLExplode(V,EVs,1.5,1.5,1.5,99,1));


text:
@testset "complex_powers" begin
    complex_powers_comparison(z, m, T=Float64) = (
        complex_powers(Complex{T}(z), m),
        z.^collect(0:m),
        eps(T)*m
    )

    for T in [Float64, Float32, Float16]
        nozpowers = Vector{Complex{T}}(undef, 0)
        fudge = one(T) + 2 * sqrt(eps(T))
        for k in 0:25
            z = cis(k*big(π)/10)
            for m in [0, 1, 2, 3, 4, 1_000]
                mine, theirs, ϵ = complex_powers_comparison(z, m, T)
                @test mine ≈ theirs rtol=2ϵ
                @test_throws DomainError complex_powers(Complex{T}(z*fudge), m)
                @test_throws DomainError complex_powers(Complex{T}(z/fudge), m)
                inplace = zeros(Complex{T}, size(mine))
                complex_powers!(inplace, Complex{T}(z))
                @test array_equal(mine, inplace)
            end
            @test length(complex_powers!(nozpowers, Complex{T}(z))) == 0
        end
    end

end


text:
using GalerkinSparseGrids
using Base.Test
using Cubature

#---------------------------------------
# Testing vector hier DG reconstruction
#---------------------------------------

print("Testing full DG reconstruction 1-D with vector coefficients... ")

for k in 1:5
    for l in 1:6
        vect = vcoeffs_DG(1, k, l, x->sin(4*x[1]); scheme="full")
        dict = V2D(1, k, l, vect; scheme="full")
        err = x->(reconstruct_DG(dict,[x[1]])-sin(4*x[1]))^2
        @test hquadrature(err,0,1; reltol=1.0e-9, abstol=1.0e-12)[1] < 1/(1<<((l+k-1)))
    end
end

println("Test Passed.")

# Testing D2V and V2D

# Full case:

print("Testing D2V and V2D Full Case 2D... ")

for k in 1:5
    for l in 1:3
        dict = coeffs_DG(2, k, l, x->sin(4*x[1]+x[2]); scheme="full")
        vect = D2V(2, k, l, dict; scheme="full")
        @test V2D(2, k, l, vect; scheme="full")==dict
    end
end

for k in 1:5
    for l in 1:3
        vect = vcoeffs_DG(2, k, l, x->sin(4*x[1]+x[2]); scheme="full")
        dict = V2D(2, k, l, vect; scheme="full")
        @test D2V(2, k, l, dict; scheme="full")==vect
    end
end

println("Test Passed.")

# Sparse case:

print("Testing D2V and V2D Sparse Case 2D... ")

for k in 1:5
    for l in 1:3
        dict = coeffs_DG(2, k, l, x->sin(4*x[1]+x[2]); scheme="sparse")
        vect = D2V(2, k, l, dict; scheme="sparse")
        @test V2D(2, k, l, vect; scheme="sparse")==dict
    end
end

for k in 1:5
    for l in 1:3
        vect = vcoeffs_DG(2, k, l, x->sin(4*x[1]+x[2]); scheme="sparse")
        dict= V2D(2, k, l, vect; scheme="sparse")
        @test D2V(2, k, l, dict; scheme="sparse")==vect
    end
end

println("Test Passed.")


text:
module CollabUtils2RainbowBerry

using ..Ahorn, Maple

@mapdef Entity "CollabUtils2/RainbowBerry" RainbowBerry(x::Integer, y::Integer, levelSet::String="SpringCollab2020/1-Beginner")

const placements = Ahorn.PlacementDict(
    "Rainbow Berry (Collab Utils 2 / READ DOCS)" => Ahorn.EntityPlacement(
        RainbowBerry
    )
)

const sprite = "CollabUtils2/rainbowBerry/rberry0030.png"

function Ahorn.selection(entity::RainbowBerry)
    x, y = Ahorn.position(entity)
    return Ahorn.getSpriteRectangle(sprite, x, y)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::RainbowBerry, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 0, 0)

end


text:
module CauViz

using Luxor
using JuMP
using GLPK

include("GetPos.jl")
include("hierViz.jl")
include("sugiViz.jl")

export getpos, hierviz, layoutTree

end