text:
import Base: isempty,
             rand,
             convert

export HPolyhedron,
       dim, σ, ∈,
       addconstraint!,
       constraints_list,
       tohrep, tovrep,
       convex_hull,
       cartesian_product,
       isempty,
       remove_redundant_constraints,
       remove_redundant_constraints!,
       constrained_dimensions,
       is_hyperplanar

"""
    HPolyhedron{N, VN<:AbstractVector{N}} <: AbstractPolyhedron{N}

Type that represents a convex polyhedron in H-representation, that is a finite intersection of half-spaces,
```math
P = \\bigcap_{i = 1}^m H_i,
```
where each ``H_i = \\{x \\in \\mathbb{R}^n : a_i^T x \\leq b_i \\}`` is a half-space,
``a_i \\in \\mathbb{R}^n`` is the normal vector of the ``i``-th half-space and ``b_i`` is the displacement.
The set ``P`` may or may not be bounded (see also [`HPolytope`](@ref), which assumes boundedness).

### Fields

- `constraints` -- vector of linear constraints
"""
struct HPolyhedron{N, VN<:AbstractVector{N}} <: AbstractPolyhedron{N}
    constraints::Vector{LinearConstraint{N, VN}}

    function HPolyhedron(constraints::Vector{LinearConstraint{N, VN}}) where {N,
                                                                              VN<:AbstractVector{N}}
        return new{N, VN}(constraints)
    end
end

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

# constructor for an HPolyhedron with no constraints
function HPolyhedron{N, VN}() where {N, VN<:AbstractVector{N}}
    HPolyhedron(Vector{LinearConstraint{N, VN}}())
end

# constructor for an HPolyhedron with no constraints and given numeric type
function HPolyhedron{N}() where {N}
    HPolyhedron(Vector{LinearConstraint{N, Vector{N}}}())
end

# constructor for an HPolyhedron without explicit numeric type, defaults to Float64
function HPolyhedron()
    HPolyhedron{Float64}()
end

# constructor for an HPolyhedron with constraints of mixed type
function HPolyhedron(constraints::Vector{<:LinearConstraint})
    HPolyhedron(_normal_Vector(constraints))
end

# constructor from a simple H-representation
HPolyhedron(A::AbstractMatrix, b::AbstractVector) =
    HPolyhedron(constraints_list(A, b))

# convenience union type
const HPoly{N} = Union{HPolytope{N}, HPolyhedron{N}}

# --- LazySet interface functions ---

"""
    dim(P::HPoly)

Return the dimension of a polyhedron in H-representation.

### Input

- `P`  -- polyhedron in H-representation

### Output

The ambient dimension of the polyhedron in H-representation.
If it has no constraints, the result is ``-1``.
"""
function dim(P::HPoly)
    return length(P.constraints) == 0 ? -1 : length(P.constraints[1].a)
end

"""
    ρ(d::AbstractVector{M}, P::HPoly{N};
      solver=default_lp_solver(M, N)) where {M, N}

Evaluate the support function of a polyhedron (in H-representation) in a given
direction.

### Input

- `d`      -- direction
- `P`      -- polyhedron in H-representation
- `solver` -- (optional, default: `default_lp_solver(M, N)`) the backend used to
              solve the linear program

### Output

The support function of the polyhedron.
If a polytope is unbounded in the given direction, we throw an error.
If a polyhedron is unbounded in the given direction, the result is `Inf`.
"""
function ρ(d::AbstractVector{M}, P::HPoly{N};
           solver=default_lp_solver(M, N)) where {M, N}
    lp, unbounded = σ_helper(d, P, solver)
    if unbounded
        if P isa HPolytope
            error("the support function in direction $(d) is undefined " *
                  "because the polytope is unbounded")
        end
        return N(Inf)
    end
    return dot(d, lp.sol)
end

"""
    σ(d::AbstractVector{M}, P::HPoly{N};
      solver=default_lp_solver(M, N) where {M, N}

Return the support vector of a polyhedron (in H-representation) in a given
direction.

### Input

- `d`      -- direction
- `P`      -- polyhedron in H-representation
- `solver` -- (optional, default: `default_lp_solver(M, N)`) the backend used to
              solve the linear program

### Output

The support vector in the given direction.
"""
function σ(d::AbstractVector{M}, P::HPoly{N};
           solver=default_lp_solver(M, N)) where {M, N}
    lp, unbounded = σ_helper(d, P, solver)
    if unbounded
        if P isa HPolytope
            error("the support vector in direction $(d) is undefined because " *
                  "the polytope is unbounded")
        end
        # construct the solution from the solver's ray result
        if lp == nothing
            ray = d
        elseif has_lp_infeasibility_ray(lp.model)
            ray = lp.sol  # infeasibility ray is stored as the solution
        else
            error("LP solver did not return an infeasibility ray")
        end
        res = Vector{N}(undef, length(ray))
        @inbounds for i in 1:length(ray)
            if ray[i] == zero(N)
                res[i] = zero(N)
            elseif ray[i] > zero(N)
                res[i] = N(Inf)
            else
                res[i] = N(-Inf)
            end
        end
        return res
    else
        return lp.sol
    end
end

function σ_helper(d::AbstractVector, P::HPoly, solver)
    # let c = -d as a Vector since GLPK does not accept sparse vectors
    # (see #1011)
    c = to_negative_vector(d)

    (A, b) = tosimplehrep(P)
    if length(b) == 0
        unbounded = true
        lp = nothing
    else
        sense = '<'
        l = -Inf
        u = Inf
        lp = linprog(c, A, sense, b, l, u, solver)
        if is_lp_infeasible(lp.status, strict=true)
            error("the support vector is undefined because the polyhedron is " *
                  "empty")
        elseif is_lp_unbounded(lp.status)
            unbounded = true
        elseif is_lp_optimal(lp.status)
            unbounded = false
        else
            error("got unknown LP status $(lp.status)")
        end
    end
    return (lp, unbounded)
end

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

Create a polyhedron.

### Input

- `HPolyhedron` -- type for dispatch
- `N`           -- (optional, default: `Float64`) numeric type
- `dim`         -- (optional, default: 2) dimension (is ignored)
- `rng`         -- (optional, default: `GLOBAL_RNG`) random number generator
- `seed`        -- (optional, default: `nothing`) seed for reseeding

### Output

A polyhedron.

### Algorithm

We first create a random polytope and then randomly remove some of the
constraints.
"""
function rand(::Type{HPolyhedron};
              N::Type{<:Real}=Float64,
              dim::Int=2,
              rng::AbstractRNG=GLOBAL_RNG,
              seed::Union{Int, Nothing}=nothing)
    rng = reseed(rng, seed)
    P = rand(HPolytope; N=N, dim=dim, rng=rng)
    constraints_P = constraints_list(P)
    constraints_Q = Vector{eltype(constraints_P)}()
    for i in 1:length(constraints_P)
        if rand(Bool)
            push!(constraints_Q, constraints_P[i])
        end
    end
    return HPolyhedron(constraints_Q)
end

# ===========================================
# HPolyhedron and HPolytope's shared methods
# ===========================================

"""
    addconstraint!(P::HPoly, constraint::LinearConstraint)

Add a linear constraint to a polyhedron in H-representation.

### Input

- `P`          -- polyhedron in H-representation
- `constraint` -- linear constraint to add

### Notes

It is left to the user to guarantee that the dimension of all linear constraints
is the same.
"""
function addconstraint!(P::HPoly, constraint::LinearConstraint)
    push!(P.constraints, constraint)
    return nothing
end

"""
    constraints_list(P::HPoly)

Return the list of constraints defining a polyhedron in H-representation.

### Input

- `P` -- polyhedron in H-representation

### Output

The list of constraints of the polyhedron.
"""
function constraints_list(P::HPoly)
    return P.constraints
end

"""
    tohrep(P::HPoly)

Return a constraint representation of the given polyhedron in constraint
representation (no-op).

### Input

- `P` -- polyhedron in constraint representation

### Output

The same polyhedron instance.
"""
function tohrep(P::HPoly)
    return P
end

"""
    normalize(P::HPoly{N}, p=N(2)) where {N}

Normalize a polyhedron in constraint representation.

### Input

- `P` -- polyhedron in constraint representation
- `p` -- (optional, default: `2`) norm

### Output

A new polyhedron in constraint representation whose normal directions ``a_i``
are normalized, i.e., such that ``‖a_i‖_p = 1`` holds.
"""
function normalize(P::HPoly{N}, p=N(2)) where {N}
    constraints = [normalize(hs, p) for hs in constraints_list(P)]
    T = basetype(P)
    return T(constraints)
end

"""
    remove_redundant_constraints(P::HPoly{N};
                                 backend=default_lp_solver(N)) where {N}

Remove the redundant constraints in a polyhedron in H-representation.

### Input

- `P`       -- polyhedron
- `backend` -- (optional, default: `default_lp_solver(N)`) the backend used to
               solve the linear program

### Output

A polyhedron equivalent to `P` but with no redundant constraints, or an empty
set if `P` is detected to be empty, which may happen if the constraints are
infeasible.

### Algorithm

See `remove_redundant_constraints!(::AbstractVector{<:LinearConstraint})` for
details.
"""
function remove_redundant_constraints(P::HPoly{N};
                                      backend=default_lp_solver(N)) where {N}
    Pred = copy(P)
    if remove_redundant_constraints!(Pred, backend=backend)
        return Pred
    else # the polyhedron P is empty
        return EmptySet{N}(dim(P))
    end
end

"""
    remove_redundant_constraints!(P::HPoly{N};
                                  backend=default_lp_solver(N)) where {N}

Remove the redundant constraints in a polyhedron in H-representation; the
polyhedron is updated in-place.

### Input

- `P`       -- polyhedron
- `backend` -- (optional, default: `default_lp_solver(N)`) the backend used to
               solve the linear program

### Output

`true` if the method was successful and the polyhedron `P` is modified by
removing its redundant constraints, and `false` if `P` is detected to be empty,
which may happen if the constraints are infeasible.

### Algorithm

See `remove_redundant_constraints!(::AbstractVector{<:LinearConstraint})` for
details.
"""
function remove_redundant_constraints!(P::HPoly{N};
                                       backend=default_lp_solver(N)) where {N}
    remove_redundant_constraints!(P.constraints, backend=backend)
end

"""
    translate(P::HPoly, v::AbstractVector; [share]::Bool=false)

Translate (i.e., shift) a polyhedron in constraint representation by a given
vector.

### Input

- `P`     -- polyhedron in constraint representation
- `v`     -- translation vector
- `share` -- (optional, default: `false`) flag for sharing unmodified parts of
             the original set representation

### Output

A translated polyhedron in constraint representation.

### Notes

The normal vectors of the constraints (vector `a` in `a⋅x ≤ b`) are shared with
the original constraints if `share == true`.

### Algorithm

We translate every constraint.
"""
function translate(P::HPoly, v::AbstractVector; share::Bool=false)
    @assert length(v) == dim(P) "cannot translate a $(dim(P))-dimensional " *
                                "set by a $(length(v))-dimensional vector"
    constraints = [translate(c, v; share=share) for c in constraints_list(P)]
    T = basetype(P)
    return T(constraints)
end

# ========================================================
# External methods that require Polyhedra.jl to be loaded
# ========================================================

"""
    tovrep(P::HPoly;
          [backend]=default_polyhedra_backend(P))

Transform a polyhedron in H-representation to a polytope in V-representation.

### Input

- `P`       -- polyhedron in constraint representation
- `backend` -- (optional, default: `default_polyhedra_backend(P)`) the
               backend for polyhedral computations

### Output

The `VPolytope` which is the vertex representation of the given polyhedron
in constraint representation.

### Notes

The conversion may not preserve the numeric type (e.g., with `N == Float32`)
depending on the backend.
For further information on the supported backends see [Polyhedra's
documentation](https://juliapolyhedra.github.io/).
"""
function tovrep(P::HPoly;
                backend=default_polyhedra_backend(P))
    require(:Polyhedra; fun_name="tovrep")
    P = polyhedron(P; backend=backend)
    return VPolytope(P)
end

"""
    isempty(P::HPoly{N}, witness::Bool=false;
            [use_polyhedra_interface]::Bool=false, [solver]=nothing,
            [backend]=nothing) where {N}

Determine whether a polyhedron is empty.

### Input

- `P`       -- polyhedron
- `witness` -- (optional, default: `false`) compute a witness if activated
- `use_polyhedra_interface` -- (optional, default: `false`) if `true`, we use
               the `Polyhedra` interface for the emptiness test
- `solver`  -- (optional, default: `nothing`) LP-solver backend, uses `default_lp_solver(N)`
               if it is not provided
- `backend` -- (optional, default: `nothing`) backend for polyhedral
               computations in `Polyhedra`, uses `default_polyhedra_backend(P)` if
               it is not provided

### Output

* If `witness` option is deactivated: `true` iff ``P = ∅``
* If `witness` option is activated:
  * `(true, [])` iff ``P = ∅``
  * `(false, v)` iff ``P ≠ ∅`` and ``v ∈ P``

### Notes

The default value of the `backend` is set internally and depends on whether the
`use_polyhedra_interface` option is set or not.
If the option is set, we use `default_polyhedra_backend(P)`.

Witness production is not supported if `use_polyhedra_interface` is `true`.

### Algorithm

The algorithm sets up a feasibility LP for the constraints of `P`.
If `use_polyhedra_interface` is `true`, we call `Polyhedra.isempty`.
Otherwise, we set up the LP internally.
"""
function isempty(P::HPoly{N},
                 witness::Bool=false;
                 use_polyhedra_interface::Bool=false,
                 solver=nothing,
                 backend=nothing) where {N}
    if length(constraints_list(P)) < 2
        # catch corner case because of problems in LP solver for Rationals
        return witness ? (false, an_element(P)) : false
    end
    if use_polyhedra_interface
        require(:Polyhedra; fun_name="isempty", explanation="with the active " *
            "option `use_polyhedra_interface`")

        if backend == nothing
            backend = default_polyhedra_backend(P)
        end

        if isnothing(solver)
            result = Polyhedra.isempty(polyhedron(P; backend=backend))
        else
            result = Polyhedra.isempty(polyhedron(P; backend=backend), solver)
        end

        if result
            return witness ? (true, N[]) : true
        elseif witness
            error("witness production is not supported yet")
        else
            return false
        end
    else
        A, b = tosimplehrep(P)
        lbounds, ubounds = -Inf, Inf
        sense = '<'
        obj = zeros(N, size(A, 2))
        if isnothing(solver)
            solver = default_lp_solver(N)
        end
        lp = linprog(obj, A, sense, b, lbounds, ubounds, solver)
        if is_lp_optimal(lp.status)
            return witness ? (false, lp.sol) : false
        elseif is_lp_infeasible(lp.status)
            return witness ? (true, N[]) : true
        end
        error("LP returned status $(lp.status) unexpectedly")
    end
end

# ==================================
# Methods that use Polyhedra.jl
# ==================================

function load_polyhedra_hpolyhedron() # function to be loaded by Requires
return quote
# see the interface file init_Polyhedra.jl for the imports

function convert(::Type{HPolyhedron}, P::HRep{N}) where {N}
    VN = Polyhedra.hvectortype(P)
    constraints = Vector{LinearConstraint{N, VN}}()
    for hi in Polyhedra.allhalfspaces(P)
        a, b = hi.a, hi.β
        if isapproxzero(norm(a))
            continue
        end
        push!(constraints, HalfSpace(a, b))
    end
    return HPolyhedron(constraints)
end

"""
     HPolyhedron(P::HRep{N}) where {N}

Return a polyhedron in H-representation given a `HRep` polyhedron
from `Polyhedra.jl`.

### Input

- `P` -- `HRep` polyhedron

### Output

An `HPolyhedron`.
"""
function HPolyhedron(P::HRep{N}) where {N}
    convert(HPolyhedron, P)
end

"""
    polyhedron(P::HPoly; [backend]=default_polyhedra_backend(P))

Return an `HRep` polyhedron from `Polyhedra.jl` given a polytope in
H-representation.

### Input

- `P`       -- polytope
- `backend` -- (optional, default: call `default_polyhedra_backend(P)`)
                the polyhedral computations backend

### Output

An `HRep` polyhedron.

### Notes

For further information on the supported backends see
[Polyhedra's documentation](https://juliapolyhedra.github.io/).
"""
function polyhedron(P::HPoly;
                    backend=default_polyhedra_backend(P))
    A, b = tosimplehrep(P)
    return Polyhedra.polyhedron(Polyhedra.hrep(A, b), backend)
end

function triangulate(X::LazySet)

    dim(X) == 3 || throw(ArgumentError("the dimension of the set should be three, got $(dim(X))"))

    poly = polyhedron(convert(HPolyhedron, X))
    mes = Mesh(poly)
    coords = Polyhedra.GeometryBasics.coordinates(mes)
    connec = Polyhedra.GeometryBasics.faces(mes)

    ntriangles = length(connec)
    npoints = 3*ntriangles
    points = Matrix{Float32}(undef, 3, npoints)

    for i in 1:npoints
        points[:, i] .= coords[i].data
    end

    connec_tup = getfield.(connec, :data)

    return points, connec_tup
end

end # quote
end # function load_polyhedra_hpolyhedron()

"""
    _isbounded_stiemke(P::HPolyhedron{N}; solver=LazySets.default_lp_solver(N),
                       check_nonempty::Bool=true) where {N}

Determine whether a polyhedron is bounded using Stiemke's theorem of alternatives.

### Input

- `P`       -- polyhedron
- `backend` -- (optional, default: `default_lp_solver(N)`) the backend used
               to solve the linear program
- `check_nonempty` -- (optional, default: `true`) if `true`, check the
                      precondition to this algorithm that `P` is non-empty

### Output

`true` iff the polyhedron is bounded

### Notes

The algorithm internally calls `isempty` to check whether the polyhedron is empty.
This computation can be avoided using the `check_nonempty` flag.

### Algorithm

The algorithm is based on Stiemke's theorem of alternatives, see e.g. [1].

Let the polyhedron ``P`` be given in constraint form ``Ax ≤ b``. We assume that
the polyhedron is not empty.

Proposition 1. If ``\\ker(A)≠\\{0\\}``, then ``P`` is unbounded.

Proposition 2. Assume that ``ker(A)={0}`` and ``P`` is non-empty.
Then ``P`` is bounded if and only if the following linear
program admits a feasible solution: ``\\min∥y∥_1`` subject to ``A^Ty=0`` and ``y≥1``.

[1] Mangasarian, Olvi L. *Nonlinear programming.*
    Society for Industrial and Applied Mathematics, 1994.
"""
function _isbounded_stiemke(P::HPolyhedron{N}; solver=LazySets.default_lp_solver(N),
                            check_nonempty::Bool=true) where {N}
    if check_nonempty && isempty(P)
        return true
    end

    A, b = tosimplehrep(P)
    m, n = size(A)

    if !isempty(nullspace(A))
        return false
    end

    At = copy(transpose(A))
    c = ones(N, m)
    lp = linprog(c, At, '=', zeros(n), one(N), Inf, solver)
    return is_lp_optimal(lp.status)
end

function is_hyperplanar(P::HPolyhedron)
    clist = P.constraints
    m = length(clist)

    # check that the number of constraints is fine
    if m > 2
        # try to remove redundant constraints
        clist = remove_redundant_constraints(clist)
        m = length(clist)
    end
    if m != 2
        return false
    end

    # check that the two half-spaces are complementary
    return @inbounds iscomplement(clist[1], clist[2])
end

# ============================================
# Functionality that requires Symbolics
# ============================================
function load_symbolics_hpolyhedron()

return quote

"""
    HPolyhedron(expr::Vector{<:Num}, vars=_get_variables(expr); N::Type{<:Real}=Float64)

Return the polyhedron in half-space representation given by a list of symbolic expressions.

### Input

- `expr` -- vector of symbolic expressions that describes each half-space
- `vars` -- (optional, default: `_get_variables(expr)`), if an array of variables is given,
            use those as the ambient variables in the set with respect to which derivations
            take place; otherwise, use only the variables which appear in the given
            expression (but be careful because the order may be incorrect; it is advicsed to always
            pass `vars` explicitly)
- `N`    -- (optional, default: `Float64`) the numeric type of the returned half-space

### Output

An `HPolyhedron`.

### Examples

```jldoctest
julia> using Symbolics

julia> vars = @variables x y
2-element Vector{Num}:
 x
 y

julia> HPolyhedron([x + y <= 1, x + y >= -1], vars)
HPolyhedron{Float64, Vector{Float64}}(HalfSpace{Float64, Vector{Float64}}[HalfSpace{Float64, Vector{Float64}}([1.0, 1.0], 1.0), HalfSpace{Float64, Vector{Float64}}([-1.0, -1.0], 1.0)])

julia> X = HPolyhedron([x == 0, y <= 0], vars)
HPolyhedron{Float64, Vector{Float64}}(HalfSpace{Float64, Vector{Float64}}[HalfSpace{Float64, Vector{Float64}}([1.0, 0.0], -0.0), HalfSpace{Float64, Vector{Float64}}([-1.0, -0.0], 0.0), HalfSpace{Float64, Vector{Float64}}([0.0, 1.0], -0.0)])
```
"""
function HPolyhedron(expr::Vector{<:Num}, vars::AbstractVector{Num}; N::Type{<:Real}=Float64)
    clist = Vector{HalfSpace{N, Vector{N}}}()
    sizehint!(clist, length(expr))
    got_hyperplane = false
    got_halfspace = false
    zeroed_vars = Dict(v => zero(N) for v in vars)
    vars_list = collect(vars)
    for ex in expr
        exval = Symbolics.value(ex)
        got_hyperplane, sexpr = _is_hyperplane(exval)
        if !got_hyperplane
            got_halfspace, sexpr = _is_halfspace(exval)
            if !got_halfspace
                throw(ArgumentError("expected an expression describing either " *
                    "a half-space of a hyperplane, got $expr"))
            end
        end

        coeffs = [N(α.val) for α in gradient(sexpr, vars_list)]
        β = -N(Symbolics.substitute(sexpr, zeroed_vars))

        push!(clist, HalfSpace(coeffs, β))
        if got_hyperplane
            push!(clist, HalfSpace(-coeffs, -β))
        end
    end
    return HPolyhedron(clist)
end

HPolyhedron(expr::Vector{<:Num}; N::Type{<:Real}=Float64) = HPolyhedron(expr, _get_variables(expr); N=N)
HPolyhedron(expr::Vector{<:Num}, vars; N::Type{<:Real}=Float64) = HPolyhedron(expr, _vec(vars); N=N)

end end  # quote / load_symbolics_hpolyhedron()


text:
@testset "Utils" begin

txt = "Sourav Ganguly was born on 8 July 1972 in Calcutta, and is the youngest son of Chandidas and Nirupa Ganguly. Chandidas ran a flourishing print business and was one of the richest men in the city. Ganguly had a luxurious childhood and was nicknamed the 'Maharaja', meaning the 'Great King'. Ganguly's father Chandidas Ganguly died at the age of 73 on 21 February 2013 after a long illness."
sd = StringDocument(lowercase(txt))
prepare!(sd, strip_punctuation)
prepare!(sd, strip_stopwords)
sd = StringDocument(join([i for i in split(text(sd), " ") if i!= ""], " "))
sd = sort(unique(split(text(sd), " ")))
wordlist = sd
convert(Array{String,1}, wordlist)
x= word_cooccurrence_matrix(txt, 5, true)
for i in range(1, length(wordlist))
    @test x[i, i] >= 1
end

for i in range(1, length(wordlist))
    for j in range(1, length(wordlist))
        print(i," ",j, "\n")
        @test x[i, j] == x[j,i]
    end
end

ext: jl


text:
module Combinations

using ACE, JuLIP, ForwardDiff, LinearAlgebra
using JuLIP.MLIPs: IPBasis
using ACE.RPI: RPIBasis

import JuLIP.Potentials: zlist, z2i, i2z, numz,
                         evaluate!, evaluate_d!,
                         alloc_temp, alloc_temp_d,
                         evaluate, evaluate_d,
                         cutoff

import JuLIP: fltype

export FitCombiPotential,
       num_features

@doc raw"""
`mutable struct FitCombiPotential`: Defines a site potential of the form
```math
\begin{aligned}
   E_i &= F(\rho_1, \dots, \rho_P), \\
   \rho_p &= \sum_k \theta_{pk} B_k
\end{aligned}
where $$B_k$$ are the basis functions and $$\theta_{pk}$$ the parameters.
```

### Fields
* `basis` : any IPBasis is allowed for now
* `params` : the parameters, a matrix; cf above.
* `F` : the nonlinear outer function.

It is assumed that the output of F is also `::T`
"""
mutable struct FitCombiPotential{T, TP, TF, NZ} <: SitePotential
   basis::RPIBasis{T}
   params::NTuple{NZ, Matrix{TP}}
   F::TF
end

zlist(V::FitCombiPotential) = zlist(V.basis)
# z2i(V::FitCombiPotential, z) = z2i(V.basis, z)
# i2z(V::FitCombiPotential, i) = i2z(V.basis, i)
# numz(V::FitCombiPotential) = length(zlist(V.basis))

fltype(V::FitCombiPotential{T, TP}) where {T, TP} = typeof( one(T) * one(TP) )

cutoff(V::FitCombiPotential) = cutoff(V.basis)

function FitCombiPotential(basis, F, numfeatures)
   T = fltype(basis)
   NZ = length(zlist(basis))
   params = ntuple( iz -> zeros(T, (numfeatures, length(basis, iz))), NZ )
   return FitCombiPotential(basis, params, F)
end

# -------------- Parameter management ...

num_features(V::FitCombiPotential) = size(V.params[1], 1)

get_params(V::FitCombiPotential) = vcat(vec.(V.params)...)

function set_params!(V::FitCombiPotential, θ)
   idx = 0
   len = length(V.params[1])
   for iz = 1:numz(V)
      V.params[iz][:] .= θ[(idx+1):(idx+len)]
      idx += len
   end
   return V
end

function set_params(V::FitCombiPotential, p::Vector{TP}) where {TP}
   new_params = ntuple( iz -> zeros(TP, (num_features(V), length(V.basis, iz))), numz(V) )
   Vnew = FitCombiPotential(V.basis, new_params, V.F)
   idx = 0
   len = length(Vnew.params[1])
   for iz = 1:numz(V)
      Vnew.params[iz][:] .= p[(idx+1):(idx+len)]
      idx += len
   end
   return Vnew
end

# -------------- Evaluation Code

function evaluate!(tmp, V::FitCombiPotential,
                   Rs::AbstractVector{JVec{T}},
                   Zs::AbstractVector{<:AtomicNumber},
                   z0::AtomicNumber) where {T}
   iz0 = z2i(V, z0)
   # Evaluate the inner basis
   #        TODO: replace with non-allocating version!
   B = evaluate(V.basis, Rs, Zs, z0)
   # only  get the part of the basis vector that matters to this site
   B_z0 = @view B[V.basis.Bz0inds[iz0]]
   # apply the parameters to get the densities
   ρ = V.params[iz0] * B_z0
   # return the nonlinear function of the densities.
   return V.F(ρ)
end

alloc_temp_d(V::FitCombiPotential, N::Integer) =
      (dV = zeros(JVec{fltype(V)}, N),
        R = zeros(JVecF, N),
        Z = zeros(AtomicNumber, N), )

function evaluate_d!(dEs, tmpd, V::FitCombiPotential,
                   Rs::AbstractVector{JVec{T}},
                   Zs::AbstractVector{<:AtomicNumber},
                   z0::AtomicNumber) where {T}
   @show eltype(dEs), fltype(V)
   iz0 = z2i(V, z0)
   # Evaluate the densities (cf evaluate! for comments)
   B = evaluate(V.basis, Rs, Zs, z0)
   B_z0 = @view B[V.basis.Bz0inds[iz0]]
   ρ = V.params[iz0] * B_z0
   # derivative of F ---  TODO : get rid of enforcing AD here!
   dF = ForwardDiff.gradient(V.F, ρ)
   # get the collapsed coefficients
   #    ∂_{rj} Ei = ∑_p ∂_{ρp} F * ∂_{rj} ρp
   #              = ∑_k [ ∑_p ∂_{ρp} F * θ_{pk} ] * ∂_{rj} B_k
   c = dF' * V.params[iz0]   # TODO: this is ridiculous ...
   dB = evaluate_d(V.basis, Rs, Zs, z0)
   dB_z0 = @view dB[V.basis.Bz0inds[iz0], :]
   dEi = (UniformScaling.(c) * dB_z0)[:]
   for n = 1:length(dEi)
      dEs[n] = dEi[n]
   end
   return dEs
end

# --------------- Parameter Gradient

# SECOND STEP: fast evaluator! -> TODO later
# mutable struct CombiPotential{TF} <: AbstractCalculator
#    ace::PIPotential
#    F::TF
# end

end


text:
module TimetableSolver

# include files
include("types.jl")
include("solver.jl")

# module name
export TimetableSolver
# from types.jl
export SubjectCounts, Division, Teacher, Period, Timetable, Schedule
# from solver.jl
export VariableData, get_model, get_solution, convertsolution, applysolution!, solve!

end


text:
using GeometricAlgebra
using Test

@basis "++++-"

const n = 1.0v4 + 1.0v5
const n̄ = 1.0v4 - 1.0v5

point(x) = 0.5 * (2x + magnitude2(x) * n - n̄)
point_pair(A, B) = A ∧ B
circle(A, B, C) = point_pair(A, B) ∧ C
line(A, B) = circle(A, B, n)
sphere(A, B, C, D) = circle(A, B, C) ∧ D
plane(A, B, C) = sphere(A, B, C, n)
radius2(X::Trivector) = - X^2 / (X ∧ n)^2
radius2(X::Quadvector) = X^2 / (X ∧ n)^2
norm(X) = - X ⋅ n
normalize(X) = X / norm(X)
center(X) = X * n * X

@testset "3D Conformal Geometric Algebra" begin
    @testset "n & n̄" begin
        @test magnitude2(n) == 0
        @test magnitude2(n̄) == 0
        @test n ⋅ n̄ == n̄ ⋅ n == 2v
        @test n̄ * n * n̄ == 4n̄
        @test n * n̄ * n == 4n
    end

    O = point(0v)

    @testset "Conformal points" begin
        @test O == -0.5n̄

        a = sum(randn(3) .* [v1, v2, v3])
        b = sum(randn(3) .* [v1, v2, v3])

        A = point(a)
        B = point(b)

        @test A ⋅ B ≈ -0.5 * (a - b)^2
        @test all(x -> -x ⋅ n ≈ 1v, [A, B])
    end

    @testset "Circles" begin
        C1 = circle(point(1v1), point(1v2), point(-1v1))
        @test C1 == 2.0v125
        @test normalize(center(C1)) == O
        @test radius2(C1) == 1v

        C2 = circle(point(4v1), point(2v1 + 2v2), point(0v1))
        @test C2 == 4(-v124 + v125 - 4v245)
        @test normalize(center(C2)) == point(2v1)
        @test radius2(C2) == 4v
    end

    @testset "Spheres" begin
        S1 = sphere(point(1v1), point(1v2), point(1v3), point(-1v1))
        @test S1 == 2.0v1235
        @test normalize(center(S1)) == O
        @test radius2(S1) == 1v

        S2 = sphere(point(0v1), point(4v1), point(2v1 + 2v2), point(2v1 + 2v3))
        @test normalize(center(S2)) == point(2v1)
        @test radius2(S2) == 4v
    end
end

# @test C1 ∨ C2 == point(1v1)
# @test line(point(0v1), point(1v1 + 1v2)) ∨ line(point(1v1), point(1v1 + 1v2)) == point(1v1 + 1v2)


text:
using FractalTools
using Makie
using Makie
import Makie.GeometryBasics: TriangleFace

# Construct delaunay
p1, p2, p3 = [-1., √3.], [-1., -√3], [2., 0.]
tridln = TriDelaunay(p1, p2, p3, addboundarypoints=true)

while npoints(tridln) ≤ 1000
    addpoint!(tridln)
end

finegrain!(tridln, 100)

# Define function
f(x, y) = x^2 + y^2

# Plot the result
fig, ax, plt = tridelaunayplotf(tridln, f, vmarkersize=10)
wireframe!(tomesh(tridln.delaunay), linewidth=3)
display(fig)


text:
using Test
using Pkg.PlatformEngines, Pkg.BinaryPlatforms
using SHA
using BinaryBuilderBase
using CodecZlib

@testset "Prefix" begin
    mktempdir() do temp_dir
        prefix = Prefix(temp_dir)

        # Test that it's taking the absolute path
        @test prefix.path == abspath(temp_dir)

        # Test that `bindir()`, `libdirs()` and `includedir()` all work
        for dir in unique([bindir(prefix), libdirs(prefix)..., includedir(prefix)])
            @test !isdir(dir)
            mkpath(dir)
        end

        # Create a little script within the bindir to ensure we can run it
        ppt_path = joinpath(bindir(prefix), "prefix_path_test.sh")
        open(ppt_path, "w") do f
            write(f, "#!/bin/sh\n")
            write(f, "echo yolo\n")
        end
        chmod(ppt_path, 0o775)

        # These tests are taken from BinaryProvider, but here we don't have
        # withenv.  Do we have anything similar?
        #
        # # Test that our `withenv()` stuff works.  :D
        # withenv(prefix) do
        #     @test startswith(ENV["PATH"], bindir(prefix))

        #     if !Sys.iswindows()
        #         envname = Sys.isapple() ? "DYLD_FALLBACK_LIBRARY_PATH" : "LD_LIBRARY_PATH"
        #         @test startswith(ENV[envname], last(libdirs(prefix)))
        #         private_libdir = abspath(joinpath(Sys.BINDIR, Base.PRIVATE_LIBDIR))
        #         @test endswith(ENV[envname], private_libdir)

        #         # Test we can run the script we dropped within this prefix.
        #         # Once again, something about Windows | busybox | Julia won't
        #         # pick this up even though the path clearly points to the file.
        #         @test success(`$sh $(ppt_path)`)
        #         @test success(`$sh -c prefix_path_test.sh`)
        #     end
        # end

        # Test that we can control libdirs() via platform arguments
        @test libdirs(prefix, Linux(:x86_64)) == [joinpath(prefix, "lib64"), joinpath(prefix, "lib")]
        @test libdirs(prefix, Linux(:i686)) == [joinpath(prefix, "lib")]
        @test libdirs(prefix, Windows(:i686)) == [joinpath(prefix, "bin")]
        @test bindir(prefix) == joinpath(prefix, "bin")
        @test includedir(prefix) == joinpath(prefix, "include")
        @test logdir(prefix) == joinpath(prefix, "logs")
    end
end

@testset "Packaging" begin
    # Clear out previous build products
    for f in readdir(".")
        if !endswith(f, ".tar.gz") && !endswith(f, ".sha256")
            continue
        end
        Base.rm(f; force=true)
    end

    # Gotta set this guy up beforehand
    tarball_path = nothing
    tarball_hash = nothing

    temp_prefix() do prefix
        # Create random files
        mkpath(bindir(prefix))
        mkpath(last(libdirs(prefix)))
        mkpath(joinpath(prefix, "etc"))
        bar_path = joinpath(bindir(prefix), "bar.sh")
        open(bar_path, "w") do f
            write(f, "#!/bin/sh\n")
            write(f, "echo yolo\n")
        end
        baz_path = joinpath(last(libdirs(prefix)), "baz.so")
        open(baz_path, "w") do f
            write(f, "this is not an actual .so\n")
        end

        qux_path = joinpath(prefix, "etc", "qux.conf")
        open(qux_path, "w") do f
            write(f, "use_julia=true\n")
        end

        # Next, package it up as a .tar.gz file
        tarball_path, tarball_hash = @test_logs (:info, r"^Tree hash of contents of") (:info, r"^SHA256 of") package(prefix, "./libfoo", v"1.0.0"; verbose=true)
        @test isfile(tarball_path)

        # Check that we are calculating the hash properly
        tarball_hash_check = open(tarball_path, "r") do f
            bytes2hex(sha256(f))
        end
        @test tarball_hash_check == tarball_hash

        # Test that packaging into a file that already exists fails
        @test_throws ErrorException package(prefix, "./libfoo", v"1.0.0")
    end

    # Test that we can inspect the contents of the tarball
    contents = list_tarball_files(tarball_path)
    libdir_name = Sys.iswindows() ? "bin" : "lib"
    @test joinpath("bin", "bar.sh") in contents
    @test joinpath(libdir_name, "baz.so") in contents
    @test joinpath("etc", "qux.conf") in contents

    # Cleanup after ourselves
    Base.rm(tarball_path; force=true)
    Base.rm("$(tarball_path).sha256"; force=true)
end

@testset "Tree symlinking" begin
    # Make sure symlink_tree works well with symlinks
    mktempdir() do tmpdir
        # Create fake source directory
        srcdir = joinpath(tmpdir, "src")
        mkdir(srcdir)

        mkdir(joinpath(srcdir, "dir"))
        open(joinpath(srcdir, "dir", "fileA"), "w") do io
            println(io, "fileA")
        end
        open(joinpath(srcdir, "dir", "fileB"), "w") do io
            println(io, "fileB")
        end
        symlink(joinpath("dir", "fileA"), joinpath(srcdir, "sym_fileA"))
        symlink("dir", joinpath(srcdir, "sym_dir"))

        dstdir = joinpath(tmpdir, "dst")

        # Set up a symlink tree inside of dstdir
        BinaryBuilderBase.symlink_tree(srcdir, dstdir)

        @test isdir(dstdir)
        @test isdir(joinpath(dstdir, "dir"))
        @test islink(joinpath(dstdir, "sym_dir"))
        @test islink(joinpath(dstdir, "sym_fileA"))
        @test islink(joinpath(dstdir, "dir", "fileA"))
        @test islink(joinpath(dstdir, "dir", "fileB"))

        @test readlink(joinpath(dstdir, "sym_dir")) == "dir"
        @test readlink(joinpath(dstdir, "sym_fileA")) == joinpath("dir", "fileA")

        @test String(read(joinpath(dstdir, "dir", "fileA"))) == "fileA\n"
        @test String(read(joinpath(dstdir, "dir", "fileB"))) == "fileB\n"
        @test String(read(joinpath(dstdir, "sym_fileA"))) == "fileA\n"
        @test String(read(joinpath(dstdir, "sym_dir", "fileB"))) == "fileB\n"

        # Create some files in `dstdir`, then unsymlink and see what happens:
        open(joinpath(dstdir, "dir", "fileC"), "w") do io
            println(io, "fileC")
        end
        symlink(joinpath("dir", "fileB"), joinpath(dstdir, "sym_fileB"))
        symlink(joinpath("dir", "fileC"), joinpath(dstdir, "sym_fileC"))
        symlink("dir", joinpath(dstdir, "sym_dir2"))

        BinaryBuilderBase.unsymlink_tree(srcdir, dstdir)

        @test isdir(dstdir)
        @test isdir(joinpath(dstdir, "dir"))
        @test !islink(joinpath(dstdir, "sym_dir"))
        @test !islink(joinpath(dstdir, "sym_fileA"))
        @test !isfile(joinpath(dstdir, "dir", "fileA"))
        @test !isfile(joinpath(dstdir, "dir", "fileB"))
        @test isfile(joinpath(dstdir, "dir", "fileC"))
        @test islink(joinpath(dstdir, "sym_dir2"))
        @test islink(joinpath(dstdir, "sym_fileB"))
        @test islink(joinpath(dstdir, "sym_fileC"))

        @test String(read(joinpath(dstdir, "dir", "fileC"))) == "fileC\n"
        @test String(read(joinpath(dstdir, "sym_fileC"))) == "fileC\n"
        @test_throws Base.IOError realpath(joinpath(dstdir, "sym_fileB"))
    end
end

@testset "Compression" begin
    lorem = "Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua."
    mktempdir() do dir
        file = joinpath(dir, "file.txt")
        file_gz = file * ".gz"
        write(file, lorem)
        @test_logs (:info, r"Compressing files in") compress_dir(dir; verbose=true)
        # Check that there is only the compressed file
        @test readdir(dir) == [basename(file_gz)]
        # Decompress it
        @test String(transcode(GzipDecompressor, read(file_gz))) == lorem
    end
end


text:
# Extend rotation and transformation methods from Rotations and CoordinateTransformations

function Rotations.AngleAxis(rx, ry, rz)
    # Note: promoting with Float64 here in order to avoid strange behavior when passing AngleAxis{Int} to AffineMap()
    #       I don't think this is a big deal, since the only time a RotMatrix{Int} will be valid is for the identity matrix
    t = promote_type(typeof(rx), typeof(ry), typeof(rz), Float64)
    theta = √(rx^2+ry^2+rz^2)
    return theta == 0 ? AngleAxis(zero(t),one(t),zero(t),zero(t)) : AngleAxis(theta, rx, ry, rz)
end

function rot_to_axis(R)
    aa = AngleAxis(R)
    return aa.theta.*(aa.axis_x, aa.axis_y, aa.axis_z)
end

function affinemap_to_params(tform::AffineMap)
    return (rot_to_axis(tform.linear)..., tform.translation)
end

CoordinateTransformations.LinearMap(rx,ry,rz) = CoordinateTransformations.LinearMap(AngleAxis(rx,ry,rz))
CoordinateTransformations.AffineMap(rx,ry,rz,tx,ty,tz) = CoordinateTransformations.AffineMap(AngleAxis(rx,ry,rz), SVector(tx,ty,tz))

# There is some concern about the inferability of the functions below. Using Test.@inferred did not throw any errors

function (tform::AffineMap)(x::AbstractIsotropicGaussian)
    T = typeof(x)
    otherfields = [getfield(x,fname) for fname in fieldnames(typeof(x))][5:end] # first 3 fields must be `μ`, `σ`, `ϕ`, and `dirs`
    return T.name.wrapper(tform(x.μ), x.σ, x.ϕ, [tform.linear*dir for dir in x.dirs], otherfields...)
end

function (tform::AffineMap)(x::AbstractIsotropicGMM)
    T = typeof(x)
    otherfields = [getfield(x,fname) for fname in fieldnames(typeof(x))][2:end] # first field must be `gaussians`
    return T.name.wrapper([tform(g) for g in x.gaussians], otherfields...)
end

function (tform::AffineMap)(x::AbstractIsotropicMultiGMM)
    T = typeof(x)
    otherfields = [getfield(x,fname) for fname in fieldnames(typeof(x))][2:end] # first field must be `gmms`
    tformgmms = [tform(x.gmms[key]) for key in keys(x.gmms)]
    gmmdict = Dict{eltype(keys(x.gmms)),eltype(tformgmms)}()
    for (i,key) in enumerate(keys(x.gmms))
        push!(gmmdict, key=>tformgmms[i])
    end
    return T.name.wrapper(gmmdict, otherfields...)
end


text:
module TestInterface
    using Test
    using CBSOData4

    include("mock_getodata.jl")

    tbls = CBSOData4.get_tables()
    @test length(tbls) > 0

end


text:
using CodeTools, LNR, Media
using CodeTools: getthing, getmodule
import REPL

using Logging: with_logger, current_logger
using .Progress: JunoProgressLogger

ends_with_semicolon(x) = REPL.ends_with_semicolon(split(x,'\n',keepempty = false)[end])

LNR.cursor(data::AbstractDict) = cursor(data["row"], data["column"])

exit_on_sigint(on) = ccall(:jl_exit_on_sigint, Nothing, (Cint,), on)

function modulenames(data, pos)
  main = haskey(data, "module") ? data["module"] :
         haskey(data, "path") ? CodeTools.filemodule(data["path"]) :
         "Main"
  main == "" && (main = "Main")
  sub = CodeTools.codemodule(data["code"], pos)
  main, sub
end

function getmodule′(args...)
  m = getmodule(args...)
  return m == nothing ? Main : m
end

handle("module") do data
  main, sub = modulenames(data, cursor(data))

  mod = getmodule(main)
  smod = getmodule(mod, sub)

  return d(:main => main,
           :sub  => sub,
           :inactive => (mod==nothing),
           :subInactive => smod==nothing)
end

handle("allmodules") do
  sort!([string(m) for m in CodeTools.allchildren(Main)])
end

isselection(data) = data["start"] ≠ data["stop"]

withpath(f, path) =
  CodeTools.withpath(f, path == nothing || isuntitled(path) ? nothing : path)

const evallock = ReentrantLock()

handle("evalshow") do data
  fixjunodisplays()
  @dynamic let Media.input = Editor()
    @destruct [text, line, path, mod] = data
    mod = getmodule′(mod)

    lock(evallock)
    result = hideprompt() do
      Base.CoreLogging.with_logger(Atom.Progress.JunoProgressLogger(Base.CoreLogging.current_logger())) do
        withpath(path) do
          try
            res = include_string(mod, text, path, line)
            res ≠ nothing && display(res)
            res
          catch e
            # should hide parts of the backtrace here
            Base.display_error(stderr, e, catch_backtrace())
          end
        end
      end
    end
    unlock(evallock)

    Base.invokelatest() do
      display = Media.getdisplay(typeof(result), Media.pool(Editor()), default = Editor())
      !isa(result, EvalError) && ends_with_semicolon(text) && (result = nothing)
      display ≠ Editor() && result ≠ nothing && render(display, result)
    end

    nothing
  end
end

handle("eval") do data
  fixjunodisplays()
  @dynamic let Media.input = Editor()
    @destruct [text, line, path, mod, displaymode || "editor"] = data
    mod = getmodule′(mod)

    lock(evallock)
    result = hideprompt() do
      Base.CoreLogging.with_logger(Atom.Progress.JunoProgressLogger(Base.CoreLogging.current_logger())) do
        withpath(path) do
          @errs include_string(mod, text, path, line)
        end
      end
    end
    unlock(evallock)

    Base.invokelatest() do
      !isa(result, EvalError) && ends_with_semicolon(text) && (result = nothing)
      displayandrender(result)
    end
  end
end

handle("evalall") do data
  fixjunodisplays()
  @dynamic let Media.input = Editor()
    @destruct [setmod = :module || nothing, path || "untitled", code] = data
    mod = if setmod ≠ nothing
       getmodule′(setmod)
    elseif isabspath(path)
      getmodule′(CodeTools.filemodule(path))
    else
      Main
    end

    lock(evallock)
    hideprompt() do
      Base.CoreLogging.with_logger(Atom.Progress.JunoProgressLogger(Base.CoreLogging.current_logger())) do
        withpath(path) do
          result = nothing
          try
            result = include_string(mod, code, path)
          catch e
            bt = catch_backtrace()
            ee = EvalError(e, stacktrace(bt))
            if isREPL()
              printstyled(stderr, "ERROR: ", color=:red)
              Base.showerror(stderr, e, bt)
              println(stderr)
            else
              render(Console(), ee)
            end
            @msg error(d(:msg => "Error evaluating $(basename(path))",
                         :detail => string(ee),
                         :dismissable => true))
          end
          Base.invokelatest() do
            displayandrender(result)
          end
        end
      end
    end
    unlock(evallock)
  end
  return
end

handle("evalrepl") do data
  fixjunodisplays()
  @dynamic let Media.input = Console()
    @destruct [mode || nothing, code, mod || "Main"] = data
    if mode == "shell"
      code = "Base.repl_cmd(`$code`, STDOUT)"
    elseif mode == "help"
      render′(@errs getdocs(mod, code))
      return
    end
    mod = getmodule′(mod)
    if isdebugging()
      render(Console(), @errs Debugger.interpret(code))
    else
      try
        lock(evallock)
        withpath(nothing) do
          result = @errs Core.eval(mod, :(ans = include_string($mod, $code, "console")))
          !isa(result,EvalError) && ends_with_semicolon(code) && (result = nothing)
          Base.invokelatest(render′, result)
        end
        unlock(evallock)
      catch e
        showerror(stderr, e, catch_stacktrace())
      end
    end
  end
  return
end

handle("docs") do data
  @destruct [mod || "Main", word] = data
  docstring = @errs getdocs(mod, word)

  docstring isa EvalError && return Dict(:error => true)

  mtable = try getmethods(mod, word)
    catch e
      []
    end

  Dict(:error    => false,
       :type     => :dom,
       :tag      => :div,
       :contents =>  map(x -> render(Inline(), x), [docstring; mtable]))
end

handle("methods") do data
  @destruct [mod || "Main", word] = data
  mtable = @errs getmethods(mod, word)
  if mtable isa EvalError
    Dict(:error => true, :items => sprint(showerror, mtable.err))
  else
    Dict(:items => [gotoitem(m) for m in mtable])
  end
end

function getmethods(mod, word)
  methods(CodeTools.getthing(getmodule′(mod), word))
end

function getdocs(mod, word)
  md = if Symbol(word) in keys(Docs.keywords)
    Core.eval(Main, :(@doc($(Symbol(word)))))
  else
    include_string(getmodule′(mod), "@doc $word")
  end
  return md_hlines(md)
end

function gotoitem(m::Method)
  _, link = view(m)
  sig = sprint(show, m)
  sig = replace(sig, r" in .* at .*$" => "")
  Dict(:text => sig,
       :file => link.file,
       :line => link.line - 1,
       :secondary => join(link.contents))
end


text:
using MeasureIR
using Unitful: s, kHz
using SampledSignals: SampleBuf, samplerate
using Suppressor
using DSP
using Test

"""
Simulate running the given stimulus through a real system, including a linear
IR, nonlinear distortion, noise, and clock skew.

If `dist` is given, it should be a multiplier that scales the stimulus before
passing into `tanh`. The distorted signal is scaled to have the same energy as
the original. A multiplier of 1 gives pretty mild distortion, and it goes up
from there.

Returns the simulated response, as well as the IR used.
"""
function simIR(stim; ir=nothing, dist=nothing, noisepow=-Inf, skew=1.0)
    if dist !== nothing
        stimenergy = sum(x->x^2, stim)
        stim .= tanh.(stim)
        distenergy = sum(x->x^2, stim)

        stim .*= sqrt(sigenergy/distenergy)
    end
    if ir === nothing
        ir = [zeros(100); 1; randn.() .* exp.(range(-1.5, -8, length=7500))]
    end
    resp = conv(stim, ir)
    if noisepow > -Inf
        resp .+= db2amp(noisepow) .* randn.()
    end
    if skew != 1.0
        resp = resample(resp, skew)
    end

    resp, ir
end

"""
Generates a broadband random signal of length `N` that's periodic with period
`P`. `P` does not need to be an integer.
"""
function randperiodic(N, P)
    f = 1/P # fundamental freq
    # we add eps() here so we don't get hit by floating-point error when f
    # divides evenly into 1
    H = Int(fld(0.5+eps(), f)) # number of harmonics we can fit below nyquist
    sum(cos.(2π.*(h.*f.*(0:N-1).+rand()))./H for h in 1:H)
end

function testmeasure(measfunc)
    meas = measfunc(8192)
    stim = stimulus(meas)
    funcname = split(string(measfunc), ".")[2]
    @testset "$funcname - pure impulse IR" begin
        ir = analyze(meas, stim)
        maxval, maxidx = findmax(ir)
        @test maxidx == 1
        @test maximum(ir[10:end].^2) < 0.01*maxval^2
    end
    @testset "$funcname - pure delay IR" begin
        delay = 150
        ir = analyze(meas, [zeros(delay); stim])
        maxval, maxidx = findmax(ir)
        @test maxidx == delay+1
        @test maximum(ir[1:delay-10].^2) < 0.01*maxval^2
        @test maximum(ir[delay+10:end].^2) < 0.01*maxval^2
    end
    @testset "$funcname - convolved IR" begin
        delay = 250
        L = 250
        testir = [zeros(delay); randn(L) .* exp.(range(0,-8, length=L))]
        # filter the impulse response so it's within the passband of the sweep.
        # We should be able to pretty much perfectly recover it
        testir = filt(digitalfilter(Bandpass(0.05/π, 0.5),
                                    FIRWindow(hanning(128))),
                      testir)
        resp = conv(stim, testir)
        ir = analyze(meas, resp)
        @test sum((ir[1:delay+L]-testir).^2)/sum(testir.^2) < 0.01
    end

    # FIXME: is this still supported at all?
    # @testset "$funcname - noisefloor" begin
    #     @test isapprox(noisefloor(meas, stim)[], 0.0, atol=0.01)
    #     nf = noisefloor(meas, stim .+ randn.() .* 0.2)[]
    #     @test isapprox(nf,  0.2^2, rtol=0.05)
    # end

    @testset "$funcname - multichannel response" begin
        delays = [150,200]
        resp1 = [zeros(delays[1]); stim[1:end-delays[1]]]
        resp2 = [zeros(delays[2]); stim[1:end-delays[2]]]
        resp = [resp1 resp2]

        ir = analyze(meas, resp)
        for ch in 1:2
            maxval, maxidx = findmax(ir[:,ch])
            @test maxidx == delays[ch]+1
            @test maximum(ir[1:delays[ch]-10, ch].^2) < 0.01*maxval^2
            @test maximum(ir[delays[ch]+10:end, ch].^2) < 0.01*maxval^2
        end
    end

    @testset "$funcname - vector/matrix type handling" begin
        meas = measfunc(16)
        stim = stimulus(meas)
        @test stim isa Vector
        @test analyze(meas, stim) isa Vector
        @test analyze(meas, [stim stim]) isa Matrix
    end

    # FIXME: this seems to be having a problem in `similar`
    # @testset "$funcname - analysis works with SampleBuf" begin
    #     meas = measfunc(16)
    #     stim = stimulus(meas)
    #     stimbuf = SampleBuf(stim, 44100)
    #     irbuf = analyze(meas, stimbuf)
    #     @test irbuf isa SampleBuf
    #     @test samplerate(irbuf) == 44100
    # end
end

@testset "MeasureIR" begin
# don't indent contents here, so we can run tests individually in Juno more
# easily

@testset "Standard Measurements" begin
    for f in [expsweep, golay, mls, rpms, impulse]
        testmeasure(f)
    end
end

@testset "Noncausal exponential sweep" begin
    L = 4096
    meas = expsweep(L)
    stim = stimulus(meas)
    @test indmax(analyze(meas, stim, warn=false)) == 1
    @test indmax(analyze(meas, stim, noncausal=true)) == L+1
end

# only the golay measurement truncates and warns - for the others theres
# no reason to truncate - though this may be reworked later
@testset "golay - truncation warning" begin
    L = 4096
    meas = golay(L)
    ir = randn(L÷2) .* exp.(linspace(0,-8, L÷2))
    stim = stimulus(meas)
    out = @capture_err analyze(meas, conv(stim, ir))
    @test out == ""
    ir = randn(3L÷2) .* exp.(linspace(0,-8, 3L÷2))
    out = @capture_err analyze(meas, conv(stim, ir))
    @test out == string("Warning: nonsilent samples past end of analysis ",
                        "window. Check your test signal is long enough for ",
                        "the response\n")
end

# not working right now. See source for `analyze(::ExpSweep, ...) for details
# @testset "Nonlinearity detection" begin
#     L = 4096
#     meas = expsweep(L)
#     stim = stimulus(meas)
#     out = @capture_err analyze(meas, stim)
#     @test out == ""
#     out = @capture_err analyze(meas, tanh.(stim*2))
#     @test out == "Warning: Energy in noncausal IR is above noise floor. Check for nonlinearity\n"
#
#     # now try with some noise and a more realistic IR
#     testir = randn(512) .* exp.(linspace(0,-8, 512)) / 2
#     resp = conv(stim, testir)[1:length(stim)]
#     out = @capture_err analyze(meas, stim.+randn.().*0.1)
#     @test out == ""
#     out = @capture_err analyze(meas, tanh.(stim*2).+randn.().*0.1)
#     @test out == "Warning: Energy in noncausal IR is above noise floor. Check for nonlinearity\n"
# end

@testset "Golay Sequences" begin
    for L in (8, 16)
        seq = golay(L)
        @test length(seq.A) == L
        @test length(seq.B) == L
        @test stimulus(seq) == [zeros(L);
                                seq.gain .* seq.A;
                                zeros(L);
                                seq.gain .* seq.B;
                                zeros(L)]
    end

    # should upgrade non-powers-of-two
    @test golay(13) == golay(16)

    # test equality
    @test golay(16) == golay(16)
    @test golay(16) != golay(8)

    @test maximum(stimulus(golay(16, gain=1.0))) ≈ 1.0

    # can set the bandwidth
    bwlim = golay(1024, upsample=2)
    spec = rfft(stimulus(bwlim))
    @test maximum(abs.(spec[round(Int, 0.6*length(spec)):end])) < 0.02
    resp = (analyze(bwlim, stimulus(bwlim)))
    # response energy should be compact into a bandlimited impulse
    @test sum(resp[20:60].^2) > 1.95
    @test maximum(resp[80:end]) < 1e16
end

@testset "Unitful Golay initialization" begin
    # unitful sampling rate is optional
    @test golay(2s, 1kHz) == golay(2000)
    @test golay(2s, 1000) == golay(2000)
end

@testset "Unitful Impulse" begin
    @test impulse(0.5s, 48_000) == impulse(24_000)
    @test impulse(0.5s, 48kHz) == impulse(24_000)
end

@testset "snr" begin
    sig = randn(100000)
    noise = randn(100000)

    @test isapprox(snr(sig*2+noise, sig), 4; rtol=0.02)
end

# sr = 48kHz
# f = 1kHz
# dur = 1s
# sig = [zeros(50000); pilot(dur, f, samplerate=sr); zeros(48000*10)]
# sig .+= sqrt(2) .* randn.()
# plt = plot(welch_pgram(sig[(1:10000) .+ 50000], 256).power)
@testset "findpilot" begin
    sr = 48kHz
    dur = 1s
    @testset "f=$f" for f in 0.99kHz:0.005kHz:1.01kHz
        sig = [zeros(50000); pilot(dur, f, samplerate=sr); zeros(48000*10)]
        sig .+= 8 .* randn.() # SNR of -21.1dB - not bad!
        onset = findpilot(sig, 1kHz, 1s; samplerate=sr)
        @test isapprox(onset, 50001, rtol=0.1)
    end
end

@testset "pilotsync" begin
    dur = 5s
    sr = 48kHz
    f = 1kHz
    L = 5*48000
    @testset "skew=$skew" for skew in 0.99:0.005:1.01
        p = pilot(dur, f; samplerate=sr)
        sig = [zeros(50000);
               isapprox(skew, 1) ? p : resample(p, skew);
               zeros(48000*15)]
        sig .+= randn.() ./ sqrt(2) # SNR of 0dB
        onset = findpilot(sig, f, dur; samplerate=sr)
        @test isapprox(pilotsync(sig[(0:L-1) .+ onset], f; samplerate=sr),
                       1/skew, rtol=5e-7)
    end
end

@testset "getperiod" begin
    N = 65535
    meas = mls(N, 20)

    # test at a variety of skews in ppm
    skews = 1 .+ (-400:80:400) ./ 1e6

    @testset "skew=$skew" for skew in skews
        resp, ir = simIR(stimulus(meas), skew=skew, noisepow=-8.2+20) # -20dB SNR
        period = getperiod(resp, N)
        @test isapprox(period, skew*N, rtol=1e-7)
    end
end

end # @testset MeasureIR


text:
using ThreeBodyTB

println("Precompiling, this may take a while...")

ThreeBodyTB.my_precompile()

println("Done Precompiling.")
println()

##
###precompile all codes used in examples
##
##println("Precompiling, this may take a while...")
##
###include(joinpath( ThreeBodyTB.TESTDIR ,"test_examples.jl"))
##
##
##ThreeBodyTB.BandStruct.set_no_display(true) #do not display plots
##
##tryusing(pkgsym) = try
##    @eval using $pkgsym
##    return true
##catch e
##    return e
##end
##
##
##if tryusing(:Suppressor) #use Suppressor to avoid excessive output
##
##    using $Suppressor
##    @suppress for t in ["calculate_energy_plot_fcc_al.jl", "load_dft_plot.jl", "load_tbc_plot.jl", "plot_Sc_P.jl", "relax_fcc_al.jl"]
##        include(joinpath( ThreeBodyTB.EXAMPLESDIR ,t))
##    end
##
##else #Suppressor is missing
##
##    for t in ["calculate_energy_plot_fcc_al.jl", "load_dft_plot.jl", "load_tbc_plot.jl", "plot_Sc_P.jl", "relax_fcc_al.jl"]
##        include(joinpath( ThreeBodyTB.EXAMPLESDIR ,t))
##    end
##
##end
##
##ThreeBodyTB.BandStruct.set_no_display(false) #default behavior


text:
"""
`module  JAC.PhotoExcitationFluores`
    ... a submodel of JAC that contains all methods for computing photo-excitation-autoionization cross
        sections and rates.
"""
module PhotoExcitationFluores

    using Printf, ..AutoIonization, ..Basics, ..ManyElectron, ..Radial, ..PhotoEmission, ..TableStrings
    #x global JAC_counter = 0

    """
    `struct  PhotoExcitationFluores.Settings`
        ... defines a type for the details and parameters of computing photon-impact excitation-autoionization
            pathways |i(N)>  --> |m(N)>  --> |f(N-1)>.

        + multipoles              ::Array{EmMultipole,1}   ... Specifies the multipoles of the radiation field that are to be included.
        + gauges                  ::Array{UseGauge,1}      ... Specifies the gauges to be included into the computations.
        + printBeforeComputation  ::Bool                   ... True, if all energies and lines are printed before their evaluation.
        + selectPathways          ::Bool                   ... True if particular pathways are selected for the computations.
        + selectedPathways        ::Array{Tuple{Int64,Int64,Int64},1}  ... List of list of pathways, given by tupels (inital, inmediate, final).
    """
    struct Settings
        multipoles                ::Array{EmMultipole,1}
        gauges                    ::Array{UseGauge,1}
        printBeforeComputation    ::Bool
        selectPathways            ::Bool
        selectedPathways          ::Array{Tuple{Int64,Int64,Int64},1}
     end

    """
    `PhotoExcitationFluores.Settings()`
        ... constructor for the default values of photon-impact excitation-autoionizaton settings.
    """
    function Settings()
        Settings( Basics.EmMultipole[], UseGauge[], false,  false, Tuple{Int64,Int64,Int64}[])
    end

    # `Base.show(io::IO, settings::PhotoExcitationFluores.Settings)`
    # 	 ... prepares a proper printout of the variable settings::PhotoExcitationFluores.Settings.
    function Base.show(io::IO, settings::PhotoExcitationFluores.Settings)
        println(io, "multipoles:              $(settings.multipoles)  ")
        println(io, "gauges:                  $(settings.gauges)  ")
        println(io, "printBeforeComputation:  $(settings.printBeforeComputation)  ")
        println(io, "selectPathways:          $(settings.selectPathways)  ")
        println(io, "selectedPathways:        $(settings.selectedPathways)  ")
    end

    #==
    """
    `struct  PhotoExcitationFluores.Channel`
        ... defines a type for a photon-impact excitaton & autoionization channel that specifies
            all quantum numbers, phases and amplitudes.

        + excitationChannel  ::PhotoEmission.Channel       ... Channel that describes the photon-impact excitation process.
        + augerChannel       ::AutoIonization.Channel      ... Channel that describes the subsequent Auger/autoionization process.
    """
    struct  Channel
        excitationChannel    ::PhotoEmission.Channel
        augerChannel         ::AutoIonization.Channel
    end ==#

    """
    `struct  PhotoExcitationFluores.Pathway`  ... defines a type for a photon-impact excitation pathway that may include the definition
                                                      of different excitation and autoionization channels and their corresponding amplitudes.

        + initialLevel        ::Level                  ... initial-(state) level
        + intermediateLevel   ::Level                  ... intermediate-(state) level
        + finalLevel          ::Level                  ... final-(state) level
        + excitEnergy         ::Float64                ... photon excitation energy of this pathway
        + fluorEnergy         ::Float64                ... photon fluorescence energy of this pathway
        + crossSection        ::EmProperty             ... total cross section of this pathway
        + hasChannels         ::Bool                   ... Determines whether the individual excitation and fluorescence channels are defined
                                                           in terms of their multipole, gauge as well as the amplitude, or not.
        + excitChannels       ::Array{PhotoEmission.Channel,1}  ... List of excitation channels of this pathway.
        + fluorChannels       ::Array{PhotoEmission.Channel,1}  ... List of fluorescence channels of this pathway.
    """
    struct  Pathway
        initialLevel          ::Level
        intermediateLevel     ::Level
        finalLevel            ::Level
        excitEnergy           ::Float64
        fluorEnergy           ::Float64
        crossSection          ::EmProperty
        hasChannels           ::Bool
        excitChannels         ::Array{PhotoEmission.Channel,1}
        fluorChannels         ::Array{PhotoEmission.Channel,1}
    end

    """
    `PhotoExcitationFluores.Pathway()`
        ... 'empty' constructor for an photon-impact excitation-autoionization pathway between a specified initial,
            intermediate and final level.
    """
    function Pathway()
        Pathway(Level(), Level(), Level(), 0., 0., 0., false, PhotoExcitationFluores.Channel[] )
    end

    # `Base.show(io::IO, pathway::PhotoExcitationFluores.Pathway)`
    #		 ... prepares a proper printout of the variable pathway::PhotoExcitationFluores.Pathway.
    function Base.show(io::IO, pathway::PhotoExcitationFluores.Pathway)
        println(io, "initialLevel:               $(pathway.initialLevel)  ")
        println(io, "intermediateLevel:          $(pathway.intermediateLevel)  ")
        println(io, "finalLevel:                 $(pathway.finalLevel)  ")
        println(io, "excitEnergy                 $(pathway.excitEnergy)  ")
        println(io, "fluorEnergy                 $(pathway.fluorEnergy)  ")
        println(io, "crossSection:               $(pathway.crossSection)  ")
        println(io, "hasChannels:                $(pathway.hasChannels)  ")
        println(io, "excitChannels:              $(pathway.excitChannels)  ")
        println(io, "fluorChannels:              $(pathway.fluorChannels)  ")
    end

    """
    `PhotoExcitationFluores.computeAmplitudesProperties(pathway::PhotoExcitationFluores.Pathway, grid::Radial.Grid,
                                                            settings::PhotoExcitationFluores.Settings)`
        ... to compute all amplitudes and properties of the given line; a pathway::PhotoExcitationFluores.Pathway is
            returned for which the amplitudes and properties have now been evaluated.
    """
    function  computeAmplitudesProperties(pathway::PhotoExcitationFluores.Pathway, grid::Radial.Grid, settings::PhotoExcitationFluores.Settings)
        # Compute all excitation channels
        neweChannels = PhotoEmission.Channel[]
        for eChannel in pathway.excitChannels
            amplitude   = PhotoEmission.amplitude("absorption", eChannel.multipole, eChannel.gauge, pathway.excitEnergy,
                                                  pathway.intermediateLevel, pathway.initialLevel, grid)
             push!( neweChannels, PhotoEmission.Channel( eChannel.multipole, eChannel.gauge, amplitude))
        end
        # Compute all fluorescence channels
        newfChannels = PhotoEmission.Channel[]
        for fChannel in pathway.fluorChannels
            amplitude   = PhotoEmission.amplitude("emission", fChannel.multipole, fChannel.gauge, pathway.fluorEnergy,
                                                  pathway.finalLevel, pathway.intermediateLevel, grid)
            push!( newfChannels, PhotoEmission.Channel( fChannel.multipole, fChannel.gauge, amplitude))
        end
        crossSection = EmProperty(-1., -1.)
        pathway = PhotoExcitationFluores.Pathway( pathway.initialLevel, pathway.intermediateLevel, pathway.finalLevel, pathway.excitEnergy,
                                                  pathway.fluorEnergy, crossSection, true, neweChannels, newfChannels)
        return( pathway )
    end

    """
    `PhotoExcitationFluores.computePathways(finalMultiplet::Multiplet, intermediateMultiplet::Multiplet, initialMultiplet::Multiplet,
                                            grid::Radial.Grid, settings::PhotoExcitationFluores.Settings; output=true)`
        ... to compute the photo-excitation-fluorescence amplitudes and all properties as requested by the given settings. A list of
            lines::Array{PhotoExcitationFluores.Lines} is returned.
    """
    function  computePathways(finalMultiplet::Multiplet, intermediateMultiplet::Multiplet, initialMultiplet::Multiplet, grid::Radial.Grid,
                              settings::PhotoExcitationFluores.Settings; output=true)
        println("")
        printstyled("PhotoExcitationFluores.computePathways(): The computation of photo-excitation-fluorescence amplitudes starts now ... \n", color=:light_green)
        printstyled("-------------------------------------------------------------------------------------------------------------------- \n", color=:light_green)
        println("")
        pathways = PhotoExcitationFluores.determinePathways(finalMultiplet, intermediateMultiplet, initialMultiplet, settings)
        # Display all selected pathways before the computations start
        if  settings.printBeforeComputation    PhotoExcitationFluores.displayPathways(pathways)    end
        # Calculate all amplitudes and requested properties
        newPathways = PhotoExcitationFluores.Pathway[]
        ##x println("computePathways: NO-pathwayS = $(length(pathways)) ")
        for  pathway in pathways
            ##x println("computePathways: pathway = $pathway ")
            newPathway = PhotoExcitationFluores.computeAmplitudesProperties(pathway, grid, settings)
            push!( newPathways, newPathway)
        end
        # Print all results to screen
        PhotoExcitationFluores.displayResults(stdout, pathways, settings)
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary   PhotoExcitationFluores.displayResults(iostream, pathways, settings)     end
        #
        if    output    return( pathways )
        else            return( nothing )
        end
    end

    """
    `PhotoExcitationFluores.determinePathways(finalMultiplet::Multiplet, intermediateMultiplet::Multiplet, initialMultiplet::Multiplet,
                                              settings::PhotoExcitationFluores.Settings)`
        ... to determine a list of dielectronic-recombination pathways between the levels from the given initial-, intermediate- and
            final-state multiplets and by taking into account the particular selections and settings for this computation;
            an Array{PhotoExcitationFluores.Pathway,1} is returned. Apart from the level specification, all physical properties are
            set to zero during the initialization process.
    """
    function  determinePathways(finalMultiplet::Multiplet, intermediateMultiplet::Multiplet, initialMultiplet::Multiplet,
                                settings::PhotoExcitationFluores.Settings)
        if    settings.selectPathways    selectPathways = true;
              selectedPathways = Basics.determineSelectedPathways(settings.selectedPathways, initialMultiplet, intermediateMultiplet, finalMultiplet)
        else                             selectPathways = false
        end

        pathways = PhotoExcitationFluores.Pathway[]
        for  i = 1:length(initialMultiplet.levels)
            for  n = 1:length(intermediateMultiplet.levels)
                for  f = 1:length(finalMultiplet.levels)
                    if  selectPathways  &&  !((i,n,f) in selectedPathways )    continue   end
                    eEnergy = intermediateMultiplet.levels[n].energy - initialMultiplet.levels[i].energy
                    fEnergy = intermediateMultiplet.levels[n].energy - finalMultiplet.levels[f].energy
                    if  eEnergy < 0.   ||   fEnergy < 0.    continue    end

                    rSettings = PhotoEmission.Settings( settings.multipoles, settings.gauges, false, false, false, Tuple{Int64,Int64}[], 0., 0., 0.)
                    eChannels = PhotoEmission.determineChannels(intermediateMultiplet.levels[n], initialMultiplet.levels[i], rSettings)
                    fChannels = PhotoEmission.determineChannels(finalMultiplet.levels[f], intermediateMultiplet.levels[n], rSettings)
                    push!( pathways, PhotoExcitationFluores.Pathway(initialMultiplet.levels[i], intermediateMultiplet.levels[i], finalMultiplet.levels[f],
                                                                    eEnergy, fEnergy, EmProperty(0., 0.), true, eChannels, fChannels) )
                end
            end
        end
        return( pathways )
    end

    """
    `PhotoExcitationFluores.displayResults(stream::IO, pathways::Array{PhotoExcitationFluores.Pathway,1},
                                           settings::PhotoExcitationFluores.Settings)`
        ... to list all results, energies, cross sections, etc. of the selected lines. A neat table is printed but nothing
            is returned otherwise.
    """
    function  displayResults(stream::IO, pathways::Array{PhotoExcitationFluores.Pathway,1}, settings::PhotoExcitationFluores.Settings)
        println(stream, " ")
        println(stream, "  Partial excitation & fluorescence cross sections:")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(150))
        sa = "    ";   sb = "    "
        sa = sa * TableStrings.center(23, "Levels"; na=2);            sb = sb * TableStrings.center(23, "i  --  e  --  f"; na=2);
        sa = sa * TableStrings.center(23, "J^P symmetries"; na=0);    sb = sb * TableStrings.center(23, "i  --  e  --  f"; na=2);
        sa = sa * TableStrings.center(30, "Energies  " * TableStrings.inUnits("energy"); na=4);
        sb = sb * TableStrings.center(30, "r--i        r--f  "; na=4)
        sa = sa * TableStrings.center(30, "Multipoles"; na=4);
        sb = sb * TableStrings.center(38, "i--r        r--f  "; na=4)
        sa = sa * TableStrings.center(36, "Cou -- cross sections -- Bab";        na=2)
        sb = sb * TableStrings.center(36, TableStrings.inUnits("cross section") * "          " *
                                              TableStrings.inUnits("cross section"); na=3)
        println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(150))
        #
        for  pathway in pathways
            sa  = " ";     isym = LevelSymmetry( pathway.initialLevel.J,      pathway.initialLevel.parity)
                           msym = LevelSymmetry( pathway.intermediateLevel.J, pathway.intermediateLevel.parity)
                           fsym = LevelSymmetry( pathway.finalLevel.J,        pathway.finalLevel.parity)
            sa = sa * TableStrings.center(23, TableStrings.levels_imf(pathway.initialLevel.index, pathway.intermediateLevel.index,
                                                                              pathway.finalLevel.index); na=3)
            sa = sa * TableStrings.center(23, TableStrings.symmetries_imf(isym, msym, fsym);  na=4)
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", pathway.excitEnergy))      * "   "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", pathway.fluorEnergy))      * "   "
            #
            multipoles = EmMultipole[]
            for  ech in pathway.excitChannels
                multipoles = push!( multipoles, ech.multipole)
            end
            multipoles = unique(multipoles);   mpString = TableStrings.multipoleList(multipoles)            * "                   "
            sa = sa * TableStrings.flushleft(16, mpString[1:16];  na=2)
            #
            multipoles = EmMultipole[]
            for  fch in pathway.fluorChannels
                multipoles = push!( multipoles, fch.multipole)
            end
            multipoles = unique(multipoles);   mpString = TableStrings.multipoleList(multipoles)             * "                   "
            sa = sa * TableStrings.flushleft(16, mpString[1:16];  na=2)
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("cross section: from atomic", pathway.crossSection.Coulomb))     * "  "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("cross section: from atomic", pathway.crossSection.Babushkin))   * "  "
            println(stream, sa)
        end
        println(stream, "  ", TableStrings.hLine(150))
        #
        return( nothing )
    end

end # module


text:
using Distributions
using Dates

using .Networks

struct UncertaintiesGenerator <: AbstractDataGenerator
    network::Networks.Network
    target_timepoints::Vector{Dates.DateTime}
    horizon_timepoints::Vector{Dates.DateTime}
    #potentially split into load_uncertainties_distribution and unit_uncertainties_distribution :
    uncertainties_distribution #ditribution info for each limitable unit and load
    nb_scenarios::Int
end

function scenario_name(index_scenario::Int)
    return "S"*string(index_scenario)
end

function generate_uncertainties(network::Networks.Network,
                                target_timepoints::Vector{Dates.DateTime},
                                horizon_timepoints::Vector{Dates.DateTime},
                                uncertainties_distribution, nb_scenarios::Int)
    generator = UncertaintiesGenerator(network, target_timepoints, horizon_timepoints, uncertainties_distribution, nb_scenarios)
    return launch(generator)
end

function launch(uncertainties_generator::UncertaintiesGenerator)
    uncertainties = Uncertainties()
    if !check_uncertainties_distribution(uncertainties_generator.network, uncertainties_generator.uncertainties_distribution)
        error("Invalid uncertainties parameters!")
    end

    scenarios = scenario_name.(1:uncertainties_generator.nb_scenarios)

    for ech in uncertainties_generator.horizon_timepoints
        for bus_or_gen in union(Networks.get_buses(uncertainties_generator.network),
                                Networks.get_generators_of_type(uncertainties_generator.network, Networks.LIMITABLE))
            nodal_injection_name = Networks.get_id(bus_or_gen)
            uncertain_distro = uncertainties_generator.uncertainties_distribution[nodal_injection_name]
            for ts in uncertainties_generator.target_timepoints
                values = generate_values(uncertain_distro, ech, ts, uncertainties_generator.nb_scenarios)
                for (scenario_name,val) in zip(scenarios,values)
                    add_uncertainty!(uncertainties, ech, nodal_injection_name, ts, scenario_name, val)
                end
            end
        end
    end

    return uncertainties
end

function generate_values(uncertain_distro::UncertaintyNDistribution,
                        ech::Dates.DateTime, ts::Dates.DateTime,
                        nb_scenarios::Int64)

    delta_time_l = max(0,
                    Dates.value(floor(ts - ech, Dates.Second)) / 3600)
    factor_l = delta_time_l * uncertain_distro.time_factor
    adjusted_sigma_l = factor_l * uncertain_distro.sigma
    random_values_l = rand(Distributions.Normal(uncertain_distro.mu, adjusted_sigma_l), nb_scenarios)
    values_l = max.(uncertain_distro.min_value, random_values_l)
    values_l = min.(uncertain_distro.max_value, values_l)

    return values_l
end

############################################################
#         Checkers
############################################################

function check_extra_rows(network, uncertainties_distribution)
    #warn if extra input
    for (id, _) in uncertainties_distribution
        gen_or_bus = Networks.get_generator_or_bus(network, id)
        if ( ismissing(gen_or_bus) ||
            ( isa(gen_or_bus, Networks.Generator) && (Networks.get_type(gen_or_bus) != Networks.LIMITABLE) )
            )
            msg = @sprintf("Uncertainty distribution input for id %s may not be needed.", id)
            @warn(msg)
        end
    end

    return true
end

function check_buses_listed(network, uncertainties_distribution)
    checks = true

    #all buses are listed
    for bus in Networks.get_buses(network)
        bus_id = Networks.get_id(bus)
        if !haskey(uncertainties_distribution, bus_id)
            msg = @sprintf("Missing uncertainty distribution for bus %s.", bus_id)
            @error(msg)
            checks = false
        end
    end

    return checks
end

function check_limitables_listed(network, uncertainties_distribution)
    checks = true

    #all limitables are listed
    for gen in Networks.get_generators_of_type(network, Networks.LIMITABLE)
        gen_id = Networks.get_id(gen)
        if !haskey(uncertainties_distribution, gen_id)
            msg = @sprintf("Missing uncertainty distribution for limitable %s.", gen_id)
            @error(msg)
            checks = false
        end
    end

    return checks
end

function check_min_max(network, uncertainties_distribution)
    checks = true

    #max_value of limitables are respected
    for (id, uncertain_distro) in uncertainties_distribution
        @assert(id == uncertain_distro.id)
        gen_or_bus = Networks.get_generator_or_bus(network, id)
        if !ismissing(gen_or_bus) && isa(gen_or_bus, Networks.Generator)
            if Networks.get_p_max(gen_or_bus) < uncertain_distro.max_value
                msg = @sprintf("Invalid Uncertainty distribution max value (%f) for id %s. \
                                Value must be less or equal to %f (the generator's Pmax).",
                                uncertain_distro.max_value, id, Networks.get_p_max(gen_or_bus))
                @error(msg)
                checks = false
            end
        end
    end

    #  max_value >= min_value >= 0
    for (id, uncertain_distro) in uncertainties_distribution
        if uncertain_distro.max_value < uncertain_distro.min_value
            msg = @sprintf("Invalid Uncertainty distribution for id %s. \
                                max_value (%f) must be greater or equal to min_value (%f).",
                                id, uncertain_distro.max_value, uncertain_distro.min_value)
                @error(msg)
                checks = false
        end
        if uncertain_distro.min_value < 0
            msg = @sprintf("Invalid Uncertainty distribution for id %s. \
                                pmin (%f) must be non-negative.",
                                id, uncertain_distro.min_value)
                @error(msg)
                checks = false
        end
    end

    return checks
end

function check_uncertainties_distribution(network, uncertainties_distribution)
    return ( check_extra_rows(network, uncertainties_distribution)
            & check_buses_listed(network, uncertainties_distribution)
            & check_limitables_listed(network, uncertainties_distribution))
end

function check_uncertainties_distribution(network, uncertainties_distribution::UncertaintiesDistribution)
    return ( check_extra_rows(network, uncertainties_distribution)
            & check_buses_listed(network, uncertainties_distribution)
            & check_limitables_listed(network, uncertainties_distribution)
            & check_min_max(network, uncertainties_distribution) )
end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
macro Σein!(ex::Expr)
    ex.head == :string && isodd(length(ex.args)) || error("")
    s = join(ex.args[1:2:end])
    tensors_rhs = parse_tensors_rhs(s)
    indices_lhs = parse_indices_lhs(s)
    pairs_index = parse_pairs_index(s)
    tree_contract = parse_tree_contract(s)
    for ((kind, klab), v) in pairs_index
        replace!(tensors_rhs[klab].indices, (kind, klab) => v)
    end
    indices_all = union([x.indices for x in tensors_rhs]...)
    indices_sym = Dict(x => gensym() for x in indices_all)

    tensors_rhs_sym = Vector{Symbol}(undef, length(tensors_rhs))
    tensor_lhs_sym = Ref{Symbol}()
    for (i, (str, sym)) in enumerate(zip(ex.args[3:2:end], ex.args[2:2:end]))
        label = match(r"\s*@??\s*(?<label>[0-9>]+)", str)[:label]
        if label ≠ ">"
            tensors_rhs_sym[parse(Int, label)] = sym
        else
            tensor_lhs_sym[] = sym
        end
    end

    tensors_rhs_ex = Vector{Expr}()
    for (i, a) in enumerate(tensors_rhs)
        ex = Expr(:ref, esc(tensors_rhs_sym[i]), [indices_sym[i] for i in a.indices]...)
        ex = a.isconj ? Expr(:call, :conj, ex) : ex
        push!(tensors_rhs_ex, ex)
    end

    ex_rhs_lhs = Expr(:ref, esc(tensor_lhs_sym[]), [indices_sym[i] for i in indices_lhs]...)
    ex_rhs_rhs = if isnothing(tree_contract)
        Expr(:call, :(*), tensors_rhs_ex...)
    else
        evaltree(i, j) = Expr(
            :call,
            :(*),
            i isa Int ? tensors_rhs_ex[i] : evaltree(i...),
            j isa Int ? tensors_rhs_ex[j] : evaltree(j...),
        )
        evaltree(i, j, k...) = evaltree((i, j), k...)
        evaltree(tree_contract...)
    end

    ex_rhs =
        macroexpand(TensorOperations, :(TensorOperations.@tensor $ex_rhs_lhs = $ex_rhs_rhs))
    return :($ex_rhs)
end

macro Σein(ex::Expr)
    ex.head == :string && isodd(length(ex.args)) || error("")
    s = join(ex.args[1:2:end])
    tensors_rhs = parse_tensors_rhs(s)
    indices_lhs = parse_indices_lhs(s)
    pairs_index = parse_pairs_index(s)
    tree_contract = parse_tree_contract(s)
    for ((kind, klab), v) in pairs_index
        replace!(tensors_rhs[klab].indices, (kind, klab) => v)
    end
    indices_all = union([x.indices for x in tensors_rhs]...)
    indices_sym = Dict(x => gensym() for x in indices_all)

    tensors_rhs_sym = Vector{Symbol}(undef, length(tensors_rhs))
    tensor_lhs_sym = gensym()
    for (i, (str, sym)) in enumerate(zip(ex.args[3:2:end], ex.args[2:2:end]))
        label = match(r"\s*@??\s*(?<label>[0-9]+)", str)[:label]
        tensors_rhs_sym[parse(Int, label)] = sym
    end

    tensors_rhs_ex = Vector{Expr}()
    for (i, a) in enumerate(tensors_rhs)
        ex = Expr(:ref, esc(tensors_rhs_sym[i]), [indices_sym[i] for i in a.indices]...)
        ex = a.isconj ? Expr(:call, :conj, ex) : ex
        push!(tensors_rhs_ex, ex)
    end

    ex_rhs_lhs = Expr(:ref, tensor_lhs_sym, [indices_sym[i] for i in indices_lhs]...)
    ex_rhs_rhs = if isnothing(tree_contract)
        Expr(:call, :(*), tensors_rhs_ex...)
    else
        evaltree(i, j) = Expr(
            :call,
            :(*),
            i isa Int ? tensors_rhs_ex[i] : evaltree(i...),
            j isa Int ? tensors_rhs_ex[j] : evaltree(j...),
        )
        evaltree(i, j, k...) = evaltree((i, j), k...)
        evaltree(tree_contract...)
    end

    ex_rhs =
        macroexpand(TensorOperations, :(TensorOperations.@tensor $ex_rhs_lhs := $ex_rhs_rhs))
    return :($ex_rhs)
end


text:
module GutzwillerMPS

using MatrixProductStateTools
using SymTensors
using LinearAlgebra
using TensorOperations

include("gutzwiller.jl")

export gutzwillerexact
export zipandgutzwiller!
export _tensorproductzip!
export _applygutzwillerzip!

end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
using Convex
using ECOS

"""
    probability_state_table(x::Tuple, y::Tuple, loc::Symbol, alleles::Dict)
Calculate the probability of observing the particular allele state given each of
the 9 Jacquard Identity States for a single locus to create Table 1 from
Milligan 2003.
"""
function probability_state_table(x::Tuple, y::Tuple, alleles::Dict)
    #TODO Improve how groups are decided based on how similar things are done with moments estimators
    #=
    Calculate Pr(Li | Sj)
    If the allele identity falls into this class (L1-L9), generate the
    probabilities of it belonging to each of the different classes and
    return that array of 9 distinct probabilities
    =#

    ## class L1 -  AᵢAᵢ AᵢAᵢ ##
    if x[1] == x[2] == y[1] == y[2]
        p = alleles[x[1]]
        [p, p^2, p^2, p^3, p^2, p^3, p^2, p^3, p^4]

    ## class L2 - AᵢAᵢ AⱼAⱼ ##
    elseif (x[1] == x[2]) & (y[1] == y[2]) & (x[1] != y[1])
        p = (alleles[x[1]], alleles[y[1]])
        [0, prod(p), 0, prod(p) * p[2], 0, prod(p) * p[1], 0, 0, prod(p) * prod(p)]

    ## class L3a - AᵢAᵢ AᵢAⱼ ## - has issues because of allele order
    elseif ((x[1] == x[2] == y[1]) & (x[1] != y[2]))
        p = (alleles[x[1]], alleles[y[2]])
        [0, 0, prod(p), 2 * prod(p) * p[1], 0, 0, 0, prod(p) * p[1], 2 * prod(p) * p[1]^2]

    ## class L3b - AᵢAᵢ AⱼAᵢ ## - has issues because of allele order
    elseif ((x[1] == x[2] == y[2]) & (x[1] != y[1]))
        p = (alleles[x[1]], alleles[y[1]])
        [0, 0, prod(p), 2 * prod(p) * p[1], 0, 0, 0, prod(p) * p[1], 2 * prod(p) * p[1]^2]

    ## class L4 - AᵢAᵢ AⱼAₖ ##
    elseif (x[1] == x[2]) & (y[1] != y[2]) & (x[1] != y[1]) & (x[1] != y[2])
        p = (alleles[x[1]], alleles[y[1]], alleles[y[2]])
        [0, 0, 0, 2 * prod(p), 0, 0, 0, 0, 2 * prod(p) * p[1]]

    ## L5a - AiAj AiAi ## - has issues because of allele order
    elseif ((x[1] == y[1] == y[2]) & (x[1] != x[2]))
        p = (alleles[x[1]], alleles[x[2]])
        [0, 0, 0, 0, prod(p), 2 * prod(p) * p[1], 0, prod(p) *p[1], 2 * prod(p) * p[1]^2]

    ## L5b - AjAi AiAi ## - has issues because of allele order
    elseif (x[2] == y[1] == y[2] & (x[1] != x[2]))
        p = (alleles[x[2]], alleles[x[1]])
        [0, 0, 0, 0, prod(p), 2 * prod(p) * p[1], 0, prod(p) *p[1], 2 * prod(p) * p[1]^2]

    ## L6 - AjAk AiAi ##
    elseif (x[1] != x[2]) & (y[1] == y[2]) & (x[1] != y[1]) & (x[2] != y[1])
        p = (alleles[y[1]], alleles[x[1]], alleles[x[2]])
        [0, 0, 0, 0, 0, 2 * prod(p), 0, 0, 2 * prod(p) * p[1]]

    ## L7 - AiAj AiAj ##
    elseif (x[1] == y[1]) & (x[2] == y[2]) & (x[1] != x[2])
        p = (alleles[x[1]], alleles[x[2]])
        [0, 0, 0, 0, 0, 0, 2 * prod(p), prod(p) * sum(p), 4 * prod(p) * prod(p)]

    ## L8a - AiAj AiAk ##  - has issues because of allele order
    elseif ((x[1] == y[1]) & (x[1] != x[2]) & (y[1] != y[2]) & (x[2] != y[2]))
        p = (alleles[x[1]], alleles[x[2]], alleles[y[2]])
        [0, 0, 0, 0, 0, 0, 0, prod(p), 4 * prod(p) * p[1]]

    ## L8b - AjAi AkAi ##  - has issues because of allele order
    elseif ((x[2] == y[2]) & (x[1] != x[2]) & (y[1] != y[2]) & (x[1] != y[1]))
        p = (alleles[x[2]], alleles[x[1]], alleles[y[1]])
        [0, 0, 0, 0, 0, 0, 0, prod(p), 4 * prod(p) * p[1]]

    ## L8c - AjAi AiAk ##  - has issues because of allele order
    elseif ((x[2] == y[1]) & (x[1] != x[2]) & (y[1] != y[2]) & (x[1] != y[2]))
        p = (alleles[x[2]], alleles[x[1]], alleles[y[2]])
        [0, 0, 0, 0, 0, 0, 0, prod(p), 4 * prod(p) * p[1]]

    ## L8d - AiAj AkAi ##  - has issues because of allele order
    elseif ((x[1] == y[2]) & (x[1] != x[2]) & (y[1] != y[2]) & (x[1] != y[1]))
        p = (alleles[x[1]], alleles[x[2]], alleles[y[1]])
        [0, 0, 0, 0, 0, 0, 0, prod(p), 4 * prod(p) * p[1]]

    ## L9 - AiAj AkAl ##
    elseif (x[1] != x[2]) & (x[1] != y[1]) & (x[1] != y[2]) & (x[2] != y[1]) & (x[2] != y[2]) & (y[1] != x[2])
        p = (alleles[x[1]], alleles[x[2]], alleles[y[1]], alleles[y[2]])
        [0, 0, 0, 0, 0, 0, 0, 0, 4 * prod(p)]
    else
        [-9, -9, -9, -9, -9, -9, -9, -9, -9]
    end
end

"""
    Δ_optim(Pr_L_S::Transpose{Float64,Array{Float64,2}}, verbose::Bool)
Takes the probability of the allelic state given the identity by descent from
all available loci (either allowing for inbreeding or not) and calculated the
maximum likelihood Δ coefficients
"""
function Δ_optim(Pr_L_S::Array{Float64,2}, verbose::Bool = false)
    #Δ is what needs to be optimized
    #consist of 9 values between 0 and 1 which must also sum to 1
    #is then used to calculate relatedness

    Δ = Variable(9)
    problem = maximize(sum(log(Pr_L_S * Δ)))
    problem.constraints += 0 <= Δ[1:9]
    problem.constraints += Δ[1:9] <= 1
    problem.constraints += sum(Δ) <= 1
    problem.constraints += 0 <= sum(Δ)
    problem.constraints += 2 * (Δ[1] + 0.5 * (Δ[3] + Δ[5] + Δ[7]) + 0.25 * Δ[8]) <= 1
    problem.constraints += 0 <= 2 * (Δ[1] + 0.5 * (Δ[3] + Δ[5] + Δ[7]) + 0.25 * Δ[8])

    Convex.solve!(problem, ECOS.Optimizer(verbose = verbose, maxit=100), verbose = verbose) #maxit=100,
    #Convex.solve!(problem, ECOSSolver(verbose = verbose, maxit=100, feastol=5e-6, reltol = 1e-3, reltol_inacc = 5e-2), verbose = verbose)
    #Convex.solve!(problem, SCSSolver(verbose = verbose, max_iters = 100), verbose = verbose)

    Δ.value, problem.status
    # Should probably include some output that confirms that it did in fact
    # converge and/or use multiple random starts to confirm not a local maxima
end

#### No inbreeding assumption
function Δ_optim_noInbreeding(Pr_L_S::Array{Float64,2}, verbose::Bool = false)
    #Δ is what needs to be optimized
    #consist of 9 values between 0 and 1 which must also sum to 1
    #is then used to calculate relatedness

    Δ = Variable(3)
    problem = maximize(sum(log(Pr_L_S * Δ)))
    problem.constraints += 0 <= Δ[1:3]
    problem.constraints += Δ[1:3] <= 1
    problem.constraints += sum(Δ) <= 1
    problem.constraints += 0 <= sum(Δ)
    problem.constraints += 2 * (0.5 * Δ[1] + 0.25 * Δ[2]) <= 1
    problem.constraints += 0 <= 2 * (0.5 * Δ[1] + 0.25 * Δ[2])

    Convex.solve!(problem, ECOS.Optimizer(verbose = verbose, maxit=100), verbose = verbose) #maxit=100,
    #Convex.solve!(problem, ECOSSolver(verbose = verbose, maxit=100, feastol=5e-6, reltol = 1e-3, reltol_inacc = 5e-2), verbose = verbose)
    #Convex.solve!(problem, SCSSolver(verbose = verbose, max_iters = 100), verbose = verbose)

    Δ.value, problem.status
    # Should probably include some output that confirms that it did in fact
    # converge and/or use multiple random starts to confirm not a local maxima
end

"""
dyadicLikelihood(ind1::GenoArray, ind2::GenoArray, locus_names::Vector{Symbol}; alleles::NamedTuple)
Calculates the maximum likelihood based relatedness using all available loci following following Milligan (2002)

-Single Locus Equation:
-How to combine multiple loci: NA inherently multi-locus

-Assumes inbreeding can be present

Milligan, B. G. (2003). Maximum-likelihood estimation of relatedness. Genetics, 163(3), 1153-1167.
"""
function dyadicLikelihood(ind1::T, ind2::T, locus_names::Vector{Symbol}, alleles::U; kwargs...) where T <: GenoArray where U <: NamedTuple
    kw_dict = Dict(kwargs...)
    Pr_Ls = Array{Float64}(undef, length(locus_names), 9)

    idx = 0
    @inbounds for (loc,gen1,gen2) in zip(locus_names, ind1, ind2)
        idx += 1
        Pr_Ls[idx, :] = probability_state_table(gen1, gen2, alleles[loc])
    end

    if kw_dict[:inbreeding] == true
        Δ = Δ_optim(Pr_Ls)
        θ = Δ[1][1] + 0.5 * (Δ[1][3] + Δ[1][5] + Δ[1][7]) + 0.25 * Δ[1][8]
    else
        Pr_Ls = Pr_Ls[:, 7:9]
        Δ = Δ_optim_noInbreeding(Pr_Ls)
        θ = (0.5 * Δ[1][1] + 0.25 * Δ[1][2])
    end
    return 2 * θ#, Δ[1], Δ[2] #secondary outputs for convergence & Δvalues if we want
end

#= should be redundant now
function dyadicLikelihood_noInbreeding(ind1::T, ind2::T, locus_names::Vector{Symbol}, alleles::U; kwargs...) where T <: GenoArray where U <: NamedTuple
    #TODO Add inbreeding toggle
    Pr_Ls = Array{Float64}(undef, length(locus_names), 9)

    idx = 0
    for (loc,gen1,gen2) in zip(locus_names, ind1, ind2)
        idx += 1
        Pr_Ls[idx, :] = probability_state_table(gen1, gen2, allele_frequencies[loc])
    end
    Pr_Ls = Pr_Ls[:, 7:9]

    Δ = Δ_optim_noInbreeding(Pr_Ls)

    θ = (0.5 * Δ[1][1] + 0.25 * Δ[1][2])
    return 2 * θ#, Δ[1], Δ[2] #secondary outputs for convergence & Δvalues if we want
end
=#


text:
export Peak_Motif

using StaticArrays: @SVector, SVector

# note:
# ----------------------------------------------------------------
# as before uses a parametrized type to limit memory allocs
#
struct Peak_Motif{P<:Abstract_Model2Fit_Peak} <: Abstract_Model2Fit_Peak
    _peak::P
    _profile::Matrix{Float64} # position, height
end

# Visit  ================
#
# This is *generic*
#
visit_submodel_size(model::Peak_Motif) = 0

# Gaussian peak specialisation ****************
# θ = (h,σ)
#
parameter_size(pm::Peak_Motif{Gaussian_Peak}) = 2

function accumulate_y!(pm::Peak_Motif{Gaussian_Peak},Y::AbstractVector,X::AbstractVector,θ::AbstractVector{T}) where {T}
    @assert length(θ) == parameter_size(pm)

    h_glob = θ[1]
    σ_glob = θ[2]

    n = size(pm._profile,1)

    for i in 1:n
        μ_loc = pm._profile[i,1]
        h_loc = pm._profile[i,2]

        accumulate_y!(pm._peak, Y, X, @SVector T[ h_glob*h_loc,μ_loc,σ_glob])
    end

    Y
end

# Other peak specialisation ****************
# θ = ...
#