text:
### Generic rand methods

"""
    rand(s::Sampleable)

Generate one sample for `s`.

    rand(s::Sampleable, n::Int)

Generate `n` samples from `s`. The form of the returned object depends on the variate form of `s`:

- When `s` is univariate, it returns a vector of length `n`.
- When `s` is multivariate, it returns a matrix with `n` columns.
- When `s` is matrix-variate, it returns an array, where each element is a sample matrix.
"""
rand(s::Sampleable)

"""
    rand!(s::Sampleable, A::AbstractArray)

Generate one or multiple samples from `s` to a pre-allocated array `A`. `A` should be in the
form as specified above. The rules are summarized as below:

- When `s` is univariate, `A` can be an array of arbitrary shape. Each element of `A` will
  be overriden by one sample.
- When `s` is multivariate, `A` can be a vector to store one sample, or a matrix with each
  column for a sample.
- When `s` is matrix-variate, `A` can be a matrix to store one sample, or an array of
  matrices with each element for a sample matrix.
"""
rand!(s::Sampleable, A::AbstractArray)

# univariate

function _rand!(s::Sampleable{Univariate}, A::AbstractArray)
    for i in 1:length(A)
        @inbounds A[i] = rand(s)
    end
    return A
end
rand!(s::Sampleable{Univariate}, A::AbstractArray) = _rand!(s, A)

rand(s::Sampleable{Univariate}, dims::Dims) =
    _rand!(s, Array{eltype(s)}(undef, dims))

rand(s::Sampleable{Univariate}, dims::Int...) =
    _rand!(s, Array{eltype(s)}(undef, dims))

# multivariate

function _rand!(s::Sampleable{Multivariate}, A::AbstractMatrix)
    for i = 1:size(A,2)
        _rand!(s, view(A,:,i))
    end
    return A
end

function rand!(s::Sampleable{Multivariate}, A::AbstractVector)
    length(A) == length(s) ||
        throw(DimensionMismatch("Output size inconsistent with sample length."))
    _rand!(s, A)
end

function rand!(s::Sampleable{Multivariate}, A::AbstractMatrix)
    size(A,1) == length(s) ||
        throw(DimensionMismatch("Output size inconsistent with sample length."))
    _rand!(s, A)
end

rand(s::Sampleable{Multivariate}) =
    _rand!(s, Vector{eltype(s)}(undef, length(s)))

rand(s::Sampleable{Multivariate}, n::Int) =
    _rand!(s, Matrix{eltype(s)}(undef, length(s), n))

# matrix-variate

function _rand!(s::Sampleable{Matrixvariate}, X::AbstractArray{M}) where M<:Matrix
    for i in 1:length(X)
        X[i] = rand(s)
    end
    return X
end

rand!(s::Sampleable{Matrixvariate}, X::AbstractArray{M}) where {M<:Matrix} =
    _rand!(s, X)

rand(s::Sampleable{Matrixvariate}, n::Int) =
    rand!(s, Vector{Matrix{eltype(s)}}(n))

"""
    sampler(d::Distribution) -> Sampleable

Samplers can often rely on pre-computed quantities (that are not parameters themselves) to improve efficiency.
If such a sampler exists, it can be provide with this `sampler` method, which would be used for batch sampling.
The general fallback is `sampler(d::Distribution) = d`.
"""
sampler(d::Distribution) = d


text:
using Makie

function target_function(a, b; x=π/2)
    return a*sin(b*x) + x/3
end

"""
select mode
mode should be `:wireframe` or `:surface`
"""
function draw_surface(mode)
    N = 100
    lspace = range(-2pi, stop=2pi, length=N)
    z = Float32[target_function(a,b) for a in lspace, b in lspace]
    r = range(0, stop=2pi, length=N)
    scene=eval(mode)(r, r, z)
    axis = scene[Axis]
    axis[:names][:axisnames] = ("a", "b", "a*sin(b*x)+x/3")
    scene
end

function main()
    draw_surface(:wireframe)
end

ext: jl

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
export LaurentPolynomial

"""
    LaurentPolynomial{T,X}(coeffs::AbstractVector, [m::Integer = 0], [var = :x])

A [Laurent](https://en.wikipedia.org/wiki/Laurent_polynomial) polynomial is of the form `a_{m}x^m + ... + a_{n}x^n` where `m,n` are  integers (not necessarily positive) with ` m <= n`.

The `coeffs` specify `a_{m}, a_{m-1}, ..., a_{n}`.
The argument `m` represents the lowest exponent of the variable in the series, and is taken to be zero by default.

Laurent polynomials and standard basis polynomials promote to  Laurent polynomials. Laurent polynomials may be  converted to a standard basis  polynomial when `m >= 0`
.

Integration will fail if there is a `x⁻¹` term in the polynomial.

!!! note
    `LaurentPolynomial` is not axis-aware by default, and it treats `coeffs` simply as a
    list of coefficients with the first index always corresponding to the constant term.
    In order to use the axis of `coeffs` as the exponents of the variable `var`,
    set `m` to `firstindex(coeff)` in the constructor.

# Examples:
```jldoctest laurent
julia> using Polynomials

julia> P = LaurentPolynomial
LaurentPolynomial

julia> p = P([1,1,1],  -1)
LaurentPolynomial(x⁻¹ + 1 + x)

julia> q = P([1,1,1])
LaurentPolynomial(1 + x + x²)

julia> pp = Polynomial([1,1,1])
Polynomial(1 + x + x^2)

julia> p + q
LaurentPolynomial(x⁻¹ + 2 + 2*x + x²)

julia> p * q
LaurentPolynomial(x⁻¹ + 2 + 3*x + 2*x² + x³)

julia> p * pp
LaurentPolynomial(x⁻¹ + 2 + 3*x + 2*x² + x³)

julia> pp - q
LaurentPolynomial(0)

julia> derivative(p)
LaurentPolynomial(-x⁻² + 1)

julia> integrate(q)
LaurentPolynomial(1.0*x + 0.5*x² + 0.3333333333333333*x³)

julia> integrate(p)  # x⁻¹  term is an issue
ERROR: ArgumentError: Can't integrate Laurent  polynomial with  `x⁻¹` term

julia> integrate(P([1,1,1], -5))
LaurentPolynomial(-0.25*x⁻⁴ - 0.3333333333333333*x⁻³ - 0.5*x⁻²)

julia> x⁻¹ = inv(variable(LaurentPolynomial)) # `inv` defined on monomials
LaurentPolynomial(1.0*x⁻¹)

julia> p = Polynomial([1,2,3])
Polynomial(1 + 2*x + 3*x^2)

julia> x = variable()
Polynomial(x)

julia> x^degree(p) * p(x⁻¹) # reverses  coefficients
LaurentPolynomial(3.0 + 2.0*x + 1.0*x²)
```
"""
struct LaurentPolynomial{T <: Number, X} <: StandardBasisPolynomial{T, X}
    coeffs::Vector{T}
    m::Base.RefValue{Int}
    n::Base.RefValue{Int}
    function LaurentPolynomial{T,X}(coeffs::AbstractVector{S},
                                    m::Union{Int, Nothing}=nothing) where {T <: Number, X, S}

        fnz = findfirst(!iszero, coeffs)
        fnz == nothing && return  new{T,X}(zeros(T,1), Ref(0), Ref(0))
        lnz = findlast(!iszero, coeffs)
        if Base.has_offset_axes(coeffs)
            # if present, use axes
            cs = convert(Vector{T}, coeffs[fnz:lnz])
            return new{T,X}(cs, Ref(fnz), Ref(lnz))
        else

            c = convert(Vector{T}, coeffs[fnz:lnz])

            m′ = fnz - 1 + (m == nothing ? 0 : m)
            n = m′ + (lnz-fnz)

            (n - m′ + 1  == length(c)) || throw(ArgumentError("Lengths do not match"))
            new{T,X}(c, Ref(m′),  Ref(n))
        end
    end
end

@register LaurentPolynomial

## constructors
function LaurentPolynomial{T}(coeffs::AbstractVector{S}, m::Int, var::SymbolLike=:x) where {
    T <: Number, S <: Number}
    LaurentPolynomial{T,Symbol(var)}(T.(coeffs), m)
end

function LaurentPolynomial{T}(coeffs::AbstractVector{T}, var::SymbolLike=:x) where {
    T <: Number}
    LaurentPolynomial{T, Symbol(var)}(coeffs, 0)
end

function LaurentPolynomial(coeffs::AbstractVector{T}, m::Int, var::SymbolLike=:x) where {T <: Number}
    LaurentPolynomial{T, Symbol(var)}(coeffs, m)
end

##
## conversion
##

# LaurentPolynomial is a wider collection than other standard basis polynomials.
Base.promote_rule(::Type{P},::Type{Q}) where {T, X, P <: LaurentPolynomial{T,X}, S, Q <: StandardBasisPolynomial{S, X}} = LaurentPolynomial{promote_type(T, S), X}

Base.promote_rule(::Type{Q},::Type{P}) where {T, X, P <: LaurentPolynomial{T,X}, S, Q <: StandardBasisPolynomial{S,X}} =
    LaurentPolynomial{promote_type(T, S),X}

function Base.convert(P::Type{<:Polynomial}, q::LaurentPolynomial)
    m,n = (extrema∘degreerange)(q)
    m < 0 && throw(ArgumentError("Can't convert a Laurent polynomial with m < 0"))
    P([q[i] for i  in 0:n], indeterminate(q))
end

# need to add p.m[], so abstract.jl method isn't sufficent
function Base.convert(::Type{P}, p::LaurentPolynomial) where {P<:LaurentPolynomial}
    S′ = _eltype(P)
    Y′ = _indeterminate(P)
    S = S′ == nothing ? eltype(p) : S′
    Y = Y′ == nothing ? indeterminate(p) : Y′
    isconstant(p) && return LaurentPolynomial{S,Y}(constantterm(p))
    LaurentPolynomial{S,Y}(p.coeffs, p.m[])
end

function Base.convert(::Type{P}, q::StandardBasisPolynomial{S}) where {T, P <:LaurentPolynomial{T},S}
    v′ = _indeterminate(P)
    X = v′ == nothing ? indeterminate(q) : v′
    ⟒(P){T,X}([q[i] for i in 0:degree(q)], 0)
end

##
## generic functions
##

function Base.inv(p::LaurentPolynomial{T, X}) where {T, X}
    m,n =  (extrema∘degreerange)(p)
    m != n && throw(ArgumentError("Only monomials can be inverted"))
    cs = [1/p for p in p.coeffs]
    LaurentPolynomial{eltype(cs), X}(cs, -m)
end

##
## changes to common.jl mostly as the range in the type is different
##
Base.:(==)(p1::LaurentPolynomial, p2::LaurentPolynomial) =
    check_same_variable(p1, p2) && (degreerange(chop!(p1)) == degreerange(chop!(p2))) && (coeffs(p1) == coeffs(p2))
Base.hash(p::LaurentPolynomial, h::UInt) = hash(indeterminate(p), hash(degreerange(p), hash(coeffs(p), h)))

isconstant(p::LaurentPolynomial) = iszero(lastindex(p)) && iszero(firstindex(p))

basis(P::Type{<:LaurentPolynomial{T,X}}, n::Int) where {T,X} = LaurentPolynomial{T,X}(ones(T,1), n)

# like that in common, only return zero if idx < firstindex(p)
function Base.getindex(p::LaurentPolynomial{T}, idx::Int) where {T}
    m,M = firstindex(p), lastindex(p)
    m <= idx <= M || return zero(T)
    p.coeffs[idx-m+1]
 end

# extend if out of bounds
function Base.setindex!(p::LaurentPolynomial{T}, value::Number, idx::Int) where {T}

    m,n = (extrema ∘ degreerange)(p)
    if idx > n
        append!(p.coeffs, zeros(T, idx-n))
        n = idx
        p.n[] = n
    elseif idx < m
        prepend!(p.coeffs, zeros(T, m-idx))
        m = idx
        p.m[] = m
    end

    i = idx - m + 1
    p.coeffs[i] = value

    return p

end

Base.firstindex(p::LaurentPolynomial) = p.m[]
Base.lastindex(p::LaurentPolynomial) = p.n[]
Base.eachindex(p::LaurentPolynomial) = degreerange(p)
degreerange(p::LaurentPolynomial) = firstindex(p):lastindex(p)

_convert(p::P, as) where {T,X,P <: LaurentPolynomial{T,X}} = ⟒(P)(as, firstindex(p), X)

## chop!
# trim  from *both* ends
function chop!(p::LaurentPolynomial{T};
               rtol::Real = Base.rtoldefault(real(T)),
               atol::Real = 0,) where {T}

    m0,n0 = m,n = (extrema ∘ degreerange)(p)
    for k in n:-1:m
        if isapprox(p[k], zero(T); rtol = rtol, atol = atol)
            n -= 1
        else
            break
        end
    end
    for k in m:n-1
        if isapprox(p[k], zero(T); rtol = rtol, atol = atol)
            m += 1
        else
            break
        end
    end

    cs = copy(coeffs(p))
    rng = m-m0+1:n-m0+1
    resize!(p.coeffs, length(rng))
    p.coeffs[:] = cs[rng]
    isempty(p.coeffs) && push!(p.coeffs,zero(T))
    p.m[], p.n[] = m, max(m,n)

    p

end

# use unicode exponents. XXX modify printexponent to always use these?
function showterm(io::IO, ::Type{<:LaurentPolynomial}, pj::T, var, j, first::Bool, mimetype) where {T}
    if iszero(pj) return false end
    pj = printsign(io, pj, first, mimetype)
    if !(pj == one(T) && !(showone(T) || j == 0))
        printcoefficient(io, pj, j, mimetype)
    end
    printproductsign(io, pj, j, mimetype)
    iszero(j) && return  true
    print(io, var)
    j ==1  && return  true
    unicode_exponent(io, j)
    return true
end

##
## ---- Conjugation has different defintions
##

"""
    conj(p)

This satisfies `conj(p(x)) = conj(p)(conj(x)) = p̄(conj(x))` or `p̄(x) = (conj ∘ p ∘ conj)(x)`

Examples

```jldoctest laurent
julia> using Polynomials;

julia> z = variable(LaurentPolynomial, :z)
LaurentPolynomial(1.0*z)

julia> p = LaurentPolynomial([im, 1+im, 2 + im], -1, :z)
LaurentPolynomial(im*z⁻¹ + 1 + im + (2 + im)z)

julia> conj(p)(conj(z)) ≈ conj(p(z))
true

julia> conj(p)(z) ≈ (conj ∘ p ∘ conj)(z)
true
```
"""
LinearAlgebra.conj(p::P) where {P <: LaurentPolynomial} = map(conj, p)

"""
    paraconj(p)

[cf.](https://ccrma.stanford.edu/~jos/filters/Paraunitary_FiltersC_3.html)

Call `p̂ = paraconj(p)` and `p̄` = conj(p)`, then this satisfies
`conj(p(z)) = p̂(1/conj(z))` or `p̂(z) = p̄(1/z) = (conj ∘ p ∘ conj ∘ inf)(z)`.

Examples:

```jldoctest laurent
julia> using Polynomials;

julia> z = variable(LaurentPolynomial, :z)
LaurentPolynomial(z)

julia> h = LaurentPolynomial([1,1], -1:0, :z)
LaurentPolynomial(z⁻¹ + 1)

julia> Polynomials.paraconj(h)(z) ≈ 1 + z ≈ LaurentPolynomial([1,1], 0:1, :z)
true

julia> h = LaurentPolynomial([3,2im,1], -2:0, :z)
LaurentPolynomial(3*z⁻² + 2im*z⁻¹ + 1)

julia> Polynomials.paraconj(h)(z) ≈ 1 - 2im*z + 3z^2 ≈ LaurentPolynomial([1, -2im, 3], 0:2, :z)
true

julia> Polynomials.paraconj(h)(z) ≈ (conj ∘ h ∘ conj ∘ inv)(z)
true
"""
function paraconj(p::LaurentPolynomial)
    cs = p.coeffs
    ds = adjoint.(cs)
    n = degree(p)
    LaurentPolynomial(reverse(ds), -n, indeterminate(p))
end

"""
    cconj(p)

Conjugation of a polynomial with respect to the imaginary axis.

The `cconj` of a polynomial, `p̃`, conjugates the coefficients and applies `s -> -s`. That is `cconj(p)(s) = conj(p)(-s)`.

This satisfies for *imaginary* `s`: `conj(p(s)) = p̃(s) = (conj ∘ p)(s) = cconj(p)(s) `

[ref](https://github.com/hurak/PolynomialEquations.jl#symmetrix-conjugate-equation-continuous-time-case)

Examples:
```jldoctest laurent
julia> using Polynomials;

julia> s = 2im
0 + 2im

julia> p = LaurentPolynomial([im,-1, -im, 1], 1, :s)
LaurentPolynomial(im*s - s² - im*s³ + s⁴)

julia> Polynomials.cconj(p)(s) ≈ conj(p(s))
true

julia> a = LaurentPolynomial([-0.12, -0.29, 1],:s)
LaurentPolynomial(-0.12 - 0.29*s + 1.0*s²)

julia> b = LaurentPolynomial([1.86, -0.34, -1.14, -0.21, 1.19, -1.12],:s)
LaurentPolynomial(1.86 - 0.34*s - 1.14*s² - 0.21*s³ + 1.19*s⁴ - 1.12*s⁵)

julia> x = LaurentPolynomial([-15.5, 50.0096551724139, 1.19], :s)
LaurentPolynomial(-15.5 + 50.0096551724139*s + 1.19*s²)

julia> Polynomials.cconj(a) * x + a * Polynomials.cconj(x) ≈ b + Polynomials.cconj(b)
true
```

"""
function cconj(p::LaurentPolynomial)
    ps = conj.(coeffs(p))
    m,n = (extrema ∘ degreerange)(p)
    for i in m:n
        if isodd(i)
            ps[i+1-m] *= -1
        end
    end
    LaurentPolynomial(ps, m, indeterminate(p))
end

##
## ----
##

# evaluation uses `evalpoly`
function evalpoly(x::S, p::LaurentPolynomial{T}) where {T,S}
    xᵐ = firstindex(p) < 0 ? inv(x)^(abs(firstindex(p))) : (x/1)^firstindex(p) # make type stable
    return EvalPoly.evalpoly(x, p.coeffs) * xᵐ
end

# scalar operations
# needed as standard-basis defn. assumes basis 1, x, x², ...
function Base.:+(p::LaurentPolynomial{T,X}, c::S) where {T, X, S <: Number}
    R = promote_type(T,S)
    q = LaurentPolynomial{R,X}(p.coeffs, firstindex(p))
    q[0] += c
    q
end

##
## Poly + and  *
##
function Base.:+(p1::P, p2::P) where {T,X,P<:LaurentPolynomial{T,X}}

    isconstant(p1) && return constantterm(p1) + p2
    isconstant(p2) && return p1 + constantterm(p2)

    m1,n1 = (extrema ∘ degreerange)(p1)
    m2,n2 = (extrema ∘ degreerange)(p2)
    m,n = min(m1,m2), max(n1, n2)

    as = zeros(T, length(m:n))
    for i in m:n
        as[1 + i-m] = p1[i] + p2[i]
    end

    q = P(as, m)
    chop!(q)

    return q

end

function Base.:*(p1::P, p2::P) where {T,X,P<:LaurentPolynomial{T,X}}

    m1,n1 = (extrema ∘ degreerange)(p1)
    m2,n2 = (extrema ∘ degreerange)(p2)
    m,n = m1 + m2, n1+n2

    as = zeros(T, length(m:n))
    for i in eachindex(p1)
        p1ᵢ = p1[i]
        for j in eachindex(p2)
            as[1 + i+j - m] = muladd(p1ᵢ, p2[j], as[1 + i + j - m])
        end
    end

    p = P(as, m)
    chop!(p)

    return p
end

##
## roots
##
"""
    roots(p)

Compute the roots of the Laurent polynomial `p`.

The roots of a function (Laurent polynomial in this case) `a(z)` are the values of `z` for which the function vanishes. A Laurent polynomial ``a(z) = a_m z^m + a_{m+1} z^{m+1} + ... + a_{-1} z^{-1} + a_0 + a_1 z + ... + a_{n-1} z^{n-1} + a_n z^n`` can equivalently be viewed as a rational function with a multiple singularity (pole) at the origin. The roots are then the roots of the numerator polynomial. For example, ``a(z) = 1/z + 2 + z`` can be written as ``a(z) = (1+2z+z^2) / z`` and the roots of `a` are the roots of ``1+2z+z^2``.

# Example

```julia
julia> using Polynomials;

julia> p = LaurentPolynomial([24,10,-15,0,1],-2:1,:z)
LaurentPolynomial(24*z⁻² + 10*z⁻¹ - 15 + z²)

julia> roots(a)
4-element Array{Float64,1}:
 -3.999999999999999
 -0.9999999999999994
  1.9999999999999998
  2.9999999999999982
```
"""
function  roots(p::P; kwargs...)  where  {T, X, P <: LaurentPolynomial{T, X}}
    c = coeffs(p)
    r = degreerange(p)
    d = r[end] - min(0, r[1]) + 1    # Length of the coefficient vector, taking into consideration
                                     # the case when the lower degree is strictly positive
                                     # (like p=3z^2).
    z = zeros(T, d)                  # Reserves space for the coefficient vector.
    z[max(0, r[1]) + 1:end] = c      # Leaves the coeffs of the lower powers as zeros.
    a = Polynomial{T,X}(z)           # The root is then the root of the numerator polynomial.
    return roots(a; kwargs...)
end

##
## d/dx, ∫
##
function derivative(p::P, order::Integer = 1) where {T, X, P<:LaurentPolynomial{T,X}}

    order < 0 && throw(ArgumentError("Order of derivative must be non-negative"))
    order == 0 && return p

    hasnan(p) && return ⟒(P)(T[NaN], 0, X)

    m,n = (extrema ∘ degreerange)(p)
    m = m - order
    n = n - order
    as =  zeros(T, length(m:n))

    for k in eachindex(p)
        idx = 1 + k - order - m
        if 0 ≤ k ≤ order - 1
            as[idx] = zero(T)
        else
            as[idx] = reduce(*, (k - order + 1):k, init = p[k])
        end
    end

    chop!(LaurentPolynomial{T,X}(as, m))

end

function integrate(p::P) where {T, X, P<: LaurentPolynomial{T, X}}

    !iszero(p[-1])  && throw(ArgumentError("Can't integrate Laurent  polynomial with  `x⁻¹` term"))
    R = eltype(one(T)/1)
    Q = ⟒(P){R, X}

    if hasnan(p)
        return Q([NaN],0)
    end

    m,n = (extrema ∘ degreerange)(p)
    if  n < 0
        n = 0
    else
        n += 1
    end
    if m < 0
        m += 1
    else
        m = 0
    end
    as = zeros(R,  length(m:n))

    for k in eachindex(p)
        as[1 + k+1-m]  =  p[k]/(k+1)
    end

    return Q(as, m)

end

function Base.gcd(p::LaurentPolynomial{T,X}, q::LaurentPolynomial{T,Y}, args...; kwargs...) where {T,X,Y}
    mp, Mp = (extrema ∘ degreerange)(p)
    mq, Mq = (extrema ∘ degreerange)(q)
    if mp < 0 || mq < 0
        throw(ArgumentError("GCD is not defined when there are `x⁻ⁿ` terms"))
    end

    degree(p) == 0 && return iszero(p) ? q : one(q)
    degree(q) == 0 && return iszero(q) ? p : one(p)
    assert_same_variable(p,q)

    pp, qq = convert(Polynomial, p), convert(Polynomial, q)
    u = gcd(pp, qq, args..., kwargs...)
    return LaurentPolynomial(coeffs(u), X)
end


text:
module TestBroadcasting

include("preamble.jl")

@testset "count_sum" begin
    @testset "$fold" for fold in [foldl, simple_reduce, foldxt_bs1, foldxt]
        foldl = nothing

        count_sum(xs, xf = Map(identity)) =
            fold(+, opcompose(xf, Map(x -> (1, x)), Broadcasting()), xs)
        @test count_sum([5, 2, 6, 8, 3]) == (5, 24)
        @test count_sum([5, 2, 6, 8, 3], Filter(isodd)) == (2, 8)
        @test_throws EmptyResultError count_sum(1:0)
        @test_throws EmptyResultError count_sum(1:10, Filter(_ -> false))
    end
end

@testset "partially started" begin
    @testset "$fold" for fold in [foldl, simple_reduce, foldxt_bs1, foldxt]
        foldl = nothing
        f(xs) =
            fold(max, opcompose(Map(x -> (x + 1, x)), Broadcasting(), Filter(isodd)), xs)
        @test f([5, 2, 6, 8, 3]) == (9, 5)
        @test_throws EmptyResultError f(2:2:8)
    end
end

@testset "init" begin
    @testset "$fold" for fold in [foldl, simple_reduce, foldxt_bs1, foldxt]
        foldl = nothing
        @test fold(+, Broadcasting(), 1:9) == 45
        @test fold(+, Broadcasting(), 1:9, init = [0]) == [45]
    end
end

sample_foldl_allocations(repeat::Integer, args...) =
    map(1:repeat) do _
        @allocated foldl(args...)
    end

@testset "allocation" begin
    args = (+, Broadcasting())
    n = 1000
    xs = [ones(n) for _ in 1:n]
    sample_foldl_allocations(2, args..., xs)
    baseline = minimum(sample_foldl_allocations(100, args..., xs))
    repeat10x = minimum(sample_foldl_allocations(100, args..., repeat(xs, 10)))
    @debug "estimated allocations" baseline repeat10x repeat10x / baseline

    # Always construct the messages to avoid introducing bugs.
    success_msg = """
    The test `repeat10x < 2 * baseline` works in Julia $(VERSION)!
    **Consider enabling this test again.**
    """
    skip_msg = """
    The test `repeat10x < 2 * baseline` does not work in Julia $(VERSION).
    Skipping this test.
    """
    if VERSION >= v"1.6-"
        if repeat10x < 2 * baseline
            @info success_msg baseline repeat10x repeat10x / baseline
        else
            @info skip_msg baseline repeat10x repeat10x / baseline
        end
        @test_skip repeat10x < 2 * baseline
    else
        @test repeat10x < 2 * baseline
    end
end

end  # module


text:
struct FilledContour{T}
    polylines::Polylines{T} # Polylines{T} defined in tricontour.jl
    lower::T
    upper::T
end
FilledContour(lower::T, upper::T) where T = FilledContour(Polylines{T}(), lower, upper)

"""
    tricontourf(x, y, z, t, levels)

Return filled contour lines of unstructured triangular data. `x`, `y`, and `z` are
one-dimensional input arrays of the same length (the number of points). `t` is an integer
array of size `(3,nt)`, where `nt` is the number of triangles. The coordinates of the `i`th
vertex of the `j`th triangle are given by `(x[t[i,j]], y[t[i,j]])`. Function values at
vertices are specified by the `z` array, and contours levels are defined by `levels.`

`tricontourf` results in regions *between* the contour levels, such that specifying `n`
levels will result in `n-1` regions. If `levels` is an array, it will be augmented to
include the minimum and maximum values of the `z` data.
"""
function tricontourf(x, y, z, t, levels)
    m = TriMesh(x, y, t)
    minz,maxz = extrema(z)
    if levels isa Integer
        levels = range(minz-eps(minz), maxz+eps(maxz), length=levels)
    else
        levels = copy(levels)
        minimum(levels) >= minz && pushfirst!(levels, minz-eps(minz))
        maximum(levels) <= maxz && push!(levels, maxz+eps(maxz))
    end
    tricontourf(m, z, levels)
end

function tricontourf(m::TriMesh, z, levels)
    @assert issorted(levels)
    nlevels = length(levels)
    lmin = minimum(levels)
    lmax = (nlevels > 1) ? maximum(levels[1:end-1]) : lmin
    filled_contours = FilledContour{eltype(levels)}[]
    for i=1:nlevels-1
        lower = levels[i]
        upper = levels[i+1]
        push!(filled_contours, generate_filled_contours(m, z, lower, upper))
    end
    filled_contours
end

function generate_filled_contours(m::TriMesh, z, lower, upper)
    reset_visited!(m)
    filled_contour = FilledContour(lower, upper)
    add_filled_bdr_lines!(filled_contour.polylines, m, z, lower, upper)
    add_interior_lines!(filled_contour.polylines, m, z, lower; on_upper=false, filled=true)
    add_interior_lines!(filled_contour.polylines, m, z, upper; on_upper=true,  filled=true)
    filled_contour
end

function follow_boundary!(polyline::Polyline, m::TriMesh, edge::TriEdge, z, lower, upper, on_upper)
    ib,ie = m.bdr_map[edge]
    m.bdr_used[ib] = true

    stop = false
    first_edge = true
    z_end = 0.0
    while !stop
        @assert !m.bdr_visited[ib][ie]
        m.bdr_visited[ib][ie] = true
        if first_edge
            z_start = z[m.t[edge.ie,edge.it]]
        else
            z_start = z_end
        end
        z_end = z[m.t[mod1(edge.ie+1,3),edge.it]]

        if z_end > z_start
            if (on_upper || !first_edge) && (z_end >= lower) && (z_start < lower)
                stop = true
                on_upper = false
            elseif (z_end >= upper) && (z_start < upper)
                stop = true
                on_upper = true
            end
        else
            if (!on_upper || !first_edge) && (z_start >= upper) && (z_end < upper)
                stop = true
                on_upper = true
            elseif (z_start >= lower) && (z_end < lower)
                stop = true
                on_upper = false
            end
        end
        first_edge = false
        if !stop
            ie = mod1(ie+1, length(m.bdr_loops[ib]))
            edge = m.bdr_loops[ib][ie]
            point = m.t[edge.ie,edge.it]
            push!(polyline, (m.x[point], m.y[point]))
        end
    end
    return on_upper,edge
end

function add_filled_bdr_lines!(polylines::Polylines{T}, m::TriMesh, z, lower, upper) where T
    for i=1:length(m.bdr_loops)
        bdr = m.bdr_loops[i]
        for j=1:length(bdr)
            m.bdr_visited[i][j] && continue
            edge = bdr[j]

            z_start = z[m.t[edge.ie,edge.it]]
            z_end = z[m.t[mod1(edge.ie+1,3),edge.it]]

            # Does this boundary edge's z increase through upper level and/or decrease
            # through lower level?
            incr_upper = (z_start < upper && z_end >= upper)
            decr_lower = (z_start >= lower && z_end < lower)

            if decr_lower || incr_upper
                start_edge = edge
                polyline = XY{T}[]
                on_upper = incr_upper
                while true
                    level_int = on_upper ? upper : lower
                    edge = follow_interior!(polyline, m, edge, z, level_int; end_on_bdr=true, on_upper)
                    on_upper,edge = follow_boundary!(polyline, m, edge, z, lower, upper, on_upper)
                    edge == start_edge && break
                end
                # Ensure closed loop
                push!(polyline, first(polyline))
                push!(polylines, polyline)
            end
        end
    end

    # Add full boundaries that lie between the lower and upper levels.  These are boundaries
    # that have not been touched by an internal contour line.
    for i=1:length(m.bdr_loops)
        m.bdr_used[i] && continue
        bdr = m.bdr_loops[i]
        zval = z[m.t[bdr[1].ie,bdr[1].it]]
        if (zval >= lower) && (zval < upper)
            polyline = XY{T}[]
            for j=1:length(bdr)
                edge = bdr[j]
                point = m.t[edge.ie,edge.it]
                push!(polyline, (m.x[point], m.y[point]))
            end
            push!(polyline, first(polyline))
            push!(polylines, polyline)
        end
    end
end


text:
#
""" use Base.ReshapedArray """
_reshape(a,dims::NTuple{D,Int}) where{D} = reshape(a,dims)
_reshape(a::Array, dims::NTuple{D,Int}) where{D} = Base.ReshapedArray(a, dims, ())

_vec(a) = vec(a)
_vec(a::AbstractVector) = a
_vec(a::AbstractArray) = _reshape(a,(length(a),))

""" check if operator(s) is square """
issquare(::UniformScaling) = true
issquare(A) = size(A,1) === size(A,2)
issquare(A...) = @. (&)(issquare(A)...)
#issquare(::Identity) = true
#


text:
include("vae_mnist.jl")

using Plots

function plot_result()
    BSON.@load "output/model.bson" encoder decoder args
    args = Args(; args...)
    device = args.cuda && CUDA.has_cuda() ? gpu : cpu
    encoder, decoder = encoder |> device, decoder |> device
    # load MNIST images
    loader = get_data(args.batch_size)

    # clustering in the latent space
    # visualize first two dims
    plt = scatter(palette=:rainbow)
    for (i, (x, y)) in enumerate(loader)
        i < 20 || break
        μ, logσ = encoder(x |> device)
        scatter!(μ[1, :], μ[2, :],
            markerstrokewidth=0, markeralpha=0.8,
            aspect_ratio=1,
            markercolor=y, label="")
    end
    savefig(plt, "output/clustering.png")

    z = range(-2.0, stop=2.0, length=11)
    len = Base.length(z)
    z1 = repeat(z, len)
    z2 = sort(z1)
    x = zeros(Float32, args.latent_dim, len^2) |> device
    x[1, :] = z1
    x[2, :] = z2
    samples = decoder(x)
    image = convert_to_image(samples, len)
    save("output/manifold.png", image)
end

if abspath(PROGRAM_FILE) == @__FILE__
    plot_result()
end


text:
using NonLinearReactionAdvectionDiffusionWithFrontData;
using NonLinearReactionAdvectionDiffusionWithFrontData: phidetermination, Φ;
using NonLinearReactionAdvectionDiffusionWithFrontData: f1, f2;
using NonLinearReactionAdvectionDiffusionWithFrontData: shishkin_mesh;
using NonLinearReactionAdvectionDiffusionWithFrontData: apply_on_dynamic_mesh;

using LinearAlgebra;

# Регрессивный тест
# На точных данных, градиент должен быть мал
@testset "Градиент на статической сетке, точные данные    " begin

    u_l(t) = -8
    u_r(t) =  4
    qf(x) = 4*sin(3 * π * x);        # Коэффициент линейного усиления, который в обратной
                                    # задаче необходимо определить, но при генерации априорной
                                    # информации мы задаем некоторый коэффициент, который,
                                    # собственно, после имея априорную информацию и будем определять.
    ε = 0.2;                        # Малый параметр при старшей производной
    a, b = 0, 1;                    # Область по X
    t₀, T = 0, 0.28;                # Область по T
    N, M = 50, 80;                  # Увеличенное Кол-во разбиений по X, T
    h = (b-a)/N;                    # шаг по X
    τ = (T-t₀)/M;                   # шаг по T
    Xₙ = [a  + n*h for n in 0:N];   # Сетка по Х
    Tₘ = [t₀ + m*τ for m in 0:M];   # Сетка по Т
    qₙ =      qf.(Xₙ);               # Сеточные значения коэффициента линейного усиления
    ulₘ=    u_l.(Tₘ);               # Сеточные значения левого  ГУ
    urₘ=    u_r.(Tₘ);               # Сеточные значения правого ГУ
    y₀ = u_init.(Xₙ);               # Начальные условия
    #
    u, XX, TP = solve(y₀, Xₙ, N, Tₘ, M, ε, ulₘ, urₘ, qₙ);
    nothing #hide

    #' ## Генерация априорной информации
    ϕl      = phidetermination(qₙ, ulₘ, Xₙ, N, Tₘ, M);                  # Левый вырожденный корень
    ϕr      = phidetermination(qₙ, urₘ, Xₙ, N, Tₘ, M, reverseX = true); # Правый вырожденный корень
    ϕ       = Φ(ϕl, ϕr, N, M);                                          # Серединный корень
    f1_data = f1(ϕ, u, Xₙ, N, M);                                       # Положение переходного слоя
    f2_data = f2(f1_data, u, Xₙ, N, M);                                 # Значение функции на переходном слое
    nothing #hide

    #' ## Решение сопряженной задачи

    Uₙₘ = u;                        # Сохраним старую матрицу
    y₀ = [0.0 for i in 1:N+1];      # Нулевые начальные условия
    ψl = [0.0 for i in 1:M+1];      # Нулевые ГУ
    ψr = [0.0 for i in 1:M+1];      # Нулевые ГУ

    ψ = solve_adjoint(y₀, Xₙ, N, Tₘ, M, ε, ψl, ψr, qₙ, Uₙₘ, f1_data, f2_data, w = 0.0005)

    grad = J_q(Uₙₘ, ψ, XX, N, Tₘ, M)

    # Значение для этих данных чисто эмпирическое
    # Мы просто проверяем, что градиент достаточно мал
    # "Достаточно" — понятие субъективное
    # Поэтому при тестировании, мы можем лишь удостовериться,
    # Что не произошло регрессии
    @test norm(grad) < 0.05             # L2 норма
    @test norm(grad, Inf) < 0.01        # Максимальный элемент в векторе max(abs.(v))

end

# Регрессивный тест
# На точных данных, градиент должен быть мал
@testset "Градиент на динамической сетке, точные данные   " begin

    u_l(t) = -8;                    # ГУ
    u_r(t) =  4;                    #
    qf(x) = 4*sin(3 * π * x);       # Коэффициент линейного усиления, который в обратной
                                    # задаче необходимо определить, но при генерации априорной
                                    # информации мы задаем некоторый коэффициент, который,
                                    # собственно, после имея априорную информацию и будем определять.
    ε = 0.03;                       # Малый параметр при старшей производной
    a, b = 0, 1;                    # Область по X
    t₀, T = 0, 0.20;                # Область по T
    x_tp = 0.22;                    # Положение переходного слоя
    M = 500;                        # Кол-во разбиений по X, T
    τ = (T-t₀)/M;                   # шаг по T
    Tₘ = [t₀ + m*τ for m in 0:M];   # Сетка по Т
    ulₘ=    u_l.(Tₘ);               # Сеточные значения левого  ГУ
    urₘ=    u_r.(Tₘ);               # Сеточные значения правого ГУ
    nothing #hide

    #' Новое поведение будет реализоваться с помощью передачи функции создания сетки внутрь `solve` в качестве аргумента.
    #' Эта функция должна принимать только один аргумент — положение переходного слоя и формировать соответствующую сетку.
    #' В пакете есть формирование кусочно-равномерной сетки со сгущениями на границе и на переходном слое.
    #' Создадим замыкание этой функции, которое будет иметь нужную сигнатуру.
    mshfrm(x_tp) = shishkin_mesh(a, b, x_tp, ε, 40, 1.0, 1.0, 0.2, 1.0);
    Xₙ = mshfrm(x_tp); ;                            # Сетка по Х
    N = length(Xₙ) - 1                              # Примем за N длину сетки, что получилась в итоге.
    qₙ =  1.0 .* qf.(Xₙ);                            # Сеточные значения коэффициента линейного усиления
                                                    # Измените линейный коэффициент, чтобы посмотреть
                                                    # на поведение сопряженной задачи (а следовательно градиента),
                                                    # когда прямая задача не соответствует наблюдаемым данным.
    u₀ = u_init.(Xₙ, ε=ε, x_tp = x_tp);             # Начальные условия

    u, XX, TP = solve(u₀, Xₙ, N, Tₘ, M, ε, ulₘ, urₘ, qₙ, create_mesh = mshfrm);
    #' ## Генерация априорной информации
    ϕl      = phidetermination(qₙ, ulₘ, Xₙ, N, Tₘ, M);                          # Левый вырожденный корень
    ϕr      = phidetermination(qₙ, urₘ, Xₙ, N, Tₘ, M, reverseX = true);         # Правый вырожденный корень
    ϕ       = Φ(ϕl, ϕr, N, M);                                                  # Полуразность вырожденных корней
    ϕ       = apply_on_dynamic_mesh(ϕ, XX, N, M);                               # Аппроксимация на переменную сетку
    ϕl      = apply_on_dynamic_mesh(ϕl, XX, N, M);                              # Аппроксимация на переменную сетку
    ϕr      = apply_on_dynamic_mesh(ϕr, XX, N, M);                              # Аппроксимация на переменную сетку
    f1_data = f1(ϕ, u, XX, N, M);                                               # Положение переходного слоя
    f2_data = f2(f1_data, u, XX, N, M);                                         # Значение функции на переходном слое
    #
    w = 0.000011;                                                               # Априорный параметр апроксимации дельта-функции

    #' ## Решение сопряженной задачи

    Uₙₘ = u;                        # Сохраним старую матрицу
    ψ₀ = [0.0 for i in 1:N+1];      # Нулевые начальные условия
    ψl = [0.0 for i in 1:M+1];      # Нулевые ГУ
    ψr = [0.0 for i in 1:M+1];      # Нулевые ГУ
    ψ = solve_adjoint(ψ₀, XX, N, Tₘ, M, ε, ψl, ψr, qₙ, Uₙₘ, f1_data, f2_data, w = w)

    # Ошибка, если не передадим вспомогательный массив
    @test_throws ArgumentError J_q(Uₙₘ, ψ, XX, N, Tₘ, M)

    X_q = [i/50.0 for i in 0:50]
    grad = J_q(Uₙₘ, ψ, XX, N, Tₘ, M, X_q)

    #=

    plot(Xₙ, grad)

    =#

    # Значение для этих данных чисто эмпирическое
    # Мы просто проверяем, что градиент достаточно мал
    # "Достаточно" — понятие субъективное
    # Поэтому при тестировании, мы можем лишь удостовериться,
    # Что не произошло регрессии
    #
    # Тесты помечены broken, потому что алгоритм неверен
    @test norm(grad) < 0.05             # L2 норма
    @test norm(grad, Inf) < 0.01        # Максимальный элемент в векторе max(abs.(v))

end


text:
using Distributions, ProgressMeter

function cross_over(survivors)
    size = length(survivors) * 2
    population = Array{Array{Float64}}(undef, size)
    for i = 1:size
        parents = rand(DiscreteUniform(1, length(survivors)), 2)
        genes = rand(DiscreteUniform(1, 2), 4)
        population[i] = [survivors[parents[genes[1]]][1],
                         survivors[parents[genes[2]]][2],
                         survivors[parents[genes[3]]][3],
                         survivors[parents[genes[4]]][4]]
    end
    return population
end

function mutate(population;
    mutation_probability = 0.25,
    qcc_σ = 0.5,
    σqcc_σ = 0.05,
    η_σ = 0.1,
    ση_σ = 0.05
)
    for i = 1:length(population)
        if rand(Uniform(0, 1)) < mutation_probability
            population[i][1] = max(0, rand(Normal(population[i][1], qcc_σ)))
        end
        if rand(Uniform(0, 1)) < mutation_probability
            population[i][2] = max(0, rand(Normal(population[i][2], σqcc_σ)))
        end
        if rand(Uniform(0, 1)) < mutation_probability
            population[i][3] = max(0, rand(Normal(population[i][3], η_σ)))
        end
        if rand(Uniform(0, 1)) < mutation_probability
            population[i][4] = max(0, rand(Normal(population[i][4], ση_σ)))
        end
    end
    return population
end

function genetic_algorithm(N, experimental; gens = 1_000)
    mins, maxes, means = zeros(gens), zeros(gens), zeros(gens)
    population = [get_quadrupolar_starting_values.(ones(Int, N)) zeros(N)]
    best = (population[1, 1], Inf)
    experimental_ecdf = get_experimental_ecdf(experimental)
    ν0 =  get_ν0(experimental, experimental_ecdf)
    @showprogress for i = 1:gens
        population[:, 2] = ols_cdf.(Quadrupolar.(population[:, 1]),
            Ref(experimental), Ref(experimental_ecdf), ν0)
        population = sortslices(population, by=x->x[2], dims = 1, rev = false)
        mins[i], maxes[i] = minimum(population[:, 2]), maximum(population[:, 2])
        if mins[i] < best[2]
            best = population[argmin(population[:, 2]), 1], mins[i]
        end
        means[i] = mean(population[:, 2])
        population[:, 1] = cross_over(population[1:Int(N / 2), 1])
        population[:, 1] = mutate(population[:, 1])
        map(x -> x[4] = max(eps(Float64), x[4]), population[:, 1])
    end
    return best, population, mins, means, maxes
end

function cross_over_csa(survivors)
    size = length(survivors) * 2
    population = Array{Array{Float64}}(undef, size)
    for i = 1:size
        parents = rand(DiscreteUniform(1, length(survivors)), 2)
        genes = rand(DiscreteUniform(1, 2), 6)
        population[i] = [survivors[parents[genes[1]]][1],
                         survivors[parents[genes[2]]][2],
                         survivors[parents[genes[3]]][3],
                         survivors[parents[genes[4]]][4],
                         survivors[parents[genes[3]]][5],
                         survivors[parents[genes[4]]][6],
                         1.0]
    end
    return population
end

function mutate_csa(population;
    mutation_probability = 0.25,
    σᵢₛₒ_σ = 50,
    σσᵢₛₒ_σ = 10,
    Δσ_σ = 50,
    σΔσ_σ = 10,
    ησ_σ = 0.1,
    σησ_σ = 0.05
)
    for i = 1:length(population)
        if rand(Uniform(0, 1)) < mutation_probability
            population[i][1] = max(0, rand(Normal(population[i][1], σᵢₛₒ_σ)))
        end
        if rand(Uniform(0, 1)) < mutation_probability
            population[i][2] = max(0, rand(Normal(population[i][2], σσᵢₛₒ_σ)))
        end
        if rand(Uniform(0, 1)) < mutation_probability
            population[i][3] = max(0, rand(Normal(population[i][3], Δσ_σ)))
        end
        if rand(Uniform(0, 1)) < mutation_probability
            population[i][4] = max(0, rand(Normal(population[i][4], σΔσ_σ)))
        end
        if rand(Uniform(0, 1)) < mutation_probability
            population[i][5] = max(0, rand(Normal(population[i][5], ησ_σ)))
        end
        if rand(Uniform(0, 1)) < mutation_probability
            population[i][6] = max(0, rand(Normal(population[i][6], σησ_σ)))
        end
    end
    return population
end

function genetic_algorithm_csa(N, experimental; gens = 1_000)
    mins, maxes, means = zeros(gens), zeros(gens), zeros(gens)
    population = [get_chemical_shift_starting_values.(ones(Int, N)) zeros(N)]
    best = (population[1, 1], Inf)
    experimental_ecdf = get_experimental_ecdf(experimental)
    @showprogress for i = 1:gens
        population[:, 2] = ols_cdf.(ChemicalShift.(population[:, 1]),
            Ref(experimental), Ref(experimental_ecdf))
        population = sortslices(population, by=x->x[2], dims = 1, rev = false)
        mins[i], maxes[i] = minimum(population[:, 2]), maximum(population[:, 2])
        if mins[i] < best[2]
            best = population[argmin(population[:, 2]), 1], mins[i]
        end
        means[i] = mean(population[:, 2])
        population[:, 1] = cross_over_csa(population[1:Int(N / 2), 1])
        population[:, 1] = mutate_csa(population[:, 1])
        map(x -> x[2] = max(eps(Float64), x[2]), population[:, 1])
        map(x -> x[4] = max(eps(Float64), x[4]), population[:, 1])
        map(x -> x[5] = max(eps(Float64), x[5]), population[:, 1])
        map(x -> x[6] = max(eps(Float64), x[6]), population[:, 1])
        map(x -> x[5] = min(1.0, x[5]), population[:, 1])
    end
    return best, population, mins, means, maxes
end


text:
include("../required.jl")

transect

##Create Adjacency matrix
#create matrix of zeros with n = colPred, m = colPrey
M = zeros(Int16, length(unique(transect.newPredID))+1, length(unique(transect.newPreyID))+1)
#assign names for cols
for y in 2:length(unique(transect.newPredID))+1
    M[y,1] = unique(transect.newPredID)[y-1]
end

#assign names for rows
for x in 2:length(unique(transect.newPreyID))+1
    M[1,x] = unique(sort!(transect.newPreyID))[x-1]
end

# Fill adjacency matrix with ones where there are interactions
u_pred_len = length(unique(transect.newPredID))
pred_len = length(transect.newPredID)
for i in 1:pred_len
    M[transect[i, 1] + 1, transect[i, 2] - u_pred_len  + 1] = 1
end
M

# # Test to make sure the adjacency matrix is correctly filled
# transectTest = [[1,6],[2,7], [3,8], [4,9],[5, 10]]
# MTest = zeros(Int16, 6, 6)
# for y in 2:6
#     MTest[y,1] = transectTest[y-1][1]
# end
#
# for x in 2:6
#     MTest[1,x] = transectTest[x-1][2]
# end
#
# for i in 1:5
#     MTest[i + 1, transectTest[i][2] - 5 + 1] = 1
# end

max_stars_repo_name: kiristern/reseaux

max_issues_repo_name: kiristern/reseaux

max_forks_repo_name: kiristern/reseaux


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule htslib_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("htslib")
JLLWrappers.@generate_main_file("htslib", UUID("f06fe41e-9474-5571-8c61-5634d2b2700c"))
end  # module htslib_jll


text:
#for now just the main functions, no wrapper
function optimize_pwr_φ_common!(φs, last_φs, opb, power_buffer, ids, scatteringMatrices, P, Ns, fmmbuf, fmmopts, mFMM)
    if φs != last_φs
        copyto!(last_φs, φs)
        #do whatever common calculations and save to shared_var
        for i in eachindex(opb)
            for ic = 1:Ns
                rng = (ic-1)*(2*P+1) .+ (1:2*P+1)
                if φs[ic] == 0.0
                    fmmbuf.rhs[rng] = scatteringMatrices[i][ids[ic]]*opb[i].α[rng]
                else
                    #rotate without matrix
                    rotateMultipole!(view(fmmbuf.rhs,rng), view(opb[i].α,rng), -φs[ic], P)
                    fmmbuf.rhs[rng] = scatteringMatrices[i][ids[ic]]*fmmbuf.rhs[rng]
                    rotateMultipole!(view(fmmbuf.rhs,rng), φs[ic], P)
                end
            end

            MVP = LinearMap{eltype(fmmbuf.rhs)}(
                        (output_, x_) -> ParticleScattering.FMM_mainMVP_pre!(output_, x_,
                                            scatteringMatrices[i], φs, ids, P, mFMM[i],
                                            fmmbuf.pre_agg, fmmbuf.trans),
                        Ns*(2*P+1), Ns*(2*P+1), ismutating = true)

            opb[i].β[:] .= 0
            opb[i].β,ch = gmres!(opb[i].β, MVP, fmmbuf.rhs,
                                    restart = Ns*(2*P+1), tol = fmmopts.tol,
                                    log = true, initially_zero = true) #no restart
            !ch.isconverged && error("""FMM process did not converge, normalized
                residual: $(norm(MVP*opb[i].β - fmmbuf.rhs)/norm(fmmbuf.rhs))""")
        end
        #calculate power for each power calculation arc
        calc_multi_pwr!(power_buffer, opb)
    end
end

function optimize_pwr_φ_f(φs, last_φs, opb, power_buffer, ids, scatteringMatrices, P, Ns, fmmbuf, fmmopts, mFMM, fobj_pwr)
    optimize_pwr_φ_common!(φs, last_φs, opb, power_buffer, ids, scatteringMatrices, P, Ns, fmmbuf, fmmopts, mFMM)

    fobj_pwr(power_buffer)
end

function optimize_pwr_φ_g!(grad_stor, φs, last_φs, opb, power_buffer, ids, scatteringMatrices, scatteringLU, P, Ns, fmmbuf, fmmopts, mFMM, gobj_pwr!)
    optimize_pwr_φ_common!(φs, last_φs, opb, power_buffer, ids, scatteringMatrices, P, Ns, fmmbuf, fmmopts, mFMM)

    fill!(grad_stor, 0) #gradient is the sum of all adjoint gradients
    #compute all ∂P/∂βᵀ
    dPdβ_pwr!.(power_buffer)
    #build rhs of all adjoint problems (-∂f/∂βᵀ)
    gobj_pwr!(power_buffer, opb)

    #for each problem, solve adjoint problem and add to total gradient
    for i in eachindex(opb)
        MVP = LinearMap{eltype(fmmbuf.rhs)}(
                (output_, x_) -> ParticleScattering.FMM_mainMVP_transpose!(output_, x_,
                                        scatteringMatrices[i], φs, ids, P, mFMM[i],
                                        fmmbuf.pre_agg, fmmbuf.trans),
                Ns*(2*P+1), Ns*(2*P+1), ismutating = true)

        #solve adjoint problem
        opb[i].λadj[:] .= 0
        opb[i].λadj, ch = gmres!(opb[i].λadj, MVP, opb[i].rhs_grad,
                            restart = Ns*(2*P+1), tol = fmmopts.tol,
                            log=true, initially_zero = true) #no restart
        !ch.isconverged && error("""FMM process did not converge for adjoint,
            normalized residual:
            $(norm(MVP*opb[i].λadj - opb[i].rhs_grad)/norm(opb[i].rhs_grad))""")
        D = -1.0im*collect(-P:P)
        v = Array{Complex{Float64}}(undef, 2*P+1) #TODO: minimize dynamic alloc
        v2 = Array{Complex{Float64}}(undef, 2*P+1)
        for n = 1:Ns
            rng = (n-1)*(2*P+1) .+ (1:2*P+1)
            rotateMultipole!(v, view(opb[i].β,rng), -φs[n], P)
            ldiv!(scatteringLU[i][ids[n]], v)
            v2[:] = scatteringMatrices[i][ids[n]]*v
            v2 .*= D
            v .*= D
            v[:] = scatteringMatrices[i][ids[n]]*v
            v2 .-= v
            rotateMultipole!(v2, φs[n], P)
            grad_stor[n] += (-2)*real(transpose(opb[i].λadj[rng])*v2)
        end
    end
end


text:
"""
A [`Message`](@ref) attachment.
More details [here](https://discordapp.com/developers/docs/resources/channel#attachment-object).
"""
struct Attachment
    id::Snowflake
    filename::String
    size::Int
    url::String
    proxy_url::String
    height::Optional{Int}
    width::Optional{Int}
end
@boilerplate Attachment :constructors :docs :lower :merge :mock


text:
@testset "Controller: MPC flamingo" begin
    # ## Simulation
    s = get_simulation("flamingo", "flat_2D_lc", "flat")
    model = s.model
    env = s.env

    # ## Reference Trajectory
    ref_traj = deepcopy(get_trajectory(s.model, s.env,
        joinpath(module_dir(), "src/dynamics/flamingo/gaits/gait_forward_36_4.jld2"),
        load_type = :split_traj_alt));

    H = ref_traj.H
    h = ref_traj.h

    # ## MPC setup
    N_sample = 5
    H_mpc = 15
    h_sim = h / N_sample
    H_sim = 1000
    κ_mpc = 2.0e-4

    obj = TrackingVelocityObjective(model, env, H_mpc,
        v = [Diagonal(1e-3 * [1e0,1,1e4,1,1,1,1,1e4,1e4]) for t = 1:H_mpc],
        q = [Diagonal(1e-1 * [3e2, 1e-6, 3e2, 1, 1, 1, 1, 0.1, 0.1]) for t = 1:H_mpc],
        u = [Diagonal(3e-1 * [0.1; 0.1; 0.3; 0.3; ones(model.nu-6); 2; 2]) for t = 1:H_mpc],
        γ = [Diagonal(1.0e-100 * ones(model.nc)) for t = 1:H_mpc],
        b = [Diagonal(1.0e-100 * ones(model.nc * friction_dim(env))) for t = 1:H_mpc]);

    p = ci_mpc_policy(deepcopy(ref_traj), s, obj,
        H_mpc = H_mpc,
        N_sample = N_sample,
        κ_mpc = κ_mpc,
        ip_opts = InteriorPointOptions(
                        undercut = 5.0,
                        κ_tol = κ_mpc,
                        r_tol = 1.0e-8,
                        diff_sol = true,
                        solver=:empty_solver,
                        max_time = 1e5),
        n_opts = NewtonOptions(
            r_tol = 3e-4,
            max_iter = 5),
        mpc_opts = CIMPCOptions(
            altitude_update=false,
            altitude_impact_threshold=0.02
        ));

    # ## Initial conditions
    q1_sim, v1_sim = initial_conditions(deepcopy(ref_traj));

    # ## Simulator
    sim = simulator(s, H_sim, h=h_sim, policy=p,
        solver_opts=InteriorPointOptions(
            undercut=Inf,
            γ_reg=0.0,
            r_tol=1e-8,
            κ_tol=1e-8,
            ϵ_min=0.05,
            diff_sol=false,
            verbose=false));

    # ## Simulate
    status = simulate!(sim, q1_sim, v1_sim);
    @test status

    @test norm(sim.traj.q[1] - ref_traj.q[2]*(1-1/N_sample) - ref_traj.q[1]/N_sample) < 1e-8
    @test norm(sim.traj.q[2] - ref_traj.q[2]) < 1e-8
    @test norm(ref_traj.q[1][2:end] - ref_traj.q[end-1][2:end]) < 1e-8
    @test norm(ref_traj.q[2][2:end] - ref_traj.q[end-0][2:end]) < 1e-8

    qerr, uerr, γerr, berr = ContactImplicitMPC.tracking_error(deepcopy(ref_traj), sim.traj, N_sample, idx_shift=[1])
    @test qerr < 0.0154 * 1.5 # 0.0154
    @test uerr < 0.0829 * 1.5 # 0.0829
    @test γerr < 0.444 * 1.5 # 0.444
    @test berr < 0.0169 * 1.5 # 0.0169
    qerr > 0.0154 * 1.2 && @warn "mild regression on q tracking: current tracking error = $qerr, nominal tracking error = 0.0154"
    uerr > 0.0829 * 1.2 && @warn "mild regression on u tracking: current tracking error = $uerr, nominal tracking error = 0.0829"
    γerr > 0.4440 * 1.2 && @warn "mild regression on γ tracking: current tracking error = $γerr, nominal tracking error = 0.444"
    berr > 0.0169 * 1.2 && @warn "mild regression on b tracking: current tracking error = $berr, nominal tracking error = 0.0169"
end


text:
## set up functions for multiplet calculations

"Return a matrix of simulated intensities of size (δH x δC x τ)"
function multiplet(param, δobsH, δobsC, τ, bfH, bfC, windowH, windowC)
  δH = param[1]
  δC = param[2]
  R2H = param[3]
  R2C = param[4]
  S2tc = param[5]
  csa = param[6]*1e-6
  J = param[7]

  # convert chemical shift axes to angular frequencies
  ωobsH = reshape(2π * δobsH * bfH, :, 1, 1)
  ωobsC = reshape(2π * δobsC * bfC, 1, :, 1)
  τobs = reshape(τ, 1, 1, :)

  # convert peak positions to angular frequencies
  ωH = 2π * δH * bfH
  ωC = 2π * δC * bfC

  # physical parameters
  μ₀ = 4π * 1e-7
  ħ = 1.055e-34
  γH = 2.675e8
  γC = 6.726e7
  rCH = 1.117e-10  # Tugarinov 2004
  rHH = sqrt(3) * rCH * sin(110.4*π/180)  # Tugarinov 2004
  ω0C = 2π * bfC * 1e6
  P2cosβ = -1/3.  # 109 degree tetrahedral geometry
  TAU = 1.0 / (2*125)  # coherence transfer period

  # relaxation rates (as function of S2tc)
  σ = 2/45 * (μ₀*ħ*γC*γH / (4π*rCH^3))^2 * S2tc * 1e-9
  η = 2/5 * (μ₀/(4π)) * P2cosβ/rCH^3 * ħ * γH * γC * ω0C * csa * S2tc * 1e-9
  ηHHHC = 1/5 * (μ₀/(4π))^2 * rCH^-3 * rHH^-3 * ħ^2 * γH^3 * γC * S2tc * 1e-9
  Δ = exp(-TAU * (9/20 * (μ₀*ħ*γH^2 / (4π*rHH^3))^2 * S2tc * 1e-9)) * cosh(-TAU * ηHHHC)
  ηHHHH = 9/40 * (μ₀/(4π))^2 * rHH^-6 * ħ^2 * γH^4 * S2tc * 1e-9

  # initial intensities (observing slow relaxing 1H coherences only)
  Iouter = 3 + 3Δ
  Iinner = 3 - Δ

  # relaxation rates for 13C transitions
  Rααα = R2C + 3σ + 2η
  Rααβ = R2C - σ + 2/3 * η
  Rαββ = R2C - σ - 2/3 * η
  Rβββ = R2C + 3σ - 2η

  return @. (
        Iouter * exp(-Rααα * τobs) * lineshape(ωC + 3π*J, Rααα, ωobsC, windowC) +
        Iinner * exp(-Rααβ * τobs) * lineshape(ωC + π*J, Rααβ, ωobsC, windowC) +
        Iinner * exp(-Rαββ * τobs) * lineshape(ωC - π*J, Rαββ, ωobsC, windowC) +
        Iouter * exp(-Rβββ * τobs) * lineshape(ωC - 3π*J, Rβββ, ωobsC, windowC)
        ) * lineshape(ωH, R2H, ωobsH, windowH)
end

## window functions
function lineshape(ω, R, ωobs, w::ExponentialWindow)
  x = @. -R + 1im*(ωobs - ω) - π*w.lb
  T = w.tmax
  return @. real((1 - exp(T*x)) / x)
end

function lineshape(ω, R, ωobs, w::CosWindow)
  x = @. -R + 1im*(ωobs - ω)
  T = w.tmax
  Tx = T * x
  return @. real(T * (π*exp(Tx) + 2*Tx) / (π^2 + 4*Tx^2))
end

function lineshape(ω, R, ωobs, w::Cos²Window)
  x = @. -R + 1im*(ωobs - ω)
  Tx = w.tmax * x
  return @. real((π^2*(1-exp(Tx)) + 2*Tx^2) / ((π^2 + Tx^2) * x))
end

Base.Broadcast.broadcastable(w::WindowFunction) = Ref(w)


text:
# =======================
# x' = Ax + U
# =======================

# Set representation: Interval
# Matrix operations: Dense
# Invariant: No
function reach_inhomog_BFFPSV18!(F, Xhat0::LazySet{N}, Φ::MT, NSTEPS, δ, X::Universe, U,
                                 ST::Type{<:Interval{N}}, vars,
                                 block_indices,
                                 row_blocks::AbstractVector{<:RBLKi},
                                 column_blocks::AbstractVector{<:CBLKj},
                                 Δt0,
                                 viewval::Val{true}) where {NM, MT<:AbstractMatrix{NM}, N, RBLKi, CBLKj}

    # initial reach set
    Δt = (zero(N) .. δ) + Δt0
    @inbounds F[1] = SparseReachSet(CartesianProductArray(Xhat0.array[block_indices]), Δt, vars)

    # cache matrix
    Φpowerk = copy(Φ)
    Φpowerk_cache = similar(Φ)

    # buffer for overapproximated Minkowski sum of each row-block
    Xhatk = Vector{ST}(undef, length(row_blocks))
    zero_int = Interval(zero(N), zero(N))

    # preallocate accumulated inputs and decompose it
    Whatk = Vector{ST}(undef, length(row_blocks))
    @inbounds for (i, bi) in enumerate(row_blocks)
        πX = Projection(U, bi)
        Whatk[i] = overapproximate(πX, ST)
    end

    @inbounds for k in 2:NSTEPS
        for (i, bi) in enumerate(row_blocks) # loop over row-blocks of interest
            Xhatk[i] = zero_int
            for (j, bj) in enumerate(column_blocks) # loop over all column-blocks
                # concrete interval-interval Minkowski sum
                aux = Φpowerk[bi[1], bj[1]] * Xhat0.array[j].dat
                Xhatk[i] = Interval(Xhatk[i].dat + aux)
            end
            Xhatk[i] = Interval(Xhatk[i].dat + Whatk[i].dat)
        end
        Δt += δ
        Xk = CartesianProductArray(copy(Xhatk))
        F[k] = SparseReachSet(Xk, Δt, vars)

        # update the input set
        for (i, bi) in enumerate(row_blocks)
            Whatk[i] = overapproximate(Whatk[i] ⊕ view(Φpowerk, bi, :) * U, ST)
        end

        mul!(Φpowerk_cache, Φpowerk, Φ)
        copyto!(Φpowerk, Φpowerk_cache)
    end
    return F
end

# Set representation: Generic
# Matrix operations: Dense
# Invariant: No
function reach_inhomog_BFFPSV18!(F, Xhat0::LazySet{N}, Φ::MT, NSTEPS, δ, X::Universe, U,
                                 ST, vars, block_indices,
                                 row_blocks::AbstractVector{<:RBLKi},
                                 column_blocks::AbstractVector{<:CBLKj},
                                 Δt0,
                                 viewval::Val{true}) where {NM, MT<:AbstractMatrix{NM}, N, RBLKi, CBLKj}

    # initial reach-set
    Δt = (zero(N) .. δ) + Δt0
    @inbounds F[1] = SparseReachSet(CartesianProductArray(Xhat0.array[block_indices]), Δt, vars)

    # cache matrix
    Φpowerk = copy(Φ)
    Φpowerk_cache = similar(Φ)

    # preallocate buffer for each row-block
    SMT = SubArray{N, 2, MT, Tuple{RBLKi, CBLKj}, false}
    buffer = Vector{LazySets.LinearMap{N, ST, N, SMT}}(undef, length(column_blocks))

    # preallocate overapproximated Minkowski sum for each row-block
    Xhatk = Vector{ST}(undef, length(row_blocks))

    # preallocate accumulated inputs and decompose it
    Whatk = Vector{ST}(undef, length(row_blocks))
    @inbounds for (i, bi) in enumerate(row_blocks)
        πX = Projection(U, bi)
        Whatk[i] = overapproximate(πX, ST)
    end

    @inbounds for k in 2:NSTEPS
        for (i, bi) in enumerate(row_blocks) # loop over row-blocks of interest
            for (j, bj) in enumerate(column_blocks) # loop over all column-blocks
                buffer[j] = view(Φpowerk, bi, bj) * Xhat0.array[j]
            end
            Xhatk[i] = overapproximate(MinkowskiSumArray(buffer) ⊕ Whatk[i], ST)
        end
        Δt += δ
        Xk = CartesianProductArray(copy(Xhatk))
        F[k] = SparseReachSet(Xk, Δt, vars)

        # update the input set
        for (i, bi) in enumerate(row_blocks)
            Whatk[i] = overapproximate(Whatk[i] ⊕ view(Φpowerk, bi, :) * U, ST)
        end

        mul!(Φpowerk_cache, Φpowerk, Φ)
        copyto!(Φpowerk, Φpowerk_cache)
    end
    return F
end

# Set representation: Generic
# Matrix operations: Sparse
# Invariant: No
function reach_inhomog_BFFPSV18!(F, Xhat0, Φ::MT, NSTEPS, δ::N,
                                 X::Universe, U,
                                 ST, vars,
                                 block_indices,
                                 row_blocks::AbstractVector{<:RBLKi},
                                 column_blocks::AbstractVector{<:CBLKj},
                                 Δt0,
                                 viewval::Val{true}) where {NM, IM, MT<:SparseMatrixCSC{NM, IM},
                                                            N, RBLKi, CBLKj}
    # store first element
    Δt = (zero(N) .. δ) + Δt0
    @inbounds F[1] = SparseReachSet(CartesianProductArray(Xhat0.array[block_indices]), Δt, vars)

    # cache matrix
    Φpowerk = copy(Φ)
    Φpowerk_cache = similar(Φ)

    # preallocate buffer for each row-block
    SMT = SubArray{N, 2, MT, Tuple{RBLKi, CBLKj}, false}
    buffer = Vector{LazySets.LinearMap{N, ST, N, SMT}}(undef, length(column_blocks))

    # preallocate overapproximated Minkowski sum for each row-block
    Xhatk = Vector{ST}(undef, length(row_blocks))

    # preallocate accumulated inputs and decompose it
    Whatk = Vector{ST}(undef, length(row_blocks))
    @inbounds for (i, bi) in enumerate(row_blocks)
        πX = Projection(U, bi)
        Whatk[i] = overapproximate(πX, ST)
    end

    @inbounds for k in 2:NSTEPS
        for (i, bi) in enumerate(row_blocks) # loop over row-blocks of interest
            for (j, bj) in enumerate(column_blocks) # loop over all column-blocks
                Φpowerk_bi_bj = Φpowerk
                if !iszero(Φpowerk_bi_bj)
                    buffer[j] = view(Φpowerk, bi, bj) * Xhat0.array[j]
                end
            end
            Xhatk[i] = overapproximate(MinkowskiSumArray(buffer) ⊕ Whatk[i], ST)
        end
        Δt += δ
        Xk = CartesianProductArray(copy(Xhatk))
        F[k] = SparseReachSet(Xk, Δt, vars)

        # update the input set
        for (i, bi) in enumerate(row_blocks)
            Whatk[i] = overapproximate(Whatk[i] ⊕ view(Φpowerk, bi, :) * U, ST)
        end

        mul!(Φpowerk_cache, Φpowerk, Φ)
        copyto!(Φpowerk, Φpowerk_cache)
    end
    return F
end

function reach_inhomog_BFFPSV18!(F, Xhat0, Φ::MT, NSTEPS, δ::N,
                                 X::Universe, U,
                                 ST, vars,
                                 block_indices,
                                 row_blocks::AbstractVector{<:RBLKi},
                                 column_blocks::AbstractVector{<:CBLKj},
                                 Δt0,
                                 viewval::Val{false}) where {NM, IM, MT<:SparseMatrixCSC{NM, IM}, N, RBLKi, CBLKj}

    # store first element
    Δt = (zero(N) .. δ) + Δt0
    @inbounds F[1] = SparseReachSet(CartesianProductArray(Xhat0.array[block_indices]), Δt, vars)

    # cache matrix
    Φpowerk = copy(Φ)
    Φpowerk_cache = similar(Φ)

    # preallocate buffer for each row-block
    SMT = SparseMatrixCSC{NM, IM}
    buffer = Vector{LazySets.LinearMap{N, ST, NM, SMT}}(undef, length(column_blocks))

    # preallocate overapproximated Minkowski sum for each row-block
    Xhatk = Vector{ST}(undef, length(row_blocks))

    # preallocate accumulated inputs and decompose it
    Whatk = Vector{ST}(undef, length(row_blocks))
    @inbounds for (i, bi) in enumerate(row_blocks)
        πX = Projection(U, bi)
        Whatk[i] = overapproximate(πX, ST)
    end

    @inbounds for k in 2:NSTEPS
        for (i, bi) in enumerate(row_blocks) # loop over row-blocks of interest
            for (j, bj) in enumerate(column_blocks) # loop over all column-blocks
                Φpowerk_bi_bj = Φpowerk
                if !iszero(Φpowerk_bi_bj)
                    buffer[j] = Φpowerk[bi, bj] * Xhat0.array[j]
                end
            end

            Xhatk[i] = overapproximate(MinkowskiSumArray(buffer) ⊕ Whatk[i], ST)
        end
        Δt += δ
        Xk = CartesianProductArray(copy(Xhatk))
        F[k] = SparseReachSet(Xk, Δt, vars)

        # update the input set
        for (i, bi) in enumerate(row_blocks)
            Whatk[i] = overapproximate(Whatk[i] ⊕ Φpowerk[bi, :] * U, ST)
        end

        mul!(Φpowerk_cache, Φpowerk, Φ)
        copyto!(Φpowerk, Φpowerk_cache)
    end
    return F
end

#######################################################################

# Set representation: Generic
# Matrix operations: Dense
# Invariant: LazySet
function reach_inhomog_BFFPSV18!(F, Xhat0::LazySet{N}, Φ::MT, NSTEPS, δ, X::LazySet, U,
                                 ST, vars, block_indices,
                                 row_blocks::AbstractVector{<:RBLKi},
                                 column_blocks::AbstractVector{<:CBLKj},
                                 Δt0,
                                 viewval::Val{true}) where {NM, MT<:AbstractMatrix{NM}, N, RBLKi, CBLKj}

    # initial reach-set
    Δt = (zero(N) .. δ) + Δt0
    @inbounds F[1] = SparseReachSet(CartesianProductArray(Xhat0.array[block_indices]), Δt, vars)

    # cache matrix
    Φpowerk = copy(Φ)
    Φpowerk_cache = similar(Φ)

    # preallocate buffer for each row-block
    SMT = SubArray{N, 2, MT, Tuple{RBLKi, CBLKj}, false}
    buffer = Vector{LazySets.LinearMap{N, ST, N, SMT}}(undef, length(column_blocks))

    # preallocate overapproximated Minkowski sum for each row-block
    Xhatk = Vector{ST}(undef, length(row_blocks))

    # preallocate accumulated inputs and decompose it
    Whatk = Vector{ST}(undef, length(row_blocks))
    @inbounds for (i, bi) in enumerate(row_blocks)
        πX = Projection(U, bi)
        Whatk[i] = overapproximate(πX, ST)
    end

    k = 2
    @inbounds while k <= NSTEPS
        for (i, bi) in enumerate(row_blocks) # loop over row-blocks of interest
            for (j, bj) in enumerate(column_blocks) # loop over all column-blocks
                buffer[j] = view(Φpowerk, bi, bj) * Xhat0.array[j]
            end
            Xhatk[i] = overapproximate(MinkowskiSumArray(buffer) ⊕ Whatk[i], ST)
        end
        Xk = CartesianProductArray(copy(Xhatk))
        _is_intersection_empty(X, Xk) && break
        Δt += δ
        F[k] = SparseReachSet(Xk, Δt, vars)

        # update the input set
        for (i, bi) in enumerate(row_blocks)
            Whatk[i] = overapproximate(Whatk[i] ⊕ view(Φpowerk, bi, :) * U, ST)
        end

        mul!(Φpowerk_cache, Φpowerk, Φ)
        copyto!(Φpowerk, Φpowerk_cache)

        k += 1
    end
    if k < NSTEPS + 1
        resize!(F, k-1)
    end
    return F
end


text:
abstract type AbstractWignerRotation end
struct TriavialWignerRotation <: AbstractWignerRotation
    k::Int
end
struct WignerRotation{N} <: AbstractWignerRotation
    k::Int
    ispositive::Bool
    iseven::Bool
end

const Arg0WignerRotation = WignerRotation{0}
const Arg2WignerRotation = WignerRotation{2}
const Arg3WignerRotation = WignerRotation{3}

WignerRotation{N}(k::Int, ispositive::Bool) where N =
    WignerRotation{N}(k, ispositive, true)

ispositive(wr::TriavialWignerRotation) = true
ispositive(wr::WignerRotation) = wr.ispositive

import Base:iseven
iseven(wr::TriavialWignerRotation) = true
iseven(wr::WignerRotation{0}) = true
iseven(wr::WignerRotation{3}) = true
iseven(wr::WignerRotation{2}) = wr.iseven

ijk(k::Int) = (k+1,k+2,k) |> x->mod.(x,Ref(Base.OneTo(3)))
ijk(wr::AbstractWignerRotation) = ijk(wr.k)
#
ij_from_k(k::Int) = ijk(k)

issequential(i,j) = (j-i) ∈ (1,-2)

"""
    wr(system_a, reference_b, particle_c=0)

Create a WignerRotation object of the right type based on provided indices.
The daughter particles are numbered 1,2,3, the mother particle is 0.
 - `system_a` tells which isobar is considered, and
 - `reference_b` tell which system is used as a reference.
For `system_a` and `reference_b` the spectator notations are used, i.e.
1 for the system (2,3), 2 for the system (3,1), and 3 for the system (1,2).
"""
function wr(system_a, reference_b, particle_c=0)
    system_a == reference_b && return TriavialWignerRotation(particle_c)
    S = issequential(system_a, reference_b)
    A,B = S ? (system_a, reference_b) : (reference_b, system_a)
    #
    particle_c == 0 && return Arg0WignerRotation(A,S)
    #
    particle_c ∉ (system_a, reference_b) && return Arg3WignerRotation(particle_c,S)
    #
    T = (particle_c == A)
    return Arg2WignerRotation(particle_c,!(S),T)
end

cosζ(wr::TriavialWignerRotation,σs,msq) = 1.0

function cosζ(wr::Arg0WignerRotation,σs,msq)
    i,j,k = ijk(wr)
    #
    s = msq[4]
    EE4s = (s+msq[i]-σs[i])*(s+msq[k]-σs[k])
    pp4s = sqrt(Kallen(s,msq[i],σs[i])*Kallen(s,msq[k],σs[k]))
    rest = σs[j]-msq[i]-msq[k]
    return (EE4s-2s*rest)/pp4s
end

function cosζ(wr::Arg2WignerRotation,σs,msq)
    i,j,k = ijk(wr)
    #
    if !(iseven(wr))
        i,j = j,i
    end
    #
    msq[k] ≈ 0 && return 1.0
    #
    s = msq[4]
    EE4mksq = (s+msq[k]-σs[k])*(σs[i]-msq[k]-msq[j])
    pp4mksq = sqrt(Kallen(s,msq[k],σs[k])*Kallen(msq[k],msq[j],σs[i]))
    rest = σs[j]-s-msq[j]
    return (2msq[k]*rest+EE4mksq)/pp4mksq
end

#
function cosζ(wr::Arg3WignerRotation,σs,msq)
    i,j,k = ijk(wr)
    #
    msq[k] ≈ 0 && return 1.0
    #
    s = msq[4]
    EE4m1sq = (σs[i]-msq[j]-msq[k])*(σs[j]-msq[k]-msq[i])
    pp4m1sq = sqrt(Kallen(σs[i],msq[j],msq[k])*Kallen(σs[j],msq[k],msq[i]))
    rest = msq[i]+msq[j]-σs[k]
    return (2msq[k]*rest+EE4m1sq)/pp4m1sq
end

# explicit
cosζ21_for1(σs, ms²) = cosζ(wr(2,1,1), σs, ms²)
cosζ21_for2(σs, ms²) = cosζ(wr(2,1,2), σs, ms²)
cosζ13_for1(σs, ms²) = cosζ(wr(1,3,1), σs, ms²)
cosζ13_for3(σs, ms²) = cosζ(wr(1,3,3), σs, ms²)
cosζ32_for3(σs, ms²) = cosζ(wr(3,2,3), σs, ms²)
cosζ32_for2(σs, ms²) = cosζ(wr(3,2,2), σs, ms²)

cosζ12_for3(σs, ms²) = cosζ(wr(1,2,3), σs, ms²)
cosζ23_for1(σs, ms²) = cosζ(wr(2,3,1), σs, ms²)
cosζ31_for2(σs, ms²) = cosζ(wr(3,1,2), σs, ms²)

cosζ12_for0(σs, ms²) = cosζ(wr(1,2,0), σs, ms²)
cosζ23_for0(σs, ms²) = cosζ(wr(2,3,0), σs, ms²)
cosζ31_for0(σs, ms²) = cosζ(wr(3,1,0), σs, ms²)

#
cosζk1_for1(k, σs, ms²) = cosζ(wr(k,1,1), σs, ms²)
cosζk2_for2(k, σs, ms²) = cosζ(wr(k,2,2), σs, ms²)
cosζk3_for3(k, σs, ms²) = cosζ(wr(k,3,3), σs, ms²)

cosθhatk1(k, σs, ms²) = cosζ(wr(k,1,0), σs, ms²)
cosθhatk2(k, σs, ms²) = cosζ(wr(k,2,0), σs, ms²)
cosθhatk3(k, σs, ms²) = cosζ(wr(k,3,0), σs, ms²)

"""
    Phase for wigner d-functions for clockwise rotations
"""
phase(two_λ1_minus_λ2) = (abs(two_λ1_minus_λ2) % 4 == 2 ? -1.0 : 1.0)
phase(two_λ1,two_λ2) = phase(two_λ1 - two_λ2)

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualNetworkUsage <: SwaggerModel
    currentValue::Any # spec type: Union{ Nothing, Float64 } # spec name: currentValue
    id::Any # spec type: Union{ Nothing, String } # spec name: id
    limit::Any # spec type: Union{ Nothing, Float64 } # spec name: limit
    name::Any # spec type: Union{ Nothing, VirtualNetworkUsageName } # spec name: name
    unit::Any # spec type: Union{ Nothing, String } # spec name: unit

    function VirtualNetworkUsage(;currentValue=nothing, id=nothing, limit=nothing, name=nothing, unit=nothing)
        o = new()
        validate_property(VirtualNetworkUsage, Symbol("currentValue"), currentValue)
        setfield!(o, Symbol("currentValue"), currentValue)
        validate_property(VirtualNetworkUsage, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(VirtualNetworkUsage, Symbol("limit"), limit)
        setfield!(o, Symbol("limit"), limit)
        validate_property(VirtualNetworkUsage, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(VirtualNetworkUsage, Symbol("unit"), unit)
        setfield!(o, Symbol("unit"), unit)
        o
    end
end # type VirtualNetworkUsage

const _property_map_VirtualNetworkUsage = Dict{Symbol,Symbol}(Symbol("currentValue")=>Symbol("currentValue"), Symbol("id")=>Symbol("id"), Symbol("limit")=>Symbol("limit"), Symbol("name")=>Symbol("name"), Symbol("unit")=>Symbol("unit"))
const _property_types_VirtualNetworkUsage = Dict{Symbol,String}(Symbol("currentValue")=>"Float64", Symbol("id")=>"String", Symbol("limit")=>"Float64", Symbol("name")=>"VirtualNetworkUsageName", Symbol("unit")=>"String")
Base.propertynames(::Type{ VirtualNetworkUsage }) = collect(keys(_property_map_VirtualNetworkUsage))
Swagger.property_type(::Type{ VirtualNetworkUsage }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_VirtualNetworkUsage[name]))}
Swagger.field_name(::Type{ VirtualNetworkUsage }, property_name::Symbol) =  _property_map_VirtualNetworkUsage[property_name]

function check_required(o::VirtualNetworkUsage)
    true
end

function validate_property(::Type{ VirtualNetworkUsage }, name::Symbol, val)
end


text:
struct TransferFunction{TE, S<:SisoTf{T} where T} <: LTISystem
    matrix::Matrix{S}
    timeevol::TE
    nu::Int
    ny::Int
    function TransferFunction{TE,S}(matrix::Matrix{S}, timeevol::TE) where {S,TE}
        # Validate size of input and output names
        ny, nu = size(matrix)
        return new{TE,S}(matrix, timeevol, nu, ny)
    end
end
function TransferFunction(matrix::Matrix{S}, timeevol::TE) where {TE<:TimeEvolution, T<:Number, S<:SisoTf{T}}
    TransferFunction{TE, S}(matrix, timeevol)
end

# # Constructor for Discrete time system
# function TransferFunction(matrix::Matrix{S}, Ts::Number) where {T<:Number, S<:SisoTf{T}}
#     return TransferFunction(matrix, Discrete(Ts))
# end
# # Constructor for Continuous time system
# function TransferFunction(matrix::Matrix{S}) where {T<:Number,S<:SisoTf{T}}
#     return TransferFunction(matrix, Continuous())
# end

noutputs(G::TransferFunction) = size(G.matrix, 1)
ninputs(G::TransferFunction) = size(G.matrix, 2)

#####################################################################
##                          Misc. Functions                        ##
#####################################################################

## INDEXING ##
Base.ndims(::TransferFunction) = 2
Base.size(G::TransferFunction) = size(G.matrix)
Base.eltype(::Type{S}) where {S<:TransferFunction} = S
Base.zero(G::TransferFunction{TE,S}) where {TE,S} = tf(zeros(numeric_type(S), size(G)), G.timeevol) # can not create a zero of a discrete system from the type alone, the sampletime is not stored.

function Base.getindex(G::TransferFunction{TE,S}, inds...) where {TE,S<:SisoTf}
    if size(inds, 1) != 2
        error("Must specify 2 indices to index TransferFunction model")
    end
    rows, cols = index2range(inds...)
    mat = Matrix{S}(undef, length(rows), length(cols))
    mat[:, :] = G.matrix[rows, cols]
    return TransferFunction(mat, G.timeevol)
end

function Base.copy(G::TransferFunction)
    return TransferFunction(copy(G.matrix), G.timeevol)
end

numvec(G::TransferFunction) = map(numvec, G.matrix)
denvec(G::TransferFunction) = map(denvec, G.matrix)

numpoly(G::TransferFunction) = map(numpoly, G.matrix)
denpoly(G::TransferFunction) = map(denpoly, G.matrix)

"""`tf = minreal(tf::TransferFunction, eps=sqrt(eps()))`

Create a minimial representation of each transfer function in `tf` by cancelling poles and zeros
will promote system to an appropriate numeric type"""
function minreal(G::TransferFunction, eps::Real=sqrt(eps()))
    matrix = similar(G.matrix, typeof(minreal(one(first(G.matrix)), eps)))
    for i = eachindex(G.matrix)
        matrix[i] = minreal(G.matrix[i], eps)
    end
    return TransferFunction(matrix, G.timeevol)
end

"""`isproper(tf)`

Returns `true` if the `TransferFunction` is proper. This means that order(den)
>= order(num))"""
function isproper(G::TransferFunction)
    return all(isproper(f) for f in G.matrix)
end

#####################################################################
##                         Math Operators                          ##
#####################################################################

## EQUALITY ##
function ==(G1::TransferFunction, G2::TransferFunction)
    fields = [:timeevol, :ny, :nu, :matrix]
    for field in fields
        if getfield(G1, field) != getfield(G2, field)
            return false
        end
    end
    return true
end

## Approximate ##
function isapprox(G1::TransferFunction, G2::TransferFunction; kwargs...)
    G1, G2 = promote(G1, G2)
    fieldsApprox = [:timeevol, :matrix]
    for field in fieldsApprox
        if !(isapprox(getfield(G1, field), getfield(G2, field); kwargs...))
            return false
        end
    end
    return true
end

function isapprox(G1::Array{T}, G2::Array{S}; kwargs...) where {T<:SisoTf, S<:SisoTf}
    all(i -> isapprox(promote(G1[i], G2[i])...; kwargs...), eachindex(G1))
end

## ADDITION ##
function +(G1::TransferFunction, G2::TransferFunction)
    if size(G1) != size(G2)
        error("Systems have different shapes.")
    end
    timeevol = common_timeevol(G1,G2)
    matrix = G1.matrix + G2.matrix
    return TransferFunction(matrix, timeevol)
end

+(G::TransferFunction, n::Number) = TransferFunction(G.matrix .+ n, G.timeevol)
+(n::Number, G::TransferFunction) = +(G, n)

## SUBTRACTION ##
-(n::Number, G::TransferFunction) = TransferFunction(n .- G.matrix, G.timeevol)
-(G1::TransferFunction, G2::TransferFunction) = +(G1, -G2)
-(G::TransferFunction, n::Number) = +(G, -n)

## NEGATION ##
-(G::TransferFunction) = TransferFunction(-G.matrix, G.timeevol)

## MULTIPLICATION ##

function *(G1::TransferFunction, G2::TransferFunction)
    # Note: G1*G2 = y <- G1 <- G2 <- u
    timeevol = common_timeevol(G1,G2)
    if G1.nu != G2.ny
        error("G1*G2: G1 must have same number of inputs as G2 has outputs")
    end
    matrix = G1.matrix * G2.matrix
    return TransferFunction(matrix, timeevol)
end

*(G::TransferFunction, n::Number) = TransferFunction(n*G.matrix, G.timeevol)
*(n::Number, G::TransferFunction) = *(G, n)

## DIVISION ##
function /(n::Number, G::TransferFunction)
    if issiso(G)
        matrix = reshape([n/G.matrix[1,1]], 1, 1)
    else
        error("MIMO TransferFunction inversion isn't implemented yet")
    end
    return TransferFunction(matrix, G.timeevol)
end
/(G::TransferFunction, n::Number) = G*(1/n)
/(G1::TransferFunction, G2::TransferFunction) = G1*(1/G2)
Base.:(/)(sys1::LTISystem, sys2::TransferFunction) = *(promote(sys1, ss(1/sys2))...) # This spcial case is needed to properly handle improper inverse transfer function (1/s)

#####################################################################
##                        Display Functions                        ##
#####################################################################

Base.print(io::IO, G::TransferFunction) = show(io, G)

function Base.show(io::IO, G::TransferFunction)
    # Compose the name vectors
    #println(io, "TransferFunction:")
    println(io, typeof(G))
    var = iscontinuous(G) ? :s : :z
    for i=1:G.nu
        for o=1:G.ny
            if !issiso(G)
                println(io, "Input $i to output $o")
            end
                print_siso(io, G.matrix[o, i], var)
            if !(i == G.nu && o == G.ny)
                print(io, "\n")
            end
        end
    end
    if iscontinuous(G)
        print(io, "\nContinuous-time transfer function model")
    elseif isdiscrete(G)
        print(io, "\nSample Time: ", G.Ts, " (seconds)")
        print(io, "\nDiscrete-time transfer function model")
    else
        print(io, "\nStatic gain transfer function model")
    end
end


text:
#
# Correctness Tests
#

const io = IOBuffer()
logging(io, kind=:warn)
logging(io, kind=:info)

using Compat
#import Compat.Sys: isapple
import Compat: stdin, @info
using Compat.MathConstants
using Compat.DelimitedFiles
using Compat.Test
using GeoEfficiency
const G = GeoEfficiency

function exec_consol_unattended(Fn::Union{Function,Type}, consol_inputs::Vector; Fn_ARGs::Vector=[])
	for input in  string.(consol_inputs)
		 write(stdin.buffer, input,"\n")
	end
	return Fn(Fn_ARGs...)
end
exec_consol_unattended(Fn::Union{Function,Type}, consol_inputs...; Fn_ARGs::Vector=[]) = exec_consol_unattended(Fn, consol_inputs; Fn_ARGs=Fn_ARGs)
exec_consol_unattended(Fn::Union{Function,Type}, consol_inputs::String; Fn_ARGs::Vector=[]) = exec_consol_unattended(Fn, split(consol_inputs); Fn_ARGs=Fn_ARGs)

tests = ["Input_Interface",
         "Physics_Model",
         "Calculations",
         "Output_Interface"]

println("\nRunning tests:")

for t in tests
	println(); @info("Begin test of $(t).....\n")
    include("test_$(t).jl")
	println(); @info("End test of $(t).....\n")
end
@test about() == nothing


text:
# Test the extensions of Base.convert
@testset "Convert Extensions" begin
    m = Model()
    @variable(m, x[1:2])
    # Array -> SparseAxisArray
    @testset "Array->Sparse" begin
        @test convert(JuMPC.SparseAxisArray,
                      x) isa JuMPC.SparseAxisArray
        @test length(convert(JuMPC.SparseAxisArray, x)) == length(x)
        inds = CartesianIndices(x)
        @test convert(JuMPC.SparseAxisArray,
                      x)[inds[1][1]] == x[inds[1]]
        @test convert(JuMPC.SparseAxisArray,
                      x)[inds[2][1]] == x[inds[2]]
    end
    # DenseAxisArray -> SparseAxisArray
    @testset "DenseArray->Sparse" begin
        @variable(m, y[2:3])
        @test convert(JuMPC.SparseAxisArray,
                      y) isa JuMPC.SparseAxisArray
        @test length(convert(JuMPC.SparseAxisArray, y)) == length(y)
        inds = collect(keys(y))
        @test convert(JuMPC.SparseAxisArray,
                      y)[inds[1][1]] == y[inds[1]]
        @test convert(JuMPC.SparseAxisArray,
                      y)[inds[2][1]] == y[inds[2]]
    end
    # Number -> AffExpr
    @testset "Number->AffExpr" begin
        expr = zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef}) + 42
        @test convert(JuMP.AbstractJuMPScalar, 42) == expr
    end
end

# Test extension of keys
@testset "_keys" begin
    m = Model()
    @variable(m, x[1:2], container = SparseAxisArray)
    @test InfiniteOpt._keys(x) == keys(x.data)
    @test InfiniteOpt._keys(ones(2)) == keys(ones(2))
end

# Test extension of map
@testset "Base.map (DenseAxisArray)" begin
    JuMPC.@container(x[3:6], 2)
    @test map(i -> iseven(i), x) isa JuMPC.DenseAxisArray
    @test map(i -> iseven(i), x).data == [true, true, true, true]
end

# Test the _make_vector functions
@testset "_make_vector" begin
    x = [2, 4, 1, 8]
    # test with JuMP containers
    @testset "AbstractArray" begin
        y = convert(JuMPC.SparseAxisArray, x)
        @test sort(InfiniteOpt._make_vector(y)) == [1, 2, 4, 8]
        y = JuMPC.DenseAxisArray(x, [1, 2, 3, 4])
        @test InfiniteOpt._make_vector(y) == x
        @test InfiniteOpt._make_vector(zeros(2, 0)) == []
    end
    # test arrays
    @testset "Array" begin
        @test InfiniteOpt._make_vector(x) == x
        @test InfiniteOpt._make_vector(ones(2, 2)) == ones(4)
    end
end

# Test the _make_ordered_vector functions
@testset "_make_ordered_vector" begin
    x = [2, 4, 1, 8]
    # test with JuMP containers
    @testset "Vector" begin
        x = ones(2)
        @test InfiniteOpt._make_ordered_vector(x) === x
    end
    # test arrays
    @testset "Array" begin
        x = ones(2, 3)
        @test InfiniteOpt._make_ordered_vector(x) == ones(6)
        JuMPC.@container(x[3:6], 2)
        @test InfiniteOpt._make_ordered_vector(x) == ones(4) * 2
    end
    # test SparseAxisArray
    @testset "SparseAxisArray" begin
        JuMPC.@container(x[i = 3:6], i, container = SparseAxisArray)
        @test InfiniteOpt._make_ordered_vector(x) == [3, 4, 5, 6]
    end
    # test fallback
    @testset "Fallback" begin
    @test InfiniteOpt._make_ordered_vector(2) == 2
    end
end

# Test allequal
@testset "_allequal" begin
    @test InfiniteOpt._allequal([1])
    @test !InfiniteOpt._allequal([1, 2])
    @test InfiniteOpt._allequal([1, 1])
    @test InfiniteOpt._allequal([1 1; 1 1])
end

# Test _make_float_info
@testset "_make_float_info" begin
    info1 = VariableInfo(false, 0, false, 0, false, 0, false, 0, false, false)
    num = Float64(0)
    info2 = VariableInfo(true, num, true, num, true, num, true, num, true, true)
    @test InfiniteOpt._make_float_info(info1) isa VariableInfo{Float64, Float64, Float64, Float64}
    @test InfiniteOpt._make_float_info(info2) == info2
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export dashboarditem

"""
    dashboarditem(;kwargs...)
    dashboarditem(children::Any;kwargs...)
    dashboarditem(children_maker::Function;kwargs...)

A DashboardItem component.
DashboardItem is a wrapper that is intended to be used with DraggableDashboard
to specify the position and size of the wrapped components.
Keyword arguments:
- `children` (Array of a list of or a singular dash component, string or numbers | a list of or a singular dash component, string or number; optional): The child or list of children wrapped by the component.
- `h` (Real; optional): The height on the of y axis (default is 4)
- `i` (String; optional): The ID used to identify this component in layout definition.
- `isBounded` (Bool; optional)
- `isDraggable` (Bool; optional): If false, will not be draggable. Overrides `static`.
- `isResizable` (Bool; optional): If false, will not be resizable. Overrides `static`.
- `maxH` (Real; optional)
- `maxW` (Real; optional)
- `minH` (Real; optional)
- `minW` (Real; optional)
- `moved` (Bool; optional)
- `resizeHandles` (String; optional)
- `static` (Bool; optional): If true, equal to `isDraggable: false, isResizable: false`.
- `w` (Real; optional): The width of the x axis (default is 6).
- `x` (Real; optional): The position on the x axis in number of columns (by default, the  max is 12).
- `y` (Real; optional): The position on the y axis (the unit is 30px, by default)
"""
function dashboarditem(; kwargs...)
        available_props = Symbol[:children, :h, :i, :isBounded, :isDraggable, :isResizable, :maxH, :maxW, :minH, :minW, :moved, :resizeHandles, :static, :w, :x, :y]
        wild_props = Symbol[]
        return Component("dashboarditem", "DashboardItem", "dash_draggable", available_props, wild_props; kwargs...)
end

dashboarditem(children::Any; kwargs...) = dashboarditem(;kwargs..., children = children)
dashboarditem(children_maker::Function; kwargs...) = dashboarditem(children_maker(); kwargs...)

max_stars_repo_name: tcgvn/dash-draggable

max_issues_repo_name: tcgvn/dash-draggable

max_forks_repo_name: tcgvn/dash-draggable


text:
Random.seed!(0)

#sig = [1 2; 2 2]
#sig = [1 1; 2 3]
sig = [2 3]
#sig = [2 2]
#sig = [3 2; 2 3]

S = random_S0Graph(sig)

w = random_bounded(S.n)

U = random_S1_unitary(sig);
SU = S0Graph(S.sig, U * S.S * U')
@time opt1 = dsw(S, w, eps=solver_eps)[1]
@time opt2 = dsw(SU, U*w*U', eps=solver_eps)[1]
@test opt1 ≈ opt2  atol=tol

max_stars_repo_name: dstahlke/NoncommutativeGraphs

max_issues_repo_name: dstahlke/NoncommutativeGraphs

max_forks_repo_name: dstahlke/NoncommutativeGraphs


text:
using Test
# Note, try run_process_test.jl
@testset "1" begin
    include("test_1.jl")
end
@testset "2" begin
    include("test_2.jl")
end
@testset "3" begin let
    include("test_3.jl")
end end
@testset "4" begin let
    include("test_4.jl")
end end
@testset "5" begin let
    include("test_5.jl")
end end
@testset "6" begin let
    include("test_6.jl")
end end
@testset "7" begin let
    include("test_7.jl")
end end
@testset "8" begin let
    include("test_8.jl")
end end
@testset "9" begin let
    include("test_9.jl")
end end
@testset "10" begin let
    include("test_10.jl")
end end
@testset "11" begin let
    include("test_11.jl")
end end
@testset "12" begin let
    include("test_12.jl")
end end
@testset "13" begin let
    include("test_13.jl")
end end
@testset "14" begin let
    include("test_14.jl")
end end
@testset "15" begin let
    include("test_15.jl")
end end
@testset "16" begin let
    include("test_16.jl")
end end
@testset "17" begin let
    include("test_17.jl")
end end
@testset "18" begin let
    include("test_18.jl")
end end
@testset "19" begin let
    include("test_19.jl")
end end
@time @testset "20" begin let
    include("test_20.jl")
end end
@time @testset "21" begin let
    include("test_21.jl")
end end
@time @testset "22" begin let
    include("test_22.jl")
end end
@time @testset "23" begin let
    include("test_23.jl")
end end
#=
@time @testset "24" begin let
    include("test_24.jl")
end end
=#
@time @testset "25" begin let
    include("test_25.jl")
end end
@time @testset "26" begin let
    include("test_26.jl")
end end
@time @testset "27" begin let
    include("test_27.jl")
end end
@testset "28" begin let
    include("test_28.jl")
end end
@testset "29" begin let
    include("test_29.jl")
end end
@testset "30" begin let
    include("test_30.jl")
end end
@testset "31" begin let
    include("test_31.jl")
end end
@testset "32" begin let
    include("test_32.jl")
end end
@testset "33" begin let
    include("test_33.jl")
end end
@testset "34" begin let
    include("test_34.jl")
end end
@testset "35" begin let
    include("test_35.jl")
end end
@testset "36" begin let
    include("test_36.jl")
end end
@testset "37" begin let
    include("test_37.jl")
end end


text:
function MullerSolver()
    return MullerSolver(())
end

function MullerSolver(arg0::jdouble)
    return MullerSolver((jdouble,), arg0)
end

function MullerSolver(arg0::jdouble, arg1::jdouble)
    return MullerSolver((jdouble, jdouble), arg0, arg1)
end


text:
"""
    AbstractCartanSchoutenConnection

Abstract type for Cartan-Schouten connections, that is connections whose geodesics
going through group identity are one-parameter subgroups. See [^Pennec2020] for details.

[^Pennec2020]:
    > X. Pennec and M. Lorenzi, “5 - Beyond Riemannian geometry: The affine connection
    > setting for transformation groups,” in Riemannian Geometric Statistics in Medical Image
    > Analysis, X. Pennec, S. Sommer, and T. Fletcher, Eds. Academic Press, 2020, pp. 169–229.
    > doi: 10.1016/B978-0-12-814725-2.00012-1.
"""
abstract type AbstractCartanSchoutenConnection <: AbstractAffineConnection end

"""
    CartanSchoutenMinus

The unique Cartan-Schouten connection such that all left-invariant vector fields are
globally defined by their value at identity. It is biinvariant with respect to the group
operation.
"""
struct CartanSchoutenMinus <: AbstractCartanSchoutenConnection end

"""
    CartanSchoutenPlus

The unique Cartan-Schouten connection such that all right-invariant vector fields are
globally defined by their value at identity. It is biinvariant with respect to the group
operation.
"""
struct CartanSchoutenPlus <: AbstractCartanSchoutenConnection end

"""
    CartanSchoutenZero

The unique torsion-free Cartan-Schouten connection. It is biinvariant with respect to the
group operation.

If the metric on the underlying manifold is bi-invariant then it is equivalent to the
Levi-Civita connection of that metric.
"""
struct CartanSchoutenZero <: AbstractCartanSchoutenConnection end

const CartanSchoutenMinusGroup{𝔽,M} = ConnectionManifold{𝔽,M,CartanSchoutenMinus}
const CartanSchoutenPlusGroup{𝔽,M} = ConnectionManifold{𝔽,M,CartanSchoutenPlus}
const CartanSchoutenZeroGroup{𝔽,M} = ConnectionManifold{𝔽,M,CartanSchoutenZero}

"""
    exp!(M::ConnectionManifold{𝔽,<:AbstractDecoratorManifold{𝔽},<:AbstractCartanSchoutenConnection}, q, p, X) where {𝔽}

Compute the exponential map on the [`ConnectionManifold`](@ref) `M` with a Cartan-Schouten
connection. See Sections 5.3.2 and 5.3.3 of [^Pennec2020] for details.

[^Pennec2020]:
    > X. Pennec and M. Lorenzi, “5 - Beyond Riemannian geometry: The affine connection
    > setting for transformation groups,” in Riemannian Geometric Statistics in Medical Image
    > Analysis, X. Pennec, S. Sommer, and T. Fletcher, Eds. Academic Press, 2020, pp. 169–229.
    > doi: 10.1016/B978-0-12-814725-2.00012-1.
"""
function exp!(
    M::ConnectionManifold{
        𝔽,
        <:AbstractDecoratorManifold{𝔽},
        <:AbstractCartanSchoutenConnection,
    },
    q,
    p,
    X,
) where {𝔽}
    Y = inverse_translate_diff(M.manifold, q, p, X)
    return compose!(M.manifold, q, p, exp_lie(M.manifold, Y))
end

"""
    log!(M::ConnectionManifold{𝔽,<:AbstractDecoratorManifold{𝔽},<:AbstractCartanSchoutenConnection}, Y, p, q) where {𝔽}

Compute the logarithmic map on the [`ConnectionManifold`](@ref) `M` with a Cartan-Schouten
connection. See Sections 5.3.2 and 5.3.3 of [^Pennec2020] for details.

[^Pennec2020]:
    > X. Pennec and M. Lorenzi, “5 - Beyond Riemannian geometry: The affine connection
    > setting for transformation groups,” in Riemannian Geometric Statistics in Medical Image
    > Analysis, X. Pennec, S. Sommer, and T. Fletcher, Eds. Academic Press, 2020, pp. 169–229.
    > doi: 10.1016/B978-0-12-814725-2.00012-1.
"""
function log!(
    M::ConnectionManifold{
        𝔽,
        <:AbstractDecoratorManifold{𝔽},
        <:AbstractCartanSchoutenConnection,
    },
    Y,
    p,
    q,
) where {𝔽}
    pinvq = compose(M.manifold, inv(M.manifold, p), q)
    log_lie!(M.manifold, Y, pinvq)
    return translate_diff!(M.manifold, Y, p, Identity(M.manifold), Y)
end

"""
    parallel_transport_to(M::CartanSchoutenMinusGroup, p, X, q)

Transport tangent vector `X` at point `p` on the group manifold `M` with the
[`CartanSchoutenMinus`](@ref) connection to point `q`. See [^Pennec2020] for details.

[^Pennec2020]:
    > X. Pennec and M. Lorenzi, “5 - Beyond Riemannian geometry: The affine connection
    > setting for transformation groups,” in Riemannian Geometric Statistics in Medical Image
    > Analysis, X. Pennec, S. Sommer, and T. Fletcher, Eds. Academic Press, 2020, pp. 169–229.
    > doi: 10.1016/B978-0-12-814725-2.00012-1.
"""
function parallel_transport_to(M::CartanSchoutenMinusGroup, p, X, q)
    return inverse_translate_diff(M.manifold, q, p, X, LeftAction())
end

function parallel_transport_to!(M::CartanSchoutenMinusGroup, Y, p, X, q)
    return inverse_translate_diff!(M.manifold, Y, q, p, X, LeftAction())
end

"""
    vector_transport_to(M::CartanSchoutenPlusGroup, p, X, q)

Transport tangent vector `X` at point `p` on the group manifold `M` with the
[`CartanSchoutenPlus`](@ref) connection to point `q`. See [^Pennec2020] for details.

[^Pennec2020]:
    > X. Pennec and M. Lorenzi, “5 - Beyond Riemannian geometry: The affine connection
    > setting for transformation groups,” in Riemannian Geometric Statistics in Medical Image
    > Analysis, X. Pennec, S. Sommer, and T. Fletcher, Eds. Academic Press, 2020, pp. 169–229.
    > doi: 10.1016/B978-0-12-814725-2.00012-1.
"""
parallel_transport_to(M::CartanSchoutenPlusGroup, p, X, q)

function parallel_transport_to!(M::CartanSchoutenPlusGroup, Y, p, X, q)
    return inverse_translate_diff!(M.manifold, Y, q, p, X, RightAction())
end

"""
    parallel_transport_direction(M::CartanSchoutenZeroGroup, ::Identity, X, d)

Transport tangent vector `X` at identity on the group manifold with the
[`CartanSchoutenZero`](@ref) connection in the direction `d`. See [^Pennec2020] for details.

[^Pennec2020]:
    > X. Pennec and M. Lorenzi, “5 - Beyond Riemannian geometry: The affine connection
    > setting for transformation groups,” in Riemannian Geometric Statistics in Medical Image
    > Analysis, X. Pennec, S. Sommer, and T. Fletcher, Eds. Academic Press, 2020, pp. 169–229.
    > doi: 10.1016/B978-0-12-814725-2.00012-1.
"""
function parallel_transport_direction(M::CartanSchoutenZeroGroup, p::Identity, X, d)
    dexp_half = exp_lie(M.manifold, d / 2)
    Y = translate_diff(M.manifold, dexp_half, p, X, RightAction())
    return translate_diff(M.manifold, dexp_half, p, Y, LeftAction())
end

function parallel_transport_direction!(M::CartanSchoutenZeroGroup, Y, p::Identity, X, d)
    dexp_half = exp_lie(M.manifold, d / 2)
    translate_diff!(M.manifold, Y, dexp_half, p, X, RightAction())
    return translate_diff!(M.manifold, Y, dexp_half, p, Y, LeftAction())
end

"""
    parallel_transport_to(M::CartanSchoutenZeroGroup, p::Identity, X, q)

Transport vector `X` at identity of group `M` equipped with the [`CartanSchoutenZero`](@ref)
connection to point `q` using parallel transport.
"""
function parallel_transport_to(M::CartanSchoutenZeroGroup, p::Identity, X, q)
    d = log_lie(M.manifold, q)
    return parallel_transport_direction(M, p, X, d)
end

function parallel_transport_to!(M::CartanSchoutenZeroGroup, Y, p::Identity, X, q)
    d = log_lie(M.manifold, q)
    return parallel_transport_direction!(M, Y, p, X, d)
end

max_stars_repo_name: JuliaNLSolvers/ManifoldMuseum

max_issues_repo_name: JuliaNLSolvers/ManifoldMuseum

max_forks_repo_name: JuliaNLSolvers/ManifoldMuseum