text:
# Tried to approximate the approach from
# https://mrzv.org/publications/geometry-helps-distances-persistence-diagrams/alenex/ and
# https://www2.cs.arizona.edu/~alon/papers/match.pdf with NearestNeighbors.jl but the
# allocations were so high, it was slower for moderately sized inputs. This approach is also
# much much simpler.
# TODO: try it again some time.

abstract type MatchingDistance end

"""
    weight(::MatchingDistance, left, right)
    weight(::Matching)

Get the weight of the matching between persistence diagrams `left` and `right`.

# See also

* [`matching`](@ref)
* [`Bottleneck`](@ref)
* [`Wasserstein`](@ref)
"""
weight(dist::MatchingDistance, left, right) = dist(left, right; matching=false)
"""
    matching(::MatchingDistance, left, right)
    matching(::Matching)

Get the matching between persistence diagrams `left` and `right`.

# See also

* [`weight`](@ref)
* [`Bottleneck`](@ref)
* [`Wasserstein`](@ref)
"""
matching(dist::MatchingDistance, left, right) = dist(left, right; matching=true)

"""
    Matching

A matching between two persistence diagrams.

# Methods

* [`weight(::Matching)`](@ref)
* [`matching(::Matching)`](@ref)
"""
struct Matching
    left::PersistenceDiagram
    right::PersistenceDiagram
    weight::Float64
    matching::Vector{Pair{Int,Int}}
    bottleneck::Bool
end

weight(match::Matching) = match.weight

Base.length(match::Matching) = length(match.matching)
Base.isempty(match::Matching) = isempty(match.matching)

function _distance(int1, int2)
    if isfinite(int1) && isfinite(int2)
        return max(abs(birth(int1) - birth(int2)), abs(death(int1) - death(int2)))
    elseif !isfinite(int1) && !isfinite(int2)
        return abs(birth(int1) - birth(int2))
    else
        return Inf
    end
end

function _distances(left, right)
    dists = zeros(length(right), length(left))
    for j in eachindex(left), i in eachindex(right)
        dists[i, j] = _distance(left[j], right[i])
    end
    return dists
end

function matching(match::Matching; bottleneck=match.bottleneck)
    result = Pair{PersistenceInterval,PersistenceInterval}[]
    n = length(match.left)
    m = length(match.right)
    for (i, j) in match.matching
        if i ≤ n && j ≤ m
            push!(result, match.left[i] => match.right[j])
        elseif i ≤ n
            # left is matched to diagonal
            l = match.left[i]
            push!(result, match.left[i] => PersistenceInterval(birth(l), birth(l)))
        elseif j ≤ m
            # right is matched to diagonal
            r = match.right[j]
            push!(result, PersistenceInterval(birth(r), birth(r)) => r)
        end
    end
    sort!(result)

    if !bottleneck
        return result
    else
        return filter!(m -> _distance(m...) == match.weight, result)
    end
end

function Base.summary(io::IO, match::Matching)
    b = match.bottleneck ? "Bottleneck " : ""
    return print(io, "$(b)Matching with weight $(match.weight)")
end
function Base.show(io::IO, match::Matching)
    return Base.summary(io, match)
end
function Base.show(io::IO, ::MIME"text/plain", match::Matching)
    print(io, match)
    if length(match) > 0
        print(io, ":")
        pairs = matching(match)
        for p in pairs
            print(io, "\n ", p)
        end
    end
end

"""
    _adjacency_matrix(left::PersistenceDiagram, right::PersistenceDiagram, power)

Get the adjacency matrix of the matching between `left` and `right`. Edge weights are equal
to distances between intervals raised to the power of `power`. Distances between diagonal
points and values that should not be matched with them are set to `Inf`. The same
holds for distances between finite and infinite intervals.

For `length(left) == n` and `length(right) == m`, it returns a ``(n m) × (m n)`` matrix.

# Example

```jldoctest
julia> left = PersistenceDiagram([(0.0, 1.0), (3.0, 4.5)]);

julia> right = PersistenceDiagram([(0.0, 1.0), (4.0, 5.0), (4.0, 7.0)]);

julia> PersistenceDiagrams._adjacency_matrix(left, right)
5×5 Array{Float64,2}:
  0.0   3.5   1.0  Inf   Inf
  4.0   1.0  Inf    1.0  Inf
  6.0   2.5  Inf   Inf    3.0
  1.0  Inf    0.0   0.0   0.0
 Inf    1.5   0.0   0.0   0.0
```
"""
function _adjacency_matrix(left, right, power=1)
    left = sort(left; by=death)
    right = sort(right; by=death)

    n = length(left)
    m = length(right)
    adj = fill(Inf, n + m, m + n)

    dists = _distances(left, right)
    adj[axes(dists)...] .= dists
    for i in (size(dists, 2) + 1):n, j in (size(dists, 1) + 1):m
        adj[j, i] = abs(birth(left[i]) - birth(right[j]))
    end
    for i in 1:n
        adj[i + m, i] = persistence(left[i])
    end
    for j in 1:m
        adj[j, j + n] = persistence(right[j])
    end
    adj[(m + 1):(m + n), (n + 1):(n + m)] .= 0.0

    if power ≠ 1
        return adj .^ power
    else
        return adj
    end
end

"""
    BottleneckGraph

Representation of the bipartite graph used for computing bottleneck distance via the
Hopcroft-Karp algorithm. In all the following functions, `left` and `right` refer to the
vertex sets of the graph. The graph has `n + m` vertices in each set corresponding to the
numbers of points in the diagrams plus the diagonals.

# Fields

* `adj::Matrix{Float64}`: the adjacency matrix.
* `match_left::Vector{Int}`: matches of left vertices.
* `match_right::Vector{Int}`: matches of right vertices.
* `edges::Vector{Float64}`: edge lengths, unique and sorted.
* `n::Int`: number of intervals in left diagram.
* `m::Int`: number of intervals in right diagram.
"""
struct BottleneckGraph
    adj::Matrix{Float64}

    match_left::Vector{Int}
    match_right::Vector{Int}

    edges::Vector{Float64}

    n_vertices::Int
end

function BottleneckGraph(left::PersistenceDiagram, right::PersistenceDiagram)
    n = length(left)
    m = length(right)
    adj = _adjacency_matrix(left, right)

    edges = filter!(isfinite, sort!(unique!(copy(vec(adj)))))

    return BottleneckGraph(adj, fill(0, n + m), fill(0, m + n), edges, n + m)
end

function _left_neighbors!(buff, graph::BottleneckGraph, vertices, ε, pred)
    empty!(buff)
    for l in vertices
        for r in axes(graph.adj, 1)
            graph.adj[r, l] ≤ ε && pred(r) && push!(buff, r)
        end
    end
    return unique!(buff)
end

function _right_neighbors!(buff, graph::BottleneckGraph, vertices)
    empty!(buff)
    for r in vertices
        push!(buff, graph.match_right[r])
    end
    return buff
end

_is_exposed_right(graph::BottleneckGraph, r) = graph.match_right[r] == 0
_exposed_left(graph::BottleneckGraph) = findall(iszero, graph.match_left)

"""
    _depth_layers(graph::BottleneckGraph, ε)

Split `graph` into layers by how deep they are from a bfs starting at exposed left
vertices in `graph` only taking into account edges of length smaller than or equal to `ε`.
Return depts of right vertices and maximum depth reached.
"""
function _depth_layers(graph::BottleneckGraph, ε)
    depths = fill(0, graph.n_vertices)
    visited = fill(false, graph.n_vertices)
    lefts = _exposed_left(graph)
    rights = Int[]
    i = 1
    while true
        _left_neighbors!(rights, graph, lefts, ε, r -> !visited[r])
        visited[rights] .= true
        depths[rights] .= i
        if isempty(rights)
            # no augmenting path exists
            return nothing, nothing
        elseif any(r -> _is_exposed_right(graph, r), rights)
            return depths, i
        else
            _right_neighbors!(lefts, graph, rights)
        end
        i += 1
    end
end

"""
    _augmenting_paths(graph::BottleneckGraph, ε)

find a maximal set of augmenting paths in graph, taking only edges with weight less than or
equal to `ε` into account.
"""
function _augmenting_paths(graph::BottleneckGraph, ε)
    depths, max_depth = _depth_layers(graph, ε)
    paths = Vector{Int}[]
    isnothing(depths) && return paths

    prev = fill(0, graph.n_vertices)
    rights = Int[]
    lefts = Int[]
    stack = Tuple{Int,Int}[]

    for l_start in _exposed_left(graph)
        empty!(stack)
        push!(stack, (l_start, 1))
        prev .= 0

        while !isempty(stack)
            l, i = pop!(stack)
            parent = graph.match_left[l]
            _left_neighbors!(rights, graph, l, ε, r -> depths[r] == i)
            if i < max_depth
                prev[rights] .= l
                _right_neighbors!(lefts, graph, rights)
                append!(stack, (l, i + 1) for l in lefts)
            else
                found_path = false
                for r in rights
                    if _is_exposed_right(graph, r)
                        prev[r] = l
                        path = Int[r]
                        depths[r] = 0

                        while (l = prev[r]) ≠ l_start
                            @assert prev[r] ≠ 0
                            r = graph.match_left[l]
                            depths[r] = 0
                            append!(path, (l, r))
                        end
                        push!(path, l_start)
                        reverse!(path)
                        push!(paths, path)
                        found_path = true
                        break
                    end
                end
                found_path && break
            end
        end
    end

    return paths
end

function _unmatch_all!(graph::BottleneckGraph)
    graph.match_left .= 0
    graph.match_right .= 0
    return (0, 0)
end

function _augment!(graph, p)
    for i in 1:2:(length(p) - 1)
        l, r = p[i], p[i + 1]
        graph.match_left[l] = r
        graph.match_right[r] = l
    end
end

function _hopcroft_karp!(graph, ε)
    _unmatch_all!(graph)
    paths = _augmenting_paths(graph, ε)
    while !isempty(paths)
        for p in paths
            _augment!(graph, p)
        end
        paths = _augmenting_paths(graph, ε)
    end
    matching =
        [i => graph.match_left[i] for i in 1:(graph.n_vertices) if graph.match_left[i] ≠ 0]
    is_perfect = length(matching) == graph.n_vertices

    return matching, is_perfect
end

"""
    Bottleneck

Use this object to find the bottleneck distance or matching between persistence diagrams.
The distance value is equal to

```math
W_\\infty(X, Y) = \\inf_{\\eta:X\\rightarrow Y} \\sup_{x\\in X} ||x-\\eta(x)||_\\infty,
```

where ``X`` and ``Y`` are the persistence diagrams and ``\\eta`` is a perfect matching
between the intervals. Note the ``X`` and ``Y`` don't need to have the same number of
points, as the diagonal points are considered in the matching as well.

# Warning

Computing the bottleneck distance requires ``\\mathcal{O}(n^2)`` space. Be careful when
computing distances between very large diagrams!

# Usage

* `Bottleneck()(left, right[; matching=false])`: find the bottleneck matching (if
  `matching=true`) or distance (if `matching=false`) between persistence diagrams `left` and
  `right`

# Example

```jldoctest
julia> left = PersistenceDiagram([(1.0, 2.0), (5.0, 8.0)]);

julia> right = PersistenceDiagram([(1.0, 2.0), (3.0, 4.0), (5.0, 10.0)]);

julia> Bottleneck()(left, right)
2.0

julia> Bottleneck()(left, right; matching=true)
Bottleneck Matching with weight 2.0:
 [5.0, 8.0) => [5.0, 10.0)

```
"""
struct Bottleneck <: MatchingDistance end

function (::Bottleneck)(left::PersistenceDiagram, right::PersistenceDiagram; matching=false)
    if count(!isfinite, left) ≠ count(!isfinite, right)
        if matching
            return Matching(left, right, Inf, Pair{Int,Int}[], true)
        else
            return Inf
        end
    end

    graph = BottleneckGraph(left, right)
    edges = graph.edges

    lo = 1
    hi = length(edges)
    while lo < hi - 1
        m = lo + ((hi - lo) >>> 0x01)
        _, succ = _hopcroft_karp!(graph, edges[m])
        if succ
            hi = m
        else
            lo = m
        end
    end
    match, succ = _hopcroft_karp!(graph, edges[lo])
    distance = edges[lo]
    if !succ
        distance = edges[hi]
        match, _ = _hopcroft_karp!(graph, edges[hi])
    end
    @assert length(match) == length(left) + length(right)
    if matching
        return Matching(left, right, distance, match, true)
    else
        return distance
    end
end

function (b::Bottleneck)(left, right; matching=false)
    if length(left) ≠ length(right)
        throw(ArgumentError("`left` and `right` must have the same length"))
    end
    results = (b(l, r; matching=matching) for (l, r) in zip(left, right))
    if matching
        return collect(results)
    else
        return maximum(results)
    end
end

"""
    Wasserstein(q=1)

Use this object to find the Wasserstein distance or matching between persistence diagrams.
The distance value is equal to

```math
W_q(X,Y)=\\left[\\inf_{\\eta:X\\rightarrow Y}\\sum_{x\\in X}||x-\\eta(x)||_\\infty^q\\right],
```

where ``X`` and ``Y`` are the persistence diagrams and ``\\eta`` is a perfect matching
between the intervals. Note the ``X`` and ``Y`` don't need to have the same number of
points, as the diagonal points are considered in the matching as well.

# Warning

Computing the Wasserstein distance requires ``\\mathcal{O}(n^2)`` space. Be careful when
computing distances between very large diagrams!

# Usage

* `Wasserstein(q=1)(left, right[; matching=false])`: find the Wasserstein matching (if
  `matching=true`) or distance (if `matching=false`) between persistence diagrams `left` and
  `right`.

# Example

```jldoctest
julia> left = PersistenceDiagram([(1.0, 2.0), (5.0, 8.0)]);

julia> right = PersistenceDiagram([(1.0, 2.0), (3.0, 4.0), (5.0, 10.0)]);

julia> Wasserstein()(left, right)
3.0

julia> Wasserstein()(left, right; matching=true)
Matching with weight 3.0:
 [1.0, 2.0) => [1.0, 2.0)
 [3.0, 3.0) => [3.0, 4.0)
 [5.0, 8.0) => [5.0, 10.0)

```
"""
struct Wasserstein <: MatchingDistance
    q::Float64

    Wasserstein(q=1) = new(Float64(q))
end

function (w::Wasserstein)(
    left::PersistenceDiagram, right::PersistenceDiagram; matching=false
)
    if count(!isfinite, left) == count(!isfinite, right)
        adj = _adjacency_matrix(right, left, w.q)
        match = collect(i => j for (i, j) in enumerate(hungarian(adj)[1]))
        distance = sum(adj[i, j] for (i, j) in match)^(1 / w.q)

        if matching
            return Matching(left, right, distance, match, false)
        else
            return distance
        end
    else
        if matching
            return Matching(left, right, Inf, Pair{Int,Int}[], false)
        else
            return Inf
        end
    end
end

function (w::Wasserstein)(left, right; matching=false)
    if length(left) ≠ length(right)
        throw(ArgumentError("`left` and `right` must have the same length"))
    end
    results = (w(l, r; matching=matching) for (l, r) in zip(left, right))
    if matching
        return collect(results)
    else
        return sum(results)
    end
end


text:
using ModelingToolkit, ModelingToolkitStandardLibrary, OrdinaryDiffEq
using ModelingToolkitStandardLibrary.Blocks

@parameters t

#=
Testing strategy:
The general strategy is to test systems using simple intputs where the solution
is known on closed form. For algebraic systems (without differential variables),
an integrator with a constant input is often used together with the system under test.
=#

@testset "Constant" begin
    @named c = Constant(; k=1)
    @named int = Integrator()
    @named iosys = ODESystem(connect(c.output, int.input), t, systems=[int, c])
    sys = structural_simplify(iosys)
    prob = ODEProblem(sys, Pair[int.x=>1.0], (0.0, 1.0))
    sol = solve(prob, Rodas4())
    @test sol[int.output.u][end] ≈ 2
end

@testset "Derivative" begin
    @named source = Sine(; frequency=1)
    @named int = Integrator(; k=1)
    @named der = Derivative(; k=1, T=0.001)
    @named iosys = ODESystem([
        connect(source.output, der.input),
        connect(der.output, int.input),
        ],
        t,
        systems=[int, source, der],
    )
    sys = structural_simplify(iosys)
    prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 10.0))
    sol = solve(prob, Rodas4())
    @test isapprox(sol[source.output.u], sol[int.output.u], atol=1e-1)
end

@testset "PT1" begin
    pt1_func(t, k, T) = k * (1 - exp(-t / T)) # Known solution to first-order system

    k, T = 1.0, 0.1
    @named c = Constant(; k=1)
    @named pt1 = FirstOrder(; k=k, T=T)
    @named iosys = ODESystem(connect(c.output, pt1.input), t, systems=[pt1, c])
    sys = structural_simplify(iosys)
    prob = ODEProblem(sys, Pair[], (0.0, 100.0))
    sol = solve(prob, Rodas4())
    @test sol[pt1.output.u] ≈ pt1_func.(sol.t, k, T) atol=1e-3
end

@testset "PT2" begin
    # Known solution to second-order system
    function pt2_func(t, k, w, d)
        y = if d==0
            -k*(-1 + cos(t*w))
        else
            d = complex(d)
            real(k*(1 + (-cosh(sqrt(-1 + d^2)*t*w) - (d*sinh(sqrt(-1 + d^2)*t*w))/sqrt(-1 + d^2))/exp(d*t*w)))
        end
    end

    k, w, d = 1.0, 1.0, 0.5
    @named c = Constant(; k=1)
    @named pt2 = SecondOrder(; k=k, w=w, d=d)
    @named iosys = ODESystem(connect(c.output, pt2.input), t, systems=[pt2, c])
    sys = structural_simplify(iosys)
    prob = ODEProblem(sys, Pair[], (0.0, 100.0))
    sol = solve(prob, Rodas4())
    @test sol[pt2.output.u] ≈ pt2_func.(sol.t, k, w, d) atol=1e-3
end

@testset "StateSpace" begin
    A = [0 1;-1 -0.5]
    B = [0, 1]
    C = [0.9 1;]
    D = [0;;]
    @named ss = StateSpace(;A,B,C,D,x_start=zeros(2))
    @named c = Constant(; k=1)
    @named model = ODESystem(
        [
           connect(c.output, ss.input),
        ],
        t,
        systems=[ss, c]
    )
    sys = structural_simplify(model)
    prob = ODEProblem(sys, Pair[], (0.0, 100.0))
    sol = solve(prob, Rodas4())
    # equilibrium point is at [1, 0]
    @test sol[ss.x[1]][end] ≈ 1 atol=1e-3
    @test sol[ss.x[2]][end] ≈ 0 atol=1e-3
end

"""Second order demo plant"""
function Plant(;name, x_start=zeros(2))
    @named input = RealInput()
    @named output = RealOutput()
    D = Differential(t)
    sts = @variables x1(t)=x_start[1] x2(t)=x_start[2]
    eqs= [
        D(x1) ~ x2
        D(x2) ~ -x1 - 0.5 * x2 + input.u
        output.u ~ 0.9 * x1 + x2
    ]
    compose(ODESystem(eqs, t, sts, []; name), [input, output])
end

@testset "PI" begin
    @named ref = Constant(; k=2)
    @named pi_controller = PI(k=1, T=1)
    @named plant = Plant()
    @named fb = Feedback()
    @named model = ODESystem(
        [
            connect(ref.output, fb.input1),
            connect(plant.output, fb.input2),
            connect(fb.output, pi_controller.err_input),
            connect(pi_controller.ctr_output, plant.input),
        ],
        t,
        systems=[pi_controller, plant, ref, fb]
    )
    sys = structural_simplify(model)
    prob = ODEProblem(sys, Pair[], (0.0, 100.0))
    sol = solve(prob, Rodas4())
    @test sol[ref.output.u - plant.output.u][end] ≈ 0 atol=1e-3 # zero control error after 100s
end

@testset "PID" begin
    @named ref = Constant(; k=2)
    @named pid_controller = PID(k=3, Ti=0.5, Td=100)
    @named plant = Plant()
    @named fb = Feedback()
    @named model = ODESystem(
        [
            connect(ref.output, fb.input1),
            connect(plant.output, fb.input2),
            connect(fb.output, pid_controller.err_input),
            connect(pid_controller.ctr_output, plant.input),
        ],
        t,
        systems=[pid_controller, plant, ref, fb]
    )
    sys = structural_simplify(model)
    prob = ODEProblem(sys, Pair[], (0.0, 100.0))
    sol = solve(prob, Rodas4())
    @test sol[ref.output.u - plant.output.u][end] ≈ 0 atol=1e-3 # zero control error after 100s
end

@test_skip begin
    k = 2
    Ti = 0.5
    Td = 0.7
    wp = 1
    wd = 1
    Ni = √(Td / Ti)
    Nd = 12
    y_max = Inf
    y_min = -Inf
    u_r = sin(t)
    u_y = 0
    function solve_with_input(; u_r, u_y,
        controller = PID(; k, Ti, Td, wp, wd, Ni, Nd, y_max, y_min, name=:controller)
    )
        @test count(ModelingToolkit.isinput, states(controller)) == 5 # 2 in PID, 1 sat, 1 I, 1 D
        @test count(ModelingToolkit.isoutput, states(controller)) == 4
        # TODO: check number of unbound inputs when available, should be 2
        @named iosys = ODESystem([controller.u_r~u_r, controller.u_y~u_y], t, systems=[controller])
        sys = structural_simplify(iosys)
        prob = ODEProblem(sys, Pair[], (0.0, 10.0))
        sol = solve(prob, Rodas4(), saveat=0:0.1:10)
        controller, sys, sol
    end

    # linearity in u_r
    controller, sys, sol1 = solve_with_input(u_r=sin(t), u_y=0)
    controller, sys, sol2 = solve_with_input(u_r=2sin(t), u_y=0)
    @test sum(abs, sol1[controller.ea]) < eps() # This is the acutator model error due to saturation
    @test 2sol1[controller.y] ≈ sol2[controller.y] rtol=1e-3 # linearity in u_r

    # linearity in u_y
    controller, sys, sol1 = solve_with_input(u_y=sin(t), u_r=0)
    controller, sys, sol2 = solve_with_input(u_y=2sin(t), u_r=0)
    @test sum(abs, sol1[controller.ea]) < eps() # This is the acutator model error due to saturation
    @test 2sol1[controller.y] ≈ sol2[controller.y] rtol=1e-3 # linearity in u_y

    # zero error
    controller, sys, sol1 = solve_with_input(u_y=sin(t), u_r=sin(t))
    @test sum(abs, sol1[controller.y]) ≈ 0 atol=sqrt(eps())

    # test saturation
    controller, sys, sol1 = solve_with_input(; u_r=10sin(t), u_y=0,
        controller = PID(; k, Ti, Td, wp, wd=0, Ni, Nd, y_max=10, y_min=-10, name=:controller)
    )
    @test extrema(sol1[controller.y]) == (-10, 10)

    # test P set-point weighting
    controller, sys, sol1 = solve_with_input(; u_r=sin(t), u_y=0,
        controller = PID(; k, Ti, Td, wp=0, wd, Ni, Nd, y_max, y_min, name=:controller)
    )
    @test sum(abs, sol1[controller.ep]) ≈ 0 atol=sqrt(eps())

    # test D set-point weighting
    controller, sys, sol1 = solve_with_input(; u_r=sin(t), u_y=0,
        controller = PID(; k, Ti, Td, wp, wd=0, Ni, Nd, y_max, y_min, name=:controller)
    )
    @test sum(abs, sol1[controller.ed]) ≈ 0 atol=sqrt(eps())

    # zero integral gain
    controller, sys, sol1 = solve_with_input(; u_r=sin(t), u_y=0,
        controller = PID(; k, Ti=false, Td, wp, wd, Ni, Nd, y_max, y_min, name=:controller)
    )
    @test isapprox(sum(abs, sol1[controller.I.y]), 0, atol=sqrt(eps()))

    # zero derivative gain
    @test_skip begin # During the resolution of the non-linear system, the evaluation of the following equation(s) resulted in a non-finite number: [5]
        controller, sys, sol1 = solve_with_input(; u_r=sin(t), u_y=0,
            controller = PID(; k, Ti, Td=false, wp, wd, Ni, Nd, y_max, y_min, name=:controller)
        )
        @test isapprox(sum(abs, sol1[controller.D.y]), 0, atol=sqrt(eps()))
    end

    # Tests below can be activated when the concept of unbound_inputs exists in MTK
    # @test isequal(Set(unbound_inputs(controller)), @nonamespace(Set([controller.u_r, controller.u_y])))
    # @test isempty(unbound_inputs(sys))
    # @test isequal(bound_inputs(sys), inputs(sys))
    # @test isequal(
    #     Set(bound_inputs(sys)),
    #     Set([controller.u_r, controller.u_y, controller.I.u, controller.D.u, controller.sat.u])
    #     )
end

@testset "LimPI" begin
    @named ref = Constant(; k=1)
    @named pi_controller_lim = LimPI(k=3, T=0.5, u_max=1.5, u_min=-1.5, Ta=0.1)
    @named pi_controller = PI(k=3, T=0.5)
    @named sat = Limiter(y_max=1.5, y_min=-1.5)
    @named plant = Plant()
    @named fb = Feedback()

    # without anti-windup measure
    sol = let
        @named model = ODESystem(
            [
                connect(ref.output, fb.input1),
                connect(plant.output, fb.input2),
                connect(fb.output, pi_controller.err_input),
                connect(pi_controller.ctr_output, sat.input),
                connect(sat.output, plant.input),
            ],
            t,
            systems=[pi_controller, plant, ref, fb, sat]
        )
        sys = structural_simplify(model)
        prob = ODEProblem(sys, Pair[], (0.0, 20.0))
        sol = solve(prob, Rodas4())
    end

    # with anti-windup measure
    sol_lim = let
        @named model = ODESystem(
            [
                connect(ref.output, fb.input1),
                connect(plant.output, fb.input2),
                connect(fb.output, pi_controller_lim.err_input),
                connect(pi_controller_lim.ctr_output, sat.input),
                connect(sat.output, plant.input),
            ],
            t,
            systems=[pi_controller_lim, plant, ref, fb, sat]
        )
        sys = structural_simplify(model)
        prob = ODEProblem(sys, Pair[], (0.0, 20.0))
        sol = solve(prob, Rodas4())
    end

    @test sol[ref.output.u - plant.output.u][end] ≈ 0 atol=1e-3 # zero control error after 100s
    @test sol_lim[ref.output.u - plant.output.u][end] ≈ 0 atol=1e-3 # zero control error after 100s

    # Plots.plot(sol; vars=[plant.output.u]) # without anti-windup measure
    # Plots.plot!(sol_lim; vars=[plant.output.u]) # with anti-windup measure
end

@testset "LimPID" begin
    @named ref = Constant(; k=1)
    @named pid_controller = LimPID(k=3, Ti=0.5, Td=100, u_max=1.5, u_min=-1.5, Ni=0.1/0.5)
    @named plant = Plant()
    @named model = ODESystem(
        [
            connect(ref.output, pid_controller.reference),
            connect(plant.output, pid_controller.measurement),
            connect(pid_controller.ctr_output, plant.input),
        ],
        t,
        systems=[pid_controller, plant, ref]
    )
    sys = structural_simplify(model)
    prob = ODEProblem(sys, Pair[], (0.0, 100.0))
    sol = solve(prob, Rodas4())

    # Plots.plot(sol, vars=[plant.output.u, plant.input.u])
    @test sol[ref.output.u - plant.output.u][end] ≈ 0 atol=1e-3 # zero control error after 100s
end

ext: jl


text:
# <editor-fold desc="definition and handling of parameters"

# XXX defines struct for handling parameter data
"""
Type including data and additional information on parameters. Fields relate to what is provided in [Parameter list](@ref) and include:
* `name::Symbol`: name of the parameter
* `dim::Tuple`: potential dimensions of parameter data
* `defVal::Union{Nothing,Float64}`: default value
* `herit::Tuple`: inheritance rules for parameter, see [Parameter overview](@ref) for details
* `data::DataFrame`: specified parameter data
"""
mutable struct ParElement
	name::Symbol
    dim::Tuple
    defVal::Union{Nothing,Float64}
    herit::Tuple
    data::DataFrame
	techPre::NamedTuple{(:preset,:mode),Tuple{Symbol,Tuple{Vararg{Symbol,N} where N}}}

	function ParElement(paraData_df::DataFrame,paraDef_ntup::NamedTuple,name::Symbol,report::Array{Tuple,1})

        setLongShort_dic = Dict(:Ts => :timestep, :R => :region, :C => :carrier, :Te => :technology, :M => :mode)
		if isempty(paraData_df) return new(name,paraDef_ntup.dim,paraDef_ntup.defVal,paraDef_ntup.herit,DataFrame()) end

        # XXX check consistency of rows in input dataframe and definition of set and rename columns according to set defintion
        # assigns array of used suffixes according to parameter defintion to each set
        splitDim_arr = map(x -> map(y -> Symbol(y), split(String(x),"_")),collect(paraDef_ntup.dim))
        setSuf_dic = Dict(x => map(y -> length(y) == 1 ? Symbol() : y[end],filter(z -> z[1] == x,splitDim_arr)) for x in unique(map(x -> x[1],splitDim_arr)))

        # loops over set columns in input dataframe and assigns them to the sets defined for the parameter
        newCol_dic = Dict(:val => :val)
        sufNum_dic = Dict(:b => 2, :c => 3, :d => 4, :e => 5, :f => 6, :g => 7)

        for colNam in setdiff(namesSym(paraData_df),[:val])
            colNam_arr = split(String(colNam),"_")
            setNam = Symbol(colNam_arr[1])

            if !haskey(setSuf_dic,setNam) # parameter provided for a set not appearing in definition (e.g. demand depending on the technology)
                push!(report,(2, "parameter assignment", string(name), "parameter data was specified for $(setLongShort_dic[setNam]) set, but it is not defined to depend on this set"))
                continue
            elseif length(setSuf_dic[setNam]) == 1 && length(colNam_arr) > 1 # they are several instances of the set provided, but it only depends on one instance (e.g. two region sets for efficiency)
                push!(report,(2, "parameter assignment", string(name), "parameter data was specified for several instances of $(setLongShort_dic[setNam]) set, but it is defined to depend only on one instance, additonal instances were ignored"))
                continue
            elseif setSuf_dic[setNam][1] == Symbol() # set has only one instance and no suffix => no change when converting from read-in dataframe to parameter element
                newCol_dic[colNam] = colNam
            elseif length(setSuf_dic[setNam]) == 1 || length(colNam_arr) == 1 # column name in dataframe has no underscore, but defintion of parameter element has one
                newCol_dic[colNam] = Symbol(setNam, "_", setSuf_dic[setNam][1])
            else
                cntRep_int = sufNum_dic[Symbol(colNam_arr[2])] # set defined for several instances
                newCol_dic[colNam] = Symbol(setNam, "_", setSuf_dic[setNam][cntRep_int])
            end
        end

        # filters only used coulumns, renames them accordingly and converts to table
        writeData_df = paraData_df[:,collect(keys(newCol_dic))]
        DataFrames.rename!(writeData_df,newCol_dic)

        new_obj = new(name,paraDef_ntup.dim,paraDef_ntup.defVal,paraDef_ntup.herit,writeData_df)
        # defines on which level parameter is presetted and which capacity restrictions are affected by different modes for all dispatch parameters, where this is specified
        if haskey(paraDef_ntup,:techPre) new_obj.techPre = paraDef_ntup.techPre end
        return new_obj
    end

	ParElement() = new()
end

# XXX specific struct for read in process of parameter data
mutable struct parEntry
	colSet::Symbol
	entry::Array{String,1}
	lvl::Array{Int,1}
	startLvl::Int
end

# </editor-fold>

# <editor-fold desc="import and extensions of base functions"

# XXX functions to copy parameter structs of parameter data
import Base.copy
function copy(par_obj::ParElement)
	out = ParElement()
	out.name = par_obj.name
	out.dim = par_obj.dim
	out.defVal = par_obj.defVal
	out.herit = par_obj.herit
	out.data = copy(par_obj.data)
	if isdefined(par_obj,:techPre) out.techPre = par_obj.techPre end
	return out
end

function copy(par_obj::ParElement,data_df::DataFrame)
	out = ParElement()
	out.name = par_obj.name
	out.dim = par_obj.dim
	out.defVal = par_obj.defVal
	out.herit = par_obj.herit
	out.data = data_df
	if isdefined(par_obj,:techPre) out.techPre = par_obj.techPre end
	return out
end

# XXX usual collect sometimes creates a mysterious error if used on dictionary keys, this command avoids this
import Base._collect
import Base.SizeUnknown

collectKeys(itr) = _collect(Symbol, itr, SizeUnknown())

# </editor-fold>

# </editor-fold>

# <editor-fold desc="struct for individual parts of the model"

# XXX defines parts of the model
abstract type AbstractModelPart end

"""
```julia
	TechPart <: AbstractModelPart
```
Type used for technology parts. Parameters, variables, and constraints are assigned as dictionaries via the fields `par`, `var`, and  `cns`, respectively. Additional fields include:
* `name::Tuple`: name of technology as a series of nodes from the technology tree
* `carrier::NamedTuple`: energy carriers by index assigned to technology by groups (e.g. generation, use, ...)
* `balLvl::NamedTuple`: temporal and spatial resolution for expansion and balance of the technology
* `capaRestr::DataFrame`: specification of capacity restrictions required for technology
* `actSt::Tuple`: actively stored carriers  altough they are not leafs by index
* `type::Tuple`: type of technology (stock, mature, or evolving)
* `disAgg::Bool`: if true, dispatch is modelled at expansion resolution instead of dispatch resolution
* `modes::Tuple`: different operational modes of technology
"""
mutable struct TechPart <: AbstractModelPart
	name::Tuple{Vararg{String,N} where N}
	par::Dict{Symbol,ParElement}
	var::Dict{Symbol,DataFrame}
	cns::Dict{Symbol,DataFrame}
	carrier::NamedTuple
	balLvl::NamedTuple{(:exp,:ref),Tuple{Tuple{Int,Int},Union{Nothing,Tuple{Int,Int}}}}
	capaRestr::DataFrame
	actSt::Tuple
	type::Symbol
	disAgg::Bool
	modes::Tuple{Vararg{Int,N} where N}
	TechPart(name::Tuple{Vararg{String,N} where N}) = new(name,Dict{Symbol,ParElement}(),Dict{Symbol,DataFrame}(),Dict{Symbol,DataFrame}())
	TechPart() = new()
end

"""
```julia
	OthPart <: AbstractModelPart
```
Type used for 'exchange', 'trade', 'balance', 'limits', and 'objectives' model parts. Parameters, variables, and constraints are assigned as dictionaries via the fields `par`, `var`, and  `cns`, respectively.
"""
mutable struct OthPart <: AbstractModelPart
	par::Dict{Symbol,ParElement}
	var::Dict{Symbol,DataFrame}
	cns::Dict{Symbol,DataFrame}
	OthPart() = new(Dict{Symbol,ParElement}(),Dict{Symbol,DataFrame}(),Dict{Symbol,DataFrame}())
end

# XXX container to store data defining a constraint (used to seperate definition and actual jump creation of constraints)
struct cnsCont
    data::DataFrame
    sign::Symbol
end

# </editor-fold>

# <editor-fold desc="structs for nodes that then make up the trees to save set data"

# XXX define nodes for set tree and tree itself
mutable struct Node
	idx::Int
	val::String
	lvl::Int
	subIdx::Int
	down::Array{Int,1}
end

mutable struct Tree
	nodes::Dict{Int,Node}
	srcTup::Dict{Tuple,Array{Int,1}}
	srcStr::Dict{Tuple{String,Int},Array{Int,1}}
	up::Dict{Int,Int}
	height::Int
	Tree() = new(Dict{Int,Node}(),Dict{Tuple,Int}(),Dict{String,Array{Int,1}}(),Dict{Int,Int}(),1)
end

# </editor-fold>

# <editor-fold desc="options for model and model itself"

# create abstract model object to reference before creation (avoid circular type definiton)
abstract type AbstractModel end

# XXX defines final model object and its options
struct modOptions
	# data in- and output
	inDir::Array{String,1}
	outDir::String
	objName::String
	csvDelim::String
	outStamp::String
	# model generation
	decomm::Symbol
	interCapa::Symbol
	supTsLvl::Int
	shortExp::Int
	redStep::Float64
	# managing numerical issues
	emissionLoss::Bool
	coefRng::NamedTuple{(:mat,:rhs),Tuple{Tuple{Float64,Float64},Tuple{Vararg{Float64,2}}}}
	scaFac::NamedTuple{(:capa,:oprCapa,:dispConv,:dispSt,:dispExc, :dispTrd, :costDisp,:costCapa,:obj),Tuple{Vararg{Float64,9}}}
	bound::NamedTuple{(:capa,:disp,:obj),Tuple{Vararg{Float64,3}}}
	avaMin::Float64
	checkRng::Float64
	# reporting related options
	reportLvl::Int
	errCheckLvl::Int
	errWrtLvl::Int
	startTime::DateTime
end

# XXX flow graph object that defines relations between technologies and carriers (and among carriers)
mutable struct flowGraph
	nodeC::Dict{Int64,Int64}
	nodeTe::Dict{Int64,Int64}
	edgeC::Array{Pair{Int,Int},1}
	edgeTe::Array{Pair{Int,Int},1}
	nodePos::Dict{Int,Array{Float64,1}}

	function flowGraph(anyM::AbstractModel)

		# creates dictionary mapping carrier id to node id
		nodeC_dic = Dict(x[2] => x[1] for x in enumerate(sort(filter(x -> x != 0,getfield.(collect(values(anyM.sets[:C].nodes)),:idx)))))

		# get all relevant technology, a technology is not relevant, where all children of a parent have the same carriers (in this case only the parent is relevant)
		t_tree = anyM.sets[:Te]

		allTech_arr =  getfield.(collect(values(t_tree.nodes)),:idx)
		tleaf_dic = Dict(x => unique(filter(y -> y in keys(anyM.parts.tech), [x,getDescendants(x,t_tree,true)...])) for x in allTech_arr)
		relTech_arr = Array{Int,1}()

		for t in keys(tleaf_dic)
		    subCar_arr = map(y -> anyM.parts.tech[y].carrier,tleaf_dic[t])
		    if length(unique(subCar_arr)) == 1
		        push!(relTech_arr,t)
		    else
		        append!(relTech_arr,collect(tleaf_dic[t]))
		    end
		end

		# creates dictionary mapping each relevant id to node id
		nodeTe_dic = Dict(x[2] => x[1] + length(nodeC_dic) for x in enumerate(filter(x -> isempty(intersect(getAncestors(x,t_tree,:int),relTech_arr)),unique(relTech_arr))))

		# creates edges between technologies
		edgeTe_arr = Array{Pair{Int,Int},1}()

		for t in keys(nodeTe_dic)
		    gotTech_boo = false; tItr = t
		    while !gotTech_boo
		        if tItr in keys(anyM.parts.tech)
		            gotTech_boo = true
		        else
		            tItr = intersect(getDescendants(t,anyM.sets[:Te],true),keys(anyM.parts.tech))[1]
		        end
		    end

		    car_ntup = anyM.parts.tech[tItr].carrier

		    for cIn in map(x -> getfield(car_ntup,x),intersect(keys(car_ntup),(:use,:stExtIn))) |> (y -> isempty(y) ? y : union(y...))
		        push!(edgeTe_arr, nodeC_dic[cIn] => nodeTe_dic[t])
		    end

		    for cOut in map(x -> getfield(car_ntup,x),intersect(keys(car_ntup),(:gen,:stExtOut))) |> (y -> isempty(y) ? y : union(y...))
		        push!(edgeTe_arr, nodeTe_dic[t] => nodeC_dic[cOut])
		    end
		end

		# creates edges between carriers
		edgeC_arr = Array{Pair{Int,Int},1}()

		for c in keys(nodeC_dic)
		    for cChild in anyM.sets[:C].nodes[c].down
		        push!(edgeC_arr, nodeC_dic[cChild] => nodeC_dic[c])
		    end
		end

		return new(nodeC_dic,nodeTe_dic,edgeC_arr,edgeTe_arr)
	end
end

# XXX specific information for graphical evaluation
mutable struct graInfo
	graph::flowGraph
	names::Dict{String,String}
	colors::Dict{String,Tuple{Float64,Float64,Float64}}

	function graInfo(anyM::AbstractModel)
		# create default options for names and colors
		graph_obj = flowGraph(anyM)

		# specificy some default names and colors used in visualisations
		namesDef_arr = ["coalPlant" => "coal plant", "gasPlant" => "gas plant", "districtHeat" => "district heat", "naturalGas" => "natural gas", "synthGas" => "synthetic gas", "fossilGas" => "fossil gas",
									"demand" => "final demand", "export" => "export", "import" => "import", "crt" => "curtailment", "lss" => "loss of load", "trdSell" => "trade sell", "trdBuy" => "trade buy"]

		# create dictionary assigning internal model names to names used within visualisations
		allVal_arr = unique(vcat(map(x -> getfield.(values(anyM.sets[x].nodes),:val) ,collect(keys(anyM.sets)))...))
		names_dic = setdiff(allVal_arr,getindex.(namesDef_arr,1))  |> (z -> Dict(vcat(namesDef_arr,Pair.(z,z))))

		# define default colors for default energy carriers
		colorsCar_arr = ["electricity" => (1.0, 0.9215, 0.2313),"heat" => (0.769,0.176,0.290),"districtHeat" => (0.6,0.0,0.169), "gas" => (1.0,0.416,0.212),
									"naturalGas" => (1.0,0.506,0.294),"fossilGas" => (0.898,0.259,0.075), "synthGas" => (0.235,0.506,0.325), "hydrogen" => (0.329,0.447,0.827),
													"coal" => (0.459,0.286,0.216),"biomass" => (0.682,0.898,0.443),"bioGas" => (0.682,0.898,0.443)]

		colors_dic = setdiff(getfield.(values(anyM.sets[:C].nodes),:val),getindex.(colorsCar_arr,1)) |> (z -> Dict(vcat(colorsCar_arr,Pair.(z,fill((0.85,0.85,0.85),length(z))))))

		return new(graph_obj,names_dic,colors_dic)
	end
end

# XXX finally, the model object itself
"""
The core model object containing all related data and subordinate objects.
#  Constructor and arguments
```julia
anyModel(inDir::Union{String,Array{String,1}},outDir::String; kwargs)
```
* `inDir::Union{String,Array{String,1}}`: directory of input files, also allows for provide multiple directories via an array
* `outDir::String`: directory of output files are written to
# Optional arguments, data handling
* `objName::String`: name of the model object, will be added to the name of output files and printed during reporting, default is an empty string
* `csvDelim::String`: specifies the delimiter used within the read-in csv files, default is a comma `,`
# Optional arguments, model generation
* `decomm::Symbol`: specifies if the model should perform endogenous decommissioning, options are:
    - `:decomm`: capacities are decommissioned endogenously, once decommissioned capacities cannot be put into operation again (default)
    - `:none`: no endogenous decommissioning, operated capacities equal installed capacities
    - `:recomm`: capacities are decommissioned endogenously and can be put back into operation
* `interCapa::Symbol`: capacity expansion can be modelled at a resolution less detailed than yearly, this options determines how capacities are distributed among the subsequent years in this case, options are:
    - `:linear`: expansion is equally distributed among years resulting in a linear increase in capacity (default)
    - `:none`: all expansion occurs in the first  year
* `supTsLvl::Int`: specifies the depth in the tree of time-steps that provides years, default is `0`
* `shortExp::Int`: intervall in years between years of capacity expansion, default is `10`
* `redStep::Float64`: scales down energy quantities within the model, can be relevant when working with reduced time-series, default is `1.0`

# Optional arguments, reporting (see [Data files](@ref) for details)
* `reportLvl::Int`: controls the frequency of writing updates to the console, default is `2`
* `errCheckLvl::Int`: controls the frequency of checking for errors, default is `2`
* `errWrtLvl::Int`: controls the frequency of writing an error report to a csv file, default is`1`
# Optional arguments, numerical issues (see [Performance and stability](@ref) for details)
* `coefRng::NamedTuple`: specifies the maximum range of coefficients in the matrix and right-hand side of the model's underlying optimization problem, default is `(mat = (1e-2,1e5), rhs = (1e-2,1e2))`
* `checkRng::Float64`: if set, reports all equations whose range exceeds the specified value, default is `NaN`
* `scaFac::NamedTuple`: scales different groups of variables within the model, default is `(capa = 1e1, oprCapa = 1e2, dispConv = 1e3, dispSt = 1e4, dispExc = 1e3, dispTrd = 1e3, costDisp = 1e1, costCapa = 1e2, obj = 1e0)`
* `bound::NamedTuple`: sets external bounds for all capacities and dispatch variables (both in GW) and for the objective value (in Mil. €), default is `(capa = NaN, disp = NaN, obj = NaN)`
* `avaMin::Float64`: availabilities smaller than this value are set to zero, since in the [Conversion capacity restriction](@ref) and [Storage capacity restriction](@ref) availabilities are inversed, this avoids high coefficients, default is `0.01`
* `emissionLoss::Bool`: determines if losses from exchange and self-discharge of storage are subject to emissions, default is `true`
"""
mutable struct anyModel <: AbstractModel

	options::modOptions
	report::Array{Tuple,1}
	optModel::Model
	lock::ReentrantLock

	supTs::NamedTuple{(:lvl,:step,:sca),Tuple{Int,Tuple{Vararg{Int,N} where N},Dict{Tuple{Int,Int},Float64}}}
	cInfo::Dict{Int,NamedTuple{(:tsDis,:tsExp,:rDis,:rExp,:eq),Tuple{Int,Int,Int,Int,Bool}}}

	sets::Dict{Symbol,Tree}
	parts::NamedTuple{(:tech,:trd,:exc,:bal,:lim,:obj),Tuple{Dict{Int,TechPart},OthPart,OthPart,OthPart,OthPart,OthPart}}

	graInfo::graInfo
	function anyModel(inDir::Union{String,Array{String,1}},outDir::String; objName = "", csvDelim = ",", decomm = :recomm, interCapa = :linear, supTsLvl = 0, shortExp = 10, redStep = 1.0, emissionLoss = true,
																										reportLvl = 2, errCheckLvl = 1, errWrtLvl = 1, coefRng = (mat = (1e-2,1e5), rhs = (1e-2,1e2)),
																											scaFac = (capa = 1e1, oprCapa = 1e2, dispConv = 1e3, dispSt = 1e4, dispExc = 1e3, dispTrd = 1e3, costDisp = 1e1, costCapa = 1e2, obj = 1e0),
																																bound = (capa = NaN, disp = NaN, obj = NaN), avaMin = 0.01, checkRng = NaN)
		anyM = new()

		# <editor-fold desc="initialize report and options"

		# XXX creates dataframe to which reporting is written
		anyM.report = Array{Tuple,1}()

		anyM.optModel = Model()
		anyM.lock = ReentrantLock()

		# XXX sets whole options object from specified directories TODO arbeite mit kwargs später
		outStamp_str = string(objName,"_",Dates.format(now(),"yyyymmddHHMM"))
		defOpt_ntup = (inDir = typeof(inDir) == String ? [inDir] : inDir, outDir = outDir, objName = objName, csvDelim = csvDelim, outStamp = outStamp_str, decomm = decomm, interCapa = interCapa,
																					supTsLvl = supTsLvl, shortExp = shortExp, redStep = redStep, emissionLoss = emissionLoss, coefRng = coefRng, scaFac = scaFac, bound = bound,
																						avaMin = avaMin, checkRng = checkRng, reportLvl = reportLvl, errCheckLvl = errCheckLvl, errWrtLvl = errWrtLvl, startTime = now())

		anyM.options = modOptions(defOpt_ntup...)

		# </editor-fold>

		# <editor-fold desc= read in set and parameter data>
		files_dic = readInputFolder(anyM.options.inDir)

		# XXX read-in sets and parameters
		setData_dic = readSets!(files_dic,anyM)
		paraTemp_dic = readParameters!(files_dic,setData_dic,anyM)

		produceMessage(anyM.options,anyM.report, 1," - Read-in all set and parameter files")

		# </editor-fold>

		# <editor-fold desc="create part objects and general mappings"
		# assign actual tech to parents
		relTech_df = setData_dic[:Te][!,Symbol.(filter(x -> occursin("technology",x) && !isnothing(tryparse(Int16,string(x[end]))), string.(namesSym(setData_dic[:Te]))))]
		relTech_df = DataFrame(filter(x -> any(collect(x) .!= ""), eachrow(relTech_df)))
		techIdx_arr = filter(z -> isempty(anyM.sets[:Te].nodes[z].down), map(x -> lookupTupleTree(tuple(collect(x)...),anyM.sets[:Te],1)[1], eachrow(relTech_df)))

		anyM.parts = (tech = Dict(x => TechPart(getUniName(x,anyM.sets[:Te])) for x in techIdx_arr), trd = OthPart(), exc = OthPart(), bal = OthPart(), lim = OthPart(), obj = OthPart())

		createCarrierMapping!(setData_dic,anyM)
		createTimestepMapping!(anyM)

		# XXX write general info about technologies
		for t in techIdx_arr createTechInfo!(t, setData_dic, anyM) end
		produceMessage(anyM.options,anyM.report, 2," - Created all mappings among sets")

		# XXX assign parameters to model parts
		parDef_dic = parameterToParts!(paraTemp_dic, techIdx_arr, anyM)
		produceMessage(anyM.options,anyM.report, 2," - Assigned parameter data to model parts")

		# XXX create object for data visualization
		anyM.graInfo = graInfo(anyM)

		produceMessage(anyM.options,anyM.report, 1," - Prepared creation of optimzation model")
		# </editor-fold>

		return anyM
	end

	anyModel() = new()
end

# </editor-fold>


text:
include("settings.jl")

settings[:doctest] = false

makedocs(sitename = sitename; settings...)


text:
## ------------------------------------------------------------------
macro new_paragraph!(ex...)
    label = _insert_label!("new_paragraph!", ex, __source__)
    _add_secobj!(() -> RBParagraph(label), __source__)
end


text:
function nlsolve(df::Union{NonDifferentiable, OnceDifferentiable},
                 initial_x::AbstractArray;
                 method::Symbol = :trust_region,
                 xtol::Real = zero(real(eltype(initial_x))),
                 ftol::Real = convert(real(eltype(initial_x)), 1e-8),
                 iterations::Integer = 1_000,
                 store_trace::Bool = false,
                 show_trace::Bool = false,
                 extended_trace::Bool = false,
                 linesearch = LineSearches.Static(),
                 linsolve=(x, A, b) -> copyto!(x, A\b),
                 factor::Real = one(real(eltype(initial_x))),
                 autoscale::Bool = true,
                 m::Integer = 10,
                 beta::Real = 1,
                 aa_start::Integer = 1,
                 droptol::Real = convert(real(eltype(initial_x)), 1e10))
    if show_trace
        @printf "Iter     f(x) inf-norm    Step 2-norm \n"
        @printf "------   --------------   --------------\n"
    end
    if method == :newton
        newton(df, initial_x, xtol, ftol, iterations,
               store_trace, show_trace, extended_trace, linesearch; linsolve=linsolve)
    elseif method == :trust_region
        trust_region(df, initial_x, xtol, ftol, iterations,
                     store_trace, show_trace, extended_trace, factor,
                     autoscale)
    elseif method == :anderson
        anderson(df, initial_x, xtol, ftol, iterations,
                 store_trace, show_trace, extended_trace, m, beta, aa_start, droptol)
    elseif method == :broyden
        broyden(df, initial_x, xtol, ftol, iterations,
                store_trace, show_trace, extended_trace, linesearch)
    else
        throw(ArgumentError("Unknown method $method"))
    end
end

function nlsolve(f,
                 initial_x::AbstractArray;
                 method::Symbol = :trust_region,
                 autodiff = :central,
                 inplace = !applicable(f, initial_x),
                 kwargs...)
    if method in (:anderson, :broyden)
        df = NonDifferentiable(f, initial_x, similar(initial_x); inplace=inplace)
    else
        df = OnceDifferentiable(f, initial_x, similar(initial_x); autodiff=autodiff, inplace=inplace)
    end

    nlsolve(df, initial_x; method = method, kwargs...)
end

function nlsolve(f,
                 j,
                 initial_x::AbstractArray;
                 inplace = !applicable(f, initial_x),
                 kwargs...)
    if inplace
        df = OnceDifferentiable(f, j, initial_x, similar(initial_x))
    else
        df = OnceDifferentiable(not_in_place(f, j)..., initial_x, similar(initial_x))
    end

    nlsolve(df, initial_x; kwargs...)
end

function nlsolve(f,
                 j,
                 fj,
                 initial_x::AbstractArray;
                 inplace = !applicable(f, initial_x),
                 kwargs...)
    if inplace
        df = OnceDifferentiable(f, j, fj, initial_x, similar(initial_x))
    else
        df = OnceDifferentiable(not_in_place(f, j, fj)..., initial_x, similar(initial_x))
    end

    nlsolve(df, initial_x; kwargs...)
end


text:
module Rotations3D

using StaticArrays
using LinearAlgebra: norm, rank, svd, Diagonal, tr

using ACE: coco_zeros, coco_init, coco_dot, coco_filter, AbstractProperty,
			  coco_type

export ClebschGordan, Rot3DCoeffs, ri_basis, rpi_basis, R3DC, Rot3DCoeffsEquiv

"""
`ClebschGordan: ` storing precomputed Clebsch-Gordan coefficients; see
`?clebschgordan` for the convention that is use.
"""
struct ClebschGordan{T}
	vals::Dict{Tuple{Int, Int, Int, Int, Int, Int}, T}
end

# -> CouplingCoeffRecursion ????
struct Rot3DCoeffs{T, TP}
   vals::Vector{Dict}  # val[N] = coeffs for correlation order N
   cg::ClebschGordan{T}
	phi::TP
end

# -----------------------------------
# iterating over an m collection
# -----------------------------------

_mvec(::CartesianIndex{0}) = SVector{0, Int}()

_mvec(mpre::CartesianIndex) = SVector(Tuple(mpre)...)

struct MRange{N, T2, TP}
   ll::SVector{N, Int}
   cartrg::T2
	phi::TP
end

Base.length(mr::MRange) =
		sum(mt -> coco_filter(mr.phi, mr.ll, _mvec(mt)), mr.cartrg)

"""
Given an l-vector `ll` iterate over all combinations of `mm` vectors  of
the same length such that `sum(mm) == 0`
"""
_mrange(phi, ll) =
	MRange(ll, CartesianIndices(ntuple(i -> -ll[i]:ll[i], length(ll))), phi)

# TODO: should we impose here that (ll, mm) are lexicographically ordered?

function Base.iterate(mr::MRange, idx::Integer=0)
	while true
		idx += 1
		if idx > length(mr.cartrg)
			return nothing
		end
		mm = _mvec(mr.cartrg[idx])
		if coco_filter(mr.phi, mr.ll, mm)
			return mm, idx
		end
	end
	error("we should never be here")
end

# ----------------------------------------------------------------------
#     ClebschGordan code
# ----------------------------------------------------------------------

cg_conditions(j1,m1, j2,m2, J,M) =
	cg_l_condition(j1, j2, J)   &&
	cg_m_condition(m1, m2, M)   &&
	(abs(m1) <= j1) && (abs(m2) <= j2) && (abs(M) <= J)

cg_l_condition(j1, j2, J) = (abs(j1-j2) <= J <= j1 + j2)

cg_m_condition(m1, m2, M) = (M == m1 + m2)

"""
`clebschgordan(j1, m1, j2, m2, J, M, T=Float64)` :

A reference implementation of Clebsch-Gordon coefficients based on

https://hal.inria.fr/hal-01851097/document
Equation (4-6)

This heavily uses BigInt and BigFloat and should therefore not be employed
for performance critical tasks, but only precomputation.

The ordering of parameters corresponds to the following convention:
```
clebschgordan(j1, m1, j2, m2, J, M) = C_{j1m1j2m2}^{JM}
```
where
```
   D_{m1k1}^{l1} D_{m2k2}^{l2}}
	=
	∑_j  C_{l1m1l2m2}^{j(m1+m2)} C_{l1k1l2k2}^{j2(k1+k2)} D_{(m1+m2)(k1+k2)}^{j}
```
"""
function clebschgordan(j1, m1, j2, m2, J, M, T=Float64)
	if !cg_conditions(j1, m1, j2, m2, J, M)
		return zero(T)
	end

   N = (2*J+1) *
       factorial(big(j1+m1)) * factorial(big(j1-m1)) *
       factorial(big(j2+m2)) * factorial(big(j2-m2)) *
       factorial(big(J+M)) * factorial(big(J-M)) /
       factorial(big( j1+j2-J)) /
       factorial(big( j1-j2+J)) /
       factorial(big(-j1+j2+J)) /
       factorial(big(j1+j2+J+1))

   G = big(0)
   # 0 ≦ k ≦ j1+j2-J
   # 0 ≤ j1-m1-k ≤ j1-j2+J   <=>   j2-J-m1 ≤ k ≤ j1-m1
   # 0 ≤ j2+m2-k ≤ -j1+j2+J  <=>   j1-J+m2 ≤ k ≤ j2+m2
   lb = (0, j2-J-m1, j1-J+m2)
   ub = (j1+j2-J, j1-m1, j2+m2)
   for k in maximum(lb):minimum(ub)
      bk = big(k)
      G += (-1)^k *
           binomial(big( j1+j2-J), big(k)) *
           binomial(big( j1-j2+J), big(j1-m1-k)) *
           binomial(big(-j1+j2+J), big(j2+m2-k))
   end

   return T(sqrt(N) * G)
end

ClebschGordan(T=Float64) =
	ClebschGordan{T}(Dict{Tuple{Int,Int,Int,Int,Int,Int}, T}())

_cg_key(j1, m1, j2, m2, J, M) = (j1, m1, j2, m2, J, M)

function (cg::ClebschGordan{T})(j1, m1, j2, m2, J, M) where {T}
	if !cg_conditions(j1,m1, j2,m2, J,M)
		return zero(T)
	end
	key = _cg_key(j1, m1, j2, m2, J, M)
	if haskey(cg.vals, key)
		return cg.vals[key]
	end
	val = clebschgordan(j1, m1, j2, m2, J, M, T)
	cg.vals[key] = val
	return val
end

# ----------------------------------------------------------------------
#     Rot3DCoeffs code: generalized cg coefficients
#
#  Note: in this section kk is a tuple of m-values, it is not
#        related to the k index in the 1-p basis (or radial basis)
# ----------------------------------------------------------------------

dicttype(N::Integer, TP) = dicttype(Val(N), TP)

dicttype(::Val{N}, TP) where {N} =
   Dict{Tuple{SVector{N,Int}, SVector{N,Int}, SVector{N,Int}}, TP}

Rot3DCoeffs(φ, T=Float64) = Rot3DCoeffs(Dict[], ClebschGordan(T), φ)

function get_vals(A::Rot3DCoeffs{T}, valN::Val{N}) where {T,N}
	# make up an ll, kk, mm and compute a dummy coupling coeff
	ll, mm, kk = SVector(0), SVector(0), SVector(0)
	cc0 = coco_zeros(A.phi, ll, mm, kk, T, A)
	TP = typeof(cc0)
	if length(A.vals) < N
		# create more dictionaries of the correct type
		for n = length(A.vals)+1:N
			push!(A.vals, dicttype(n, TP)())
		end
	end
   return (A.vals[N])::(dicttype(valN, TP))
end

_key(ll::StaticVector{N}, mm::StaticVector{N}, kk::StaticVector{N}) where {N} =
      (SVector{N, Int}(ll), SVector{N, Int}(mm), SVector{N, Int}(kk))

function (A::Rot3DCoeffs{T})(ll::StaticVector{N},
                             mm::StaticVector{N},
                             kk::StaticVector{N}) where {T, N}
   vals = get_vals(A, Val(N))  # this should infer the type!
   key = _key(ll, mm, kk)
   if haskey(vals, key)
      val  = vals[key]
   else
      val = _compute_val(A, key...)
      vals[key] = val
   end
   return val
end

# the recursion has two steps so we need to define the
# coupling coefficients for N = 1, 2
# TODO: actually this seems false; it is only one recursion step, and a bit
#       or reshuffling should allow us to get rid of the {N = 2} case.

(A::Rot3DCoeffs{T})(ll::StaticVector{1},
                 mm::StaticVector{1},
                 kk::StaticVector{1}) where {T} =
		coco_init(A.phi, ll[1], mm[1], kk[1], T, A)

function _compute_val(A::Rot3DCoeffs{T}, ll::StaticVector{N},
                                         mm::StaticVector{N},
                                         kk::StaticVector{N}) where {T, N}
	val = coco_zeros(A.phi, ll, mm, kk, T, A)
	TV = typeof(val)

	tmp = zero(MVector{N-1, Int})

	function _get_pp(aa, ap)
		for i = 1:N-2
			@inbounds tmp[i] = aa[i]
		end
		tmp[N-1] = ap
		return SVector(tmp)
	end

	jmin = maximum( ( abs(ll[N-1]-ll[N]),
				         abs(kk[N-1]+kk[N]),
						   abs(mm[N-1]+mm[N]) ) )
   jmax = ll[N-1]+ll[N]
   for j = jmin:jmax
		cgk = A.cg(ll[N-1], kk[N-1], ll[N], kk[N], j, kk[N-1]+kk[N])
		cgm = A.cg(ll[N-1], mm[N-1], ll[N], mm[N], j, mm[N-1]+mm[N])
		if cgk * cgm  != 0
			llpp = _get_pp(ll, j) # SVector(llp..., j)
			mmpp = _get_pp(mm, mm[N-1]+mm[N]) # SVector(mmp..., mm[N-1]+mm[N])
			kkpp = _get_pp(kk, kk[N-1]+kk[N]) # SVector(kkp..., kk[N-1]+kk[N])
			a = A(llpp, mmpp, kkpp)::TV
			val += cgk * cgm * a
		end
   end
   return val
end

# ----------------------------------------------------------------------
#   construction of a possible set of generalised CG coefficient;
#   numerically via SVD; this could be done analytically which might
#   be more efficient.
# ----------------------------------------------------------------------

function re_basis(A::Rot3DCoeffs{T}, ll::SVector) where {T}
	TP = typeof(A.phi)
	TCC = coco_type(TP)
	CC, Mll = compute_Al(A, ll)  # CC::Vector{Vector{...}}
	G = [ sum( coco_dot(CC[a][i], CC[b][i]) for i = 1:length(Mll) )
			for a = 1:length(CC), b = 1:length(CC) ]
	svdC = svd(G)
	rk = rank(Diagonal(svdC.S), rtol = 1e-7)
	# Diagonal(sqrt.(svdC.S[1:rk])) * svdC.U[:, 1:rk]' * CC
	# construct the new basis
	Ured = Diagonal(sqrt.(svdC.S[1:rk])) * svdC.U[:, 1:rk]'
	Ure = Matrix{TCC}(undef, rk, length(Mll))
	for i = 1:rk
		Ure[i, :] = sum(Ured[i, j] * CC[j]  for j = 1:length(CC))
	end
	return Ure, Mll
end

# function barrier
function compute_Al(A::Rot3DCoeffs, ll::SVector)
	Mll = collect(_mrange(A.phi, ll))
	TP = typeof(A.phi)
	if length(Mll) == 0
		return Vector{TP}[], Mll
	end

	TA = typeof(A(ll, Mll[1], Mll[1]))
	return __compute_Al(A, ll, Mll, TP, TA)
end

# TODO: what was TA for? Can we get rid of it via coco_type?

function __compute_Al(A::Rot3DCoeffs{T}, ll, Mll, TP, TA) where {T}
	lenMll = length(Mll)
	# each element of CC will be one row of the coupling coefficients
	TCC = coco_type(TP)
	CC = Vector{TCC}[]
	# some utility funcions to allow coco_init to return either a property
	# or a vector of properties
	function __into_cc!(cc, cc0::AbstractProperty, im)
		@assert length(cc) == 1
		cc[1][im] = cc0
	end
	function __into_cc!(cc, cc0::AbstractVector, im)
		@assert length(cc) == length(cc0)
		for p = 1:length(cc)
			cc[p][im] = cc0[p]
		end
	end

	for (ik, kk) in enumerate(Mll)  # loop over possible basis functions
		# do a dummy calculation to determine how many coefficients we will get
		cc0 = A(ll, Mll[1], kk)::TA
		numcc = (cc0 isa AbstractProperty ? 1 : length(cc0))
		# allocate the right number of vectors to store basis function coeffs
		cc = [ Vector{TCC}(undef, lenMll) for _=1:numcc ]
		for (im, mm) in enumerate(Mll) # loop over possible indices
			if !coco_filter(A.phi, ll, mm, kk)
				cc00 = zeros(TP, length(cc))::TA
				__into_cc!(cc, cc00, im)
			else
				# get all possible coupling coefficients
				cc0 = A(ll, mm, kk)::TA
				__into_cc!(cc, cc0, im)
			end
		end
		# and now push them onto the big stack.
		append!(CC, cc)
	end

	return CC, Mll
end

end


text:
# Note: construction of subviews involve several steps:
#
#  - compute view offset (w.r.t. the parent)
#  - compute view shape
#  - compute strides (for non-contiguous views only)
#  - decide contiguous rank (statically)
#  - make the view
#

##### Compute offset #####

_offset(i::Colon) = 0
_offset(i::Int) = i - 1
_offset(i::Real) = to_index(i) - 1
_offset(i::Range) = to_index(first(i)) - 1

_step(i::Real) = 1
_step(i::Colon) = 1
_step(i::Range) = step(i)

# aoffset: offset w.r.t. the underlying array (i.e. parent)

@compat aoffset(a::Union{Array, UnsafeArrayView}, i::Subs) = roffset(a, i)
@compat aoffset(a::Union{Array, UnsafeArrayView}, i1::Subs, i2::Subs) = roffset(a, i1, i2)
@compat aoffset(a::Union{Array, UnsafeArrayView}, i1::Subs, i2::Subs, i3::Subs) = roffset(a, i1, i2, i3)
@compat aoffset(a::Union{Array, UnsafeArrayView}, i1::Subs, i2::Subs, i3::Subs, i4::Subs) =
    roffset(a, i1, i2, i3, i4)
@compat aoffset(a::Union{Array, UnsafeArrayView}, i1::Subs, i2::Subs, i3::Subs, i4::Subs, i5::Subs, I::Subs...) =
    roffset(a, i1, i2, i3, i4, i5, I...)

aoffset(a::ArrayView, i::Subs) = a.offset + roffset(a, i)
aoffset(a::ArrayView, i1::Subs, i2::Subs) = a.offset + roffset(a, i1, i2)
aoffset(a::ArrayView, i1::Subs, i2::Subs, i3::Subs) =
    a.offset + roffset(a, i1, i2, i3)
aoffset(a::ArrayView, i1::Subs, i2::Subs, i3::Subs, i4::Subs) =
    a.offset + roffset(a, i1, i2, i3, i4)
aoffset(a::ArrayView, i1::Subs, i2::Subs, i3::Subs, i4::Subs, i5::Subs, I::Subs...) =
    a.offset + roffset(a, i1, i2, i3, i4, i5, I...)

# roffset: offset w.r.t. the first element of the view

# 1D

roffset(a::ContiguousArray, i::Colon) = 0
roffset(a::ContiguousArray, i::SubsNC) = _offset(i)

roffset(a::StridedVector, i::Colon) = 0
roffset(a::StridedVector, i::SubsNC) = _offset(i) * stride(a,1)

# 2D

roffset(a::ContiguousArray, i1::Colon, i2::Colon) = 0
roffset(a::ContiguousArray, i1::Colon, i2::SubsNC) = size(a,1) * _offset(i2)
roffset(a::ContiguousArray, i1::SubsNC, i2::Colon) = _offset(i1)
roffset(a::ContiguousArray, i1::SubsNC, i2::SubsNC) =
    _offset(i1) + size(a,1) * _offset(i2)

roffset(a::StridedMatrix, i1::Colon, i2::Colon) = 0
roffset(a::StridedMatrix, i1::Colon, i2::SubsNC) = _offset(i2) * stride(a,2)
roffset(a::StridedMatrix, i1::SubsNC, i2::Colon) = _offset(i1) * stride(a,1)
roffset(a::StridedMatrix, i1::SubsNC, i2::SubsNC) =
    _offset(i1) * stride(a,1) + _offset(i2) * stride(a,2)

# 3D

roffset(a::ContiguousArray, i1::Colon, i2::Colon, i3::Colon) = 0
roffset(a::ContiguousArray, i1::Colon, i2::Colon, i3::SubsNC) =
    size(a,1) * size(a,2) * _offset(i3)
roffset(a::ContiguousArray, i1::Colon, i2::SubsNC, i3::Colon) =
    size(a,1) * _offset(i2)
roffset(a::ContiguousArray, i1::Colon, i2::SubsNC, i3::SubsNC) =
    size(a,1) * (_offset(i2) + size(a,2) * _offset(i3))
roffset(a::ContiguousArray, i1::SubsNC, i2::Colon, i3::Colon) = _offset(i1)
roffset(a::ContiguousArray, i1::SubsNC, i2::Colon, i3::SubsNC) =
    _offset(i1) + (size(a,1) * size(a,2) * _offset(i3))
roffset(a::ContiguousArray, i1::SubsNC, i2::SubsNC, i3::Colon) =
    _offset(i1) + size(a,1) * _offset(i2)
roffset(a::ContiguousArray, i1::SubsNC, i2::SubsNC, i3::SubsNC) =
    _offset(i1) + size(a,1) * (_offset(i2) + size(a,2) * _offset(i3))

roffset{T}(a::StridedArray{T,3}, i1::Colon, i2::Colon, i3::Colon) = 0
roffset{T}(a::StridedArray{T,3}, i1::Colon, i2::Colon, i3::SubsNC) =
    _offset(i3) * stride(a,3)
roffset{T}(a::StridedArray{T,3}, i1::Colon, i2::SubsNC, i3::Colon) =
    _offset(i2) * stride(a,2)
roffset{T}(a::StridedArray{T,3}, i1::Colon, i2::SubsNC, i3::SubsNC) =
    _offset(i2) * stride(a,2) + _offset(i3) * stride(a,3)
roffset{T}(a::StridedArray{T,3}, i1::SubsNC, i2::Colon, i3::Colon) =
    _offset(i1) * stride(a,1)
roffset{T}(a::StridedArray{T,3}, i1::SubsNC, i2::Colon, i3::SubsNC) =
    _offset(i1) * stride(a,1) + _offset(i3) * stride(a,3)
roffset{T}(a::StridedArray{T,3}, i1::SubsNC, i2::SubsNC, i3::Colon) =
    _offset(i1) * stride(a,1) + _offset(i2) * stride(a,2)
roffset{T}(a::StridedArray{T,3}, i1::SubsNC, i2::SubsNC, i3::SubsNC) =
    _offset(i1) * stride(a,1) + _offset(i2) * stride(a,2) + _offset(i3) * stride(a,3)

# 4D (partial)
function roffset(a::ContiguousArray, i1::Subs, i2::Subs, i3::Subs, i4::Subs)
    o = _offset(i1)
    s = size(a,1)
    o += s * _offset(i2)
    o += (s *= size(a,2)) * _offset(i3)
    o += (s *= size(a,3)) * _offset(i4)
    return o
end

function roffset(a::ContiguousArray, i1::Colon, i2::Subs, i3::Subs, i4::Subs)
    s = size(a,1)
    o = s * _offset(i2)
    o += (s *= size(a,2)) * _offset(i3)
    o += (s *= size(a,3)) * _offset(i4)
    return o
end

function roffset(a::ContiguousArray, i1::Colon, i2::Colon, i3::Subs, i4::Subs)
    s = size(a,1)
    o = (s *= size(a,2)) * _offset(i3)
    o += (s *= size(a,3)) * _offset(i4)
    return o
end

# General

roffset(a::ContiguousArray, i1::Colon, i2::Colon, i3::Colon, i4::Colon, I::Colon...) = 0

function roffset(a::ContiguousArray, i1::Subs, i2::Subs, i3::Subs, i4::Subs, I::Subs...)
    o = _offset(i1)
    s = size(a,1)
    o += s * _offset(i2)
    o += (s *= size(a,2)) * _offset(i3)
    o += (s *= size(a,3)) * _offset(i4)
    for i = 1:length(I)
        o += (s *= size(a,i+3)) * _offset(I[i])
    end
    return o
end

roffset(a::StridedArray, i1::Subs, i2::Subs, i3::Subs, I::Subs...) =
    _roffset(strides(a), tuple(i1, i2, i3, I...))::Int

function _roffset{N}(ss::NTuple{N,Int}, subs::NTuple{N})
    o = _offset(subs[1]) * ss[1]
    for i = 2:N
        o += _offset(subs[i]) * ss[i]
    end
    return o
end

##### Compute view shape #####

_dim(a::AbstractArray, d::Int, r::Colon) = size(a, d)
_dim(a::AbstractArray, d::Int, r::Range) = length(r)
_dim(a::AbstractArray, d::Int, r::Real) = 1

_dim{N}(siz::NTuple{N,Int}, d::Int, r::Colon) = d <= N ? siz[d] : 1
_dim(siz::Tuple, d::Int, r::Range) = length(r)
_dim(siz::Tuple, d::Int, r::Real) = 1

# 1D
vshape(a::DenseArray, i::Real) = ()
vshape(a::DenseArray, i::Colon) = (length(a),)
vshape(a::DenseArray, i::Range) = (length(i),)

# 2D

_succlen2{T}(a::DenseArray{T,1}) = 1
_succlen2{T}(a::DenseArray{T,2}) = size(a, 2)
_succlen2{T}(a::DenseArray{T,3}) = size(a, 2) * size(a, 3)
_succlen2(a::DenseArray) = prod(size(a)[2:end])::Int

vshape(a::DenseArray, i1::Real, i2::Real) = ()
vshape(a::DenseArray, i1::SubsRange, i2::Real) = (_dim(a,1,i1),)
vshape(a::DenseArray, i1::Subs, i2::Colon) = (_dim(a,1,i1), _succlen2(a))
vshape(a::DenseArray, i1::Subs, i2::Range) = (_dim(a,1,i1), length(i2))

# 3D

_succlen3{T}(a::DenseArray{T,1}) = 1
_succlen3{T}(a::DenseArray{T,2}) = 1
_succlen3{T}(a::DenseArray{T,3}) = size(a, 3)
_succlen3{T}(a::DenseArray{T,4}) = size(a, 3) * size(a, 4)
_succlen3(a::DenseArray) = prod(size(a)[3:end])::Int

vshape(a::DenseArray, i1::Real, i2::Real, i3::Real) = ()
vshape(a::DenseArray, i1::SubsRange, i2::Real, i3::Real) = (_dim(a,1,i1),)
vshape(a::DenseArray, i1::Subs, i2::SubsRange, i3::Real) =
    (_dim(a,1,i1), _dim(a,2,i2))

vshape(a::DenseArray, i1::Subs, i2::Subs, i3::Colon) =
    (_dim(a,1,i1), _dim(a,2,i2), _succlen3(a))
vshape(a::DenseArray, i1::Subs, i2::Subs, i3::Range) =
    (_dim(a,1,i1), _dim(a,2,i2), length(i3))

# 4D

_succlen4{T}(a::DenseArray{T,1}) = 1
_succlen4{T}(a::DenseArray{T,2}) = 1
_succlen4{T}(a::DenseArray{T,3}) = 1
_succlen4{T}(a::DenseArray{T,4}) = size(a, 4)
_succlen4{T}(a::DenseArray{T,5}) = size(a, 4) * size(a, 5)
_succlen4(a::DenseArray) = prod(size(a)[4:end])::Int

vshape(a::DenseArray, i1::Real, i2::Real, i3::Real, i4::Real) = ()
vshape(a::DenseArray, i1::SubsRange, i2::Real, i3::Real, i4::Real) = (_dim(a,1,i1),)
vshape(a::DenseArray, i1::Subs, i2::SubsRange, i3::Real, i4::Real) =
    (_dim(a,1,i1), _dim(a,2,i2))
vshape(a::DenseArray, i1::Subs, i2::Subs, i3::SubsRange, i4::Real) =
    (_dim(a,1,i1), _dim(a,2,i2), _dim(a,3,i3))

vshape(a::DenseArray, i1::Subs, i2::Subs, i3::Subs, i4::Colon) =
    (_dim(a,1,i1), _dim(a,2,i2), _dim(a,3,i3), _succlen4(a))
vshape(a::DenseArray, i1::Subs, i2::Subs, i3::Subs, i4::Range) =
    (_dim(a,1,i1), _dim(a,2,i2), _dim(a,3,i3), length(i4))

# multi-dimensional

vshape{T,N}(a::DenseArray{T,N}, i1::Subs, i2::Subs, i3::Subs, i4::Subs, i5::Subs, I::Subs...) =
    _vshape(size(a), i1, i2, i3, i4, i5, I...)

_vshape{N}(siz::NTuple{N,Int}, i1::Real) = ()
_vshape{N}(siz::NTuple{N,Int}, i1::Real, i2::Real...) = ()

_vshape{N}(siz::NTuple{N,Int}, i1::Colon) = (prod(siz),)
_vshape{N}(siz::NTuple{N,Int}, i1::Range) = (length(i1),)
_vshape{N}(siz::NTuple{N,Int}, i1::Subs, i2::Subs...) = tuple(_dim(siz,1,i1), _vshape(siz[2:N], i2...)...)

##### Compute strides #####

# 1D

vstrides(a::ContiguousArray, i::Subs) = (_step(i),)
vstrides(a::DenseArray, i::Subs) = (stride(a,1) * _step(i),)

# 2D

vstrides(a::ContiguousArray, i1::Subs, i2::Real) = (_step(i1),)
@compat vstrides(a::ContiguousArray, i1::Subs, i2::Union{Colon,UnitRange}) = (_step(i1), stride(a,2))
vstrides(a::ContiguousArray, i1::Subs, i2::Range) = (_step(i1), stride(a,2) * _step(i2))

vstrides(a::DenseArray, i1::Subs, i2::Real) = (stride(a,1) * _step(i1),)
vstrides(a::DenseArray, i1::Subs, i2::Subs) =
    (stride(a,1) * _step(i1), stride(a,2) * _step(i2))

# 3D

vstrides(a::ContiguousArray, i1::Subs, i2::Real, i3::Real) =
    (_step(i1),)
vstrides(a::ContiguousArray, i1::Subs, i2::Subs, i3::Real) =
    (_step(i1), stride(a,2) * _step(i2))
vstrides(a::ContiguousArray, i1::Subs, i2::Subs, i3::Subs) =
    (_step(i1), stride(a,2) * _step(i2), stride(a,3) * _step(i3))

vstrides(a::DenseArray, i1::Subs, i2::Real, i3::Real) =
    (stride(a,1) * _step(i1),)
vstrides(a::DenseArray, i1::Subs, i2::Subs, i3::Real) =
    (stride(a,1) * _step(i1), stride(a,2) * _step(i2))
vstrides(a::DenseArray, i1::Subs, i2::Subs, i3::Subs) =
    (stride(a,1) * _step(i1), stride(a,2) * _step(i2), stride(a,3) * _step(i3))

# multi-dimensional array

vstrides(a::DenseArray, i1::Subs, i2::Subs, i3::Subs, i4::Subs, I::Subs...) =
    _vstrides(strides(a), 1, i1, i2, i3, i4, I...)

_vstrides{N}(ss::NTuple{N,Int}, k::Int, i1::Real, i2::Real) = ()
_vstrides{N}(ss::NTuple{N,Int}, k::Int, i1::Subs, i2::Real) =
    (ss[k] * _step(i1),)
_vstrides{N}(ss::NTuple{N,Int}, k::Int, i1::Subs, i2::Subs) =
    (ss[k] * _step(i1), ss[k+1] * _step(i2))

_vstrides{N}(ss::NTuple{N,Int}, k::Int, i1::Real, i2::Real, i3::Real, I::Real...) = ()
_vstrides{N}(ss::NTuple{N,Int}, k::Int, i1::Subs, i2::Subs, i3::Subs, I::Subs...) =
    tuple(ss[k] * _step(i1), _vstrides(ss, k+1, i2, i3, I...)...)

##### View construction ######

make_view{N}(a::DenseArray, cr::Type{ContRank{N}}, shp::NTuple{N,Int}, i::Subs) =
    ContiguousView(parent(a), aoffset(a, i), shp)

make_view{N}(a::DenseArray, cr::Type{ContRank{N}}, shp::NTuple{N,Int}, i1::Subs, i2::Subs) =
    ContiguousView(parent(a), aoffset(a, i1, i2), shp)

make_view{N}(a::DenseArray, cr::Type{ContRank{N}}, shp::NTuple{N,Int}, i1::Subs, i2::Subs, i3::Subs) =
    ContiguousView(parent(a), aoffset(a, i1, i2, i3), shp)

make_view{N}(a::DenseArray, cr::Type{ContRank{N}}, shp::NTuple{N,Int}, i1::Subs, i2::Subs, i3::Subs, i4::Subs) =
    ContiguousView(parent(a), aoffset(a, i1, i2, i3, i4), shp)

make_view{N}(a::DenseArray, cr::Type{ContRank{N}}, shp::NTuple{N,Int}, i1::Subs, i2::Subs, i3::Subs, i4::Subs, i5::Subs, I::Subs...) =
    ContiguousView(parent(a), aoffset(a, i1, i2, i3, i4, i5, I...), shp)

make_view{M,N}(a::DenseArray, cr::Type{ContRank{M}}, shp::NTuple{N,Int}, i::Subs) =
    StridedView(parent(a), aoffset(a, i), shp, cr, vstrides(a, i))

make_view{M,N}(a::DenseArray, cr::Type{ContRank{M}}, shp::NTuple{N,Int}, i1::Subs, i2::Subs) =
    StridedView(parent(a), aoffset(a, i1, i2), shp, cr, vstrides(a, i1, i2))

make_view{M,N}(a::DenseArray, cr::Type{ContRank{M}}, shp::NTuple{N,Int}, i1::Subs, i2::Subs, i3::Subs) =
    StridedView(parent(a), aoffset(a, i1, i2, i3), shp, cr, vstrides(a, i1, i2, i3))

make_view{M,N}(a::DenseArray, cr::Type{ContRank{M}}, shp::NTuple{N,Int}, i1::Subs, i2::Subs, i3::Subs, i4::Subs) =
    StridedView(parent(a), aoffset(a, i1, i2, i3, i4), shp, cr, vstrides(a, i1, i2, i3, i4))

make_view{M,N}(a::DenseArray, cr::Type{ContRank{M}}, shp::NTuple{N,Int}, i1::Subs, i2::Subs, i3::Subs, i4::Subs, i5::Subs, I::Subs...) =
    StridedView(parent(a), aoffset(a, i1, i2, i3, i4, i5, I...), shp, cr, vstrides(a, i1, i2, i3, i4, i5, I...))

make_unsafe_view{N}(a::DenseArray, cr::Type{ContRank{N}}, shp::NTuple{N,Int}, i::Subs) =
    UnsafeContiguousView(parent_or_ptr(a), aoffset(a, i), shp)

make_unsafe_view{N}(a::DenseArray, cr::Type{ContRank{N}}, shp::NTuple{N,Int}, i1::Subs, i2::Subs) =
    UnsafeContiguousView(parent_or_ptr(a), aoffset(a, i1, i2), shp)

make_unsafe_view{N}(a::DenseArray, cr::Type{ContRank{N}}, shp::NTuple{N,Int}, i1::Subs, i2::Subs, i3::Subs) =
    UnsafeContiguousView(parent_or_ptr(a), aoffset(a, i1, i2, i3), shp)

make_unsafe_view{N}(a::DenseArray, cr::Type{ContRank{N}}, shp::NTuple{N,Int}, i1::Subs, i2::Subs, i3::Subs, i4::Subs) =
    UnsafeContiguousView(parent_or_ptr(a), aoffset(a, i1, i2, i3, i4), shp)

make_unsafe_view{N}(a::DenseArray, cr::Type{ContRank{N}}, shp::NTuple{N,Int}, i1::Subs, i2::Subs, i3::Subs, i4::Subs, i5::Subs, I::Subs...) =
    UnsafeContiguousView(parent_or_ptr(a), aoffset(a, i1, i2, i3, i4, i5, I...), shp)

make_unsafe_view{M,N}(a::DenseArray, cr::Type{ContRank{M}}, shp::NTuple{N,Int}, i::Subs) =
    UnsafeStridedView(parent_or_ptr(a), aoffset(a, i), shp, cr, vstrides(a, i))

make_unsafe_view{M,N}(a::DenseArray, cr::Type{ContRank{M}}, shp::NTuple{N,Int}, i1::Subs, i2::Subs) =
    UnsafeStridedView(parent_or_ptr(a), aoffset(a, i1, i2), shp, cr, vstrides(a, i1, i2))

make_unsafe_view{M,N}(a::DenseArray, cr::Type{ContRank{M}}, shp::NTuple{N,Int}, i1::Subs, i2::Subs, i3::Subs) =
    UnsafeStridedView(parent_or_ptr(a), aoffset(a, i1, i2, i3), shp, cr, vstrides(a, i1, i2, i3))

make_unsafe_view{M,N}(a::DenseArray, cr::Type{ContRank{M}}, shp::NTuple{N,Int}, i1::Subs, i2::Subs, i3::Subs, i4::Subs) =
    UnsafeStridedView(parent_or_ptr(a), aoffset(a, i1, i2, i3, i4), shp, cr, vstrides(a, i1, i2, i3, i4))

make_unsafe_view{M,N}(a::DenseArray, cr::Type{ContRank{M}}, shp::NTuple{N,Int}, i1::Subs, i2::Subs, i3::Subs, i4::Subs, i5::Subs, I::Subs...) =
    UnsafeStridedView(parent_or_ptr(a), aoffset(a, i1, i2, i3, i4, i5, I...), shp, cr, vstrides(a, i1, i2, i3, i4, i5, I...))

##### Interface

view(a::Array) = ContiguousView(a, size(a))
view(a::ArrayView) = a

view(a::DenseArray, i::Subs) =
    (shp = vshape(a, i); make_view(a, restrict_crank(acontrank(a, i), shp), shp, i))

view(a::DenseArray, i1::Subs, i2::Subs) =
    (shp = vshape(a, i1, i2); make_view(a, restrict_crank(acontrank(a, i1, i2), shp), shp, i1, i2))

view(a::DenseArray, i1::Subs, i2::Subs, i3::Subs) =
    (shp = vshape(a, i1, i2, i3); make_view(a, restrict_crank(acontrank(a, i1, i2, i3), shp), shp, i1, i2, i3))

view(a::DenseArray, i1::Subs, i2::Subs, i3::Subs, i4::Subs) =
    (shp = vshape(a, i1, i2, i3, i4); make_view(a, restrict_crank(acontrank(a, i1, i2, i3, i4), shp), shp, i1, i2, i3, i4))

view(a::DenseArray, i1::Subs, i2::Subs, i3::Subs, i4::Subs, i5::Subs, I::Subs...) =
    (shp = vshape(a, i1, i2, i3, i4, i5, I...);
     make_view(a, restrict_crank(acontrank(a, i1, i2, i3, i4, i5, I...), shp), shp, i1, i2, i3, i4, i5, I...))

unsafe_view(a::Array) = ContiguousView(a, size(a))
unsafe_view(a::ArrayView) = a

unsafe_view(a::DenseArray, i::Subs) =
    (shp = vshape(a, i); make_unsafe_view(a, restrict_crank(acontrank(a, i), shp), shp, i))

unsafe_view(a::DenseArray, i1::Subs, i2::Subs) =
    (shp = vshape(a, i1, i2); make_unsafe_view(a, restrict_crank(acontrank(a, i1, i2), shp), shp, i1, i2))

unsafe_view(a::DenseArray, i1::Subs, i2::Subs, i3::Subs) =
    (shp = vshape(a, i1, i2, i3); make_unsafe_view(a, restrict_crank(acontrank(a, i1, i2, i3), shp), shp, i1, i2, i3))

unsafe_view(a::DenseArray, i1::Subs, i2::Subs, i3::Subs, i4::Subs) =
    (shp = vshape(a, i1, i2, i3, i4); make_unsafe_view(a, restrict_crank(acontrank(a, i1, i2, i3, i4), shp), shp, i1, i2, i3, i4))

unsafe_view(a::DenseArray, i1::Subs, i2::Subs, i3::Subs, i4::Subs, i5::Subs, I::Subs...) =
    (shp = vshape(a, i1, i2, i3, i4, i5, I...);
     make_unsafe_view(a, restrict_crank(acontrank(a, i1, i2, i3, i4, i5, I...), shp), shp, i1, i2, i3, i4, i5, I...))


text:
module TheAlgorithms

# Usings/Imports (keep sorted)
using LinearAlgebra
using Random

## Exports
# Please keep the folders/functions sorted

# Exports: data_structures
export AbstractBinaryTree
export AbstractBinaryTree_arr
export BinaryTree
export ch
export depth
export DisjointSet
export find
export height
export insert!
export isleaf
export left
export merge!
export right

# Exports: knapsack
export complete_pack!
export zero_one_pack!

# Exports: math
export SIR # TODO: make the name lowercase if possible
export abs_max
export abs_min
export abs_val
export area_circle
export area_ellipse
export area_heron_triangle
export area_parallelogram
export area_rectangle
export area_rhombus
export area_square
export area_trapezium
export area_triangle
export ceil_val
export collatz_sequence
export eratosthenes
export euler_method
export factorial_iterative
export factorial_recursive
export floor_val
export is_armstrong
export line_length
export krishnamurthy
export mean
export median
export mode
export monte_carlo_integration
export prime_check
export prime_factors
export perfect_cube
export perfect_number
export perfect_square
export prime_check
export prime_factors
export sum_ap
export sum_gp
export surfarea_cube
export surfarea_sphere
export trapazoidal_area
export verlet_integration
export vol_cube
export vol_cuboid
export vol_cone
export vol_right_circ_cone
export vol_prism
export vol_pyramid
export vol_sphere
export vol_circular_cylinder

# Exports: matrix
export determinant
export lu_decompose
export rotation_matrix

# Exports: project-rosalind
export count_nucleotides
export dna2rna
export reverse_complement

# Exports: searches
export binary_search
export exponential_search
export interpolation_search
export jump_search
export linear_search

# Exports: sorts
export BubbleSort!,InsertionSort!,SelectionSort!

# Exports: statistics
export OLSbeta # TODO: make the name lowercase if possible
export pearson_correlation
export variance

# Exports: strings
export is_palindrome
export detect_anagrams

# Exports: scheduling
export fcfs

# Exports: conversions
export celsius_to_fahrenheit
export celsius_to_kelvin
export fahrenheit_to_celsius
export fahrenheit_to_kelvin
export kelvin_to_celsius
export kelvin_to_fahrenheit
export weight_conversion

## Includes
# Please keep the folders/files sorted (by dependencies then alphabetical order)

# Includes: data_structures
include("data_structures/binary_tree/basic_binary_tree.jl")
include("data_structures/disjoint_set/disjoint_set.jl")

# Includes: knapsack
include("knapsack/knapsack.jl")

# Includes: math
include("math/abs.jl")
include("math/area.jl")
include("math/area_under_curve.jl")
include("math/armstrong_number.jl")
include("math/average_mean.jl")
include("math/average_median.jl")
include("math/average_mode.jl")
include("math/ceil_floor.jl") # needed by average_median
include("math/collatz_sequence.jl")
include("math/euler_method.jl")
include("math/factorial.jl")
include("math/krishnamurthy_number.jl")
include("math/line_length.jl")
include("math/krishnamurthy_number.jl")
include("math/monte_carlo_integration.jl")
include("math/prime_check.jl")
include("math/prime_factors.jl")
include("math/perfect_cube.jl")
include("math/perfect_number.jl")
include("math/perfect_square.jl")
include("math/prime_check.jl")
include("math/prime_factors.jl")
include("math/sieve_of_eratosthenes.jl")
include("math/sir_model.jl")
include("math/sum_of_arithmetic_series.jl")
include("math/sum_of_geometric_progression.jl")
include("math/verlet.jl")
include("math/volume.jl")

# Includes: matrix
include("matrix/lu_decompose.jl") # used by determinant.jl
include("matrix/determinant.jl")
include("matrix/rotation-matrix.jl")

# Includes: project-rosalind
include("project-rosalind/count_nucleotide.jl")
include("project-rosalind/dna2rna.jl")
include("project-rosalind/reverse_complement.jl")

# Includes: scheduling
include("scheduling/fcfs.jl")

# Includes: search
include("searches/binary_search.jl")
include("searches/exponential_search.jl")
include("searches/interpolation_search.jl")
include("searches/jump_search.jl")
include("searches/linear_search.jl")

# Includes: sorts
include("sorts/bubble_sort.jl")
include("sorts/insertion_sort.jl")
include("sorts/selection_sort.jl")

# Includes: statistics
include("statistics/ordinary_least_squares.jl")
include("statistics/pearson_correlation.jl")
include("statistics/variance.jl")

# Includes: strings
include("strings/is_palindrome.jl")
include("strings/detect_anagrams.jl")

# Includes: scheduling
include("scheduling/fcfs.jl")

# Includes: conversions
include("conversions/weight_conversion.jl")
include("conversions/temparature_conversion.jl")

end

max_stars_repo_name: mayurdahibhate/Julia

max_issues_repo_name: mayurdahibhate/Julia

max_forks_repo_name: mayurdahibhate/Julia


text:
using DataStructures
using DataStructures: OrderedDict, Deque
using Base

mutable struct PriceLevel
    orders::OrderedDict{UInt128, LimitOrder}
    size::UInt128
    function PriceLevel()
        new(OrderedDict{UInt128, LimitOrder}(), UInt128(0))
    end
end

function orders(pl::PriceLevel)::Array{LimitOrder}
    return pl.orders.vals
end

function Base.size(pl::PriceLevel)::UInt128
    return pl.size
end

function add!(pl::PriceLevel, o::LimitOrder)
    pl.orders[id(o)] = o
    pl.size += size(o)
end

function remove!(pl::PriceLevel, order_id::UInt128)
    o = pop!(pl.orders, order_id)
    pl.size -= size(o)

    if pl.size == UInt128(0)
        pl.orders = OrderedDict{UInt128, LimitOrder}()
    end
end

function DataStructures.first(pl::PriceLevel)
    return DataStructures.first(pl.orders).second
end

function match!(pl::PriceLevel, size::UInt128, trades::Deque{LimitOrder})::Tuple{Deque{LimitOrder}, UInt128}
    # Matches the size to the price level
    #for order in values(pl.orders)
    while pl.size > 0
        order = first(values(pl.orders))
        diff = Int128(order.size) - Int128(size)
        # Current limit order is larger than ordered size
        if diff > 0
            order.size = order.size - size
            pl.size = pl.size - size
            push!(trades, LimitOrder(order.price, size, order.trader_id,
                                     order.side, order.order_id))
            return trades, 0
        # Current limit order is smaller or equal to than ordered size
        elseif diff <= 0
            push!(trades, order)
            size = size - order.size
            pl.size = pl.size - order.size
            pop!(pl.orders)
            if (pl.size == 0) | (diff == 0)
                return trades, size
            end
        end
    end
end

function amount(trades::Deque{LimitOrder})::UInt128
    a::UInt128 = 0
    for order in trades
        a += order.size*order.price
    end
    return a
end


text:
#based on the optimizing angle tutorial
@testset "optimize angle" begin
    λ0 = 1 #doesn't matter since everything is normalized to λ0
    k0 = 2π/λ0
    kin = 0.5k0
    pw = PlaneWave(0.0)
    M = 20

    sfun(N) = rounded_star(0.35λ0, 0.1λ0, 4, N)
    # N, errN = minimumN(k0, kin, sfun; tol=1e-6, N_start = 200, N_min = 100,
                # N_max = 400)
    shapes = [sfun(202)]
    P = 12#P, errP = minimumP(k0, kin, shapes[1]; P_max = 30, tol = 1e-7)

    centers =  rect_grid(2, div(M,2), λ0, λ0) #2xM/2 grid with distance λ0
    ids = ones(Int, M)
    φs0 = zeros(M)
    sp = ScatteringProblem(shapes, ids, centers, φs0)

    points = [-0.05 0.0; 0.05 0.0; 0.0 0.05; 0.0 -0.05]
    @assert verify_min_distance(sp, points)

    fmm_options = FMMoptions(true, acc = 6, dx = 2λ0)

    optim_options = Optim.Options(f_tol = 1e-6, iterations = 50,
                        store_trace = true, show_trace = false)
    optim_method = Optim.BFGS(;linesearch = LineSearches.BackTracking())
    res_min = optimize_φ(φs0, points, P, pw, k0, kin, shapes, centers, ids, fmm_options,
            optim_options, optim_method; minimize = true)
    @test res_min.f_converged
    res_max = optimize_φ(φs0, points, P, pw, k0, kin, shapes, centers, ids, fmm_options,
            optim_options, optim_method; minimize = false)
    @test res_max.f_converged
    sp_min = ScatteringProblem(shapes, ids, centers, res_min.minimizer)
    sp_max = ScatteringProblem(shapes, ids, centers, res_max.minimizer)

    u_bef = calc_near_field(k0, kin, P, sp, points, pw;
            opt = fmm_options, verbose = false)
    u_min = calc_near_field(k0, kin, P, sp_min, points, pw; opt = fmm_options,
            verbose = false)
    u_max = calc_near_field(k0, kin, P, sp_max, points, pw; opt = fmm_options,
            verbose = false)

    fobj0 = sum(abs2.(u_bef))
    fobj_min = sum(abs2.(u_min))
    fobj_max = sum(abs2.(u_max))
    @test fobj_max > fobj0 && fobj_min < fobj0
    @test isapprox(res_min.trace[1].value, fobj0)
    @test isapprox(res_min.minimum, fobj_min)
    @test isapprox(-res_max.minimum, fobj_max)

    border = find_border(sp, points)
end


text:
@testset "HuffmanTree and allpaths function" begin

    using DataStructures

    vocab_hash = Dict("this" => 1,
                "is" => 2,
                "the" => 3,
                "most" => 4,
                "WICKED" => 5,
                "AWESOME" =>6,
                "test" => 7,
                "evar!" => 8)

    words = ["this", "is", "the", "most", "WICKED", "AWESOME", "test", "evar!"]

    freqs = [128,64,32,16,8,4,2,1]

    pq = PriorityQueue(Dict(zip(words, freqs)))
    ht = HuffmanTree(pq, vocab_hash)

    test_nodes, test_branches = allpaths(ht, vocab_hash, words)

    true_nodes = [[7], [6, 7], [5, 6, 7], [4, 5, 6, 7], [3, 4, 5, 6, 7],
                  [2, 3, 4, 5, 6, 7], [1, 2, 3, 4, 5, 6, 7],
                  [1, 2, 3, 4, 5, 6, 7]]

    true_branches = [[0], [0, 1], [0, 1, 1], [0, 1, 1, 1], [0, 1, 1, 1, 1],
                     [0, 1, 1, 1, 1, 1], [0, 1, 1, 1, 1, 1, 1],
                     [1, 1, 1, 1, 1, 1, 1]]

   @test test_nodes == true_nodes
   @test test_branches == true_branches
end


text:
# This file is a part of TextSearch.jl

module TextSearch

import Base: broadcastable
import StatsBase: fit, predict
using SimilaritySearch, InvertedFiles, LinearAlgebra, SparseArrays
using SimilaritySearch: getminbatch
using Polyester

include("dvec.jl")

const BOW = DVEC{UInt32,Int32}

struct TextSearchBuffer
    normtext::Vector{Char}
    tokens::Vector{String}
    unigrams::Vector{String}
    io::IOBuffer
    bow::BOW

    function TextSearchBuffer(n=128)
        normtext = Char[]
        tokens = UInt64[]
        unigrams = String[]
        bow = BOW()
        sizehint!(normtext, n)
        sizehint!(tokens, n)
        sizehint!(unigrams, n)
        sizehint!(bow, n)

        new(normtext, tokens, unigrams, IOBuffer(), bow)
    end
end

const CACHES = Channel{TextSearchBuffer}(Inf)

function Base.empty!(buff::TextSearchBuffer)
    empty!(buff.normtext)
    empty!(buff.tokens)
    empty!(buff.unigrams)
    empty!(buff.bow)
end

function __init__()
    for _ in 1:Threads.nthreads()
        put!(CACHES, TextSearchBuffer())
    end
end

@inline function textbuffer(f)
    buff = take!(CACHES)
    try
        f(buff)
    finally
        put!(CACHES, buff)
    end
end

include("textconfig.jl")
include("normalize.jl")
include("tokenize.jl")
include("voc.jl")
include("bow.jl")
include("sparseconversions.jl")
include("vmodel.jl")
include("emodel.jl")

end

max_stars_repo_name: sadit/TextModel

max_issues_repo_name: sadit/TextModel

max_forks_repo_name: sadit/TextModel


text:
module SwoLeios

export cycle_set!, next_set, output_set

using DataFrames
using CSV

include("selection.jl")
include("update.jl")

function output_set(superset)
    (set, reps) = next_set(superset)
    output_set(set, reps)
end

function output_set(set, reps)
    set_str = string(reps)*' '*string(set[:exercise])
    return set_str
end

function next_set(superset)
    set = copy(superset[1][1][1,:])
    reps = select_reps(set)
    return (set, reps)
end

function cycle_set!(superset)
    set = select_set!(superset)
    reps = select_reps(set)
    println(output_set(set, reps))
    return (set, reps)
end

end # module

max_stars_repo_name: leios/SwoLeios

max_issues_repo_name: leios/SwoLeios

max_forks_repo_name: leios/SwoLeios


text:
include("utilities.jl")

using SEAL
using Printf

"""
    example_serialization()

Show example for how to use serialization in a client-server setup, where only the client knows the
secret key and the servers does not have knowledge of any unencrypted data.
This function is based on the file `native/examples/6_serialization.cpp` of the original
SEAL library and should yield the exact same output, except for differences in compression ratios.

* [SEAL](https://github.com/microsoft/SEAL)
* [native/examples/6_serialization.cpp](https://github.com/microsoft/SEAL/blob/master/native/examples/6_serialization.cpp)

See also: [`example_bfv_basics`](@ref), [`example_ckks_basics`](@ref)
"""
function example_serialization()
  print_example_banner("Example: Serialization")

  parms_stream = UInt8[]
  data_stream1 = UInt8[]
  data_stream2 = UInt8[]
  data_stream3 = UInt8[]
  data_stream4 = UInt8[]
  sk_stream = UInt8[]

  # Use `let` to create new variable scope to mimic curly braces-delimited blocks in C++

  # Server
  let
    parms = EncryptionParameters(SchemeType.ckks)
    poly_modulus_degree = 8192
    set_poly_modulus_degree!(parms, poly_modulus_degree)
    set_coeff_modulus!(parms, coeff_modulus_create(poly_modulus_degree, [50, 20, 50]))

    resize!(parms_stream, save_size(parms))
    out_bytes = save!(parms_stream, parms)
    resize!(parms_stream, out_bytes)

    print_line(@__LINE__)
    println("EncryptionParameters: wrote ", out_bytes, " bytes")

    print_line(@__LINE__)
    println("EncryptionParameters: data size upper bound (compr_mode_type::none): ",
            save_size(ComprModeType.none, parms))
    println("             EncryptionParameters: data size upper bound (compression): ",
            save_size(parms))

    byte_buffer = Vector{UInt8}(undef, save_size(parms))
    save!(byte_buffer, length(byte_buffer), parms)

    parms2 = EncryptionParameters()
    load!(parms2, byte_buffer, length(byte_buffer))

    print_line(@__LINE__)
    println("EncryptionParameters: parms == parms2: ", parms == parms2)
  end

  # Client
  let
    parms = EncryptionParameters()
    load!(parms, parms_stream)

    context = SEALContext(parms)

    keygen = KeyGenerator(context)
    sk = secret_key(keygen)
    pk = PublicKey()
    create_public_key!(pk, keygen)

    resize!(sk_stream, save_size(sk))
    out_bytes = save!(sk_stream, sk)
    resize!(sk_stream, out_bytes)

    rlk = create_relin_keys(keygen)

    rlk_big = RelinKeys()
    create_relin_keys!(rlk_big, keygen)

    resize!(data_stream1, save_size(rlk))
    size_rlk = save!(data_stream1, rlk)
    resize!(data_stream1, size_rlk)
    resize!(data_stream2, save_size(rlk_big))
    size_rlk_big = save!(data_stream2, rlk_big)
    resize!(data_stream2, size_rlk_big)

    print_line(@__LINE__)
    println("Serializable<RelinKeys>: wrote ", size_rlk, " bytes")
    println("             ", "RelinKeys wrote ", size_rlk_big, " bytes")

    initial_scale = 2.0^20
    encoder = CKKSEncoder(context)
    plain1 = Plaintext()
    plain2 = Plaintext()
    encode!(plain1, 2.3, initial_scale, encoder)
    encode!(plain2, 4.5, initial_scale, encoder)

    encryptor = Encryptor(context, pk)

    encrypted1 = Ciphertext()
    encrypt!(encrypted1, plain1, encryptor)
    resize!(data_stream2, save_size(encrypted1))
    size_encrypted1 = save!(data_stream2, encrypted1)
    resize!(data_stream2, size_encrypted1)

    set_secret_key!(encryptor, sk)
    sym_encrypted2 = encrypt_symmetric(plain2, encryptor)
    resize!(data_stream3, save_size(sym_encrypted2))
    size_sym_encrypted2 = save!(data_stream3, sym_encrypted2)
    resize!(data_stream3, size_sym_encrypted2)

    print_line(@__LINE__)
    println("Serializable<Ciphertext> (public-key): wrote ", size_encrypted1, " bytes")
    println("             ",
            "Serializable<Ciphertext> (seeded secret-key): wrote ", size_sym_encrypted2, " bytes")
  end

  # Server
  let
    parms = EncryptionParameters()
    load!(parms, parms_stream)
    context = SEALContext(parms)

    evaluator = Evaluator(context)

    rlk = RelinKeys()
    encrypted1 = Ciphertext()
    encrypted2 = Ciphertext()

    load!(rlk, context, data_stream1)
    load!(encrypted1, context, data_stream2)
    load!(encrypted2, context, data_stream3)

    encrypted_prod = Ciphertext()
    multiply!(encrypted_prod, encrypted1, encrypted2, evaluator)
    relinearize_inplace!(encrypted_prod, rlk, evaluator)
    rescale_to_next_inplace!(encrypted_prod, evaluator)

    resize!(data_stream4, save_size(encrypted_prod))
    size_encrypted_prod = save!(data_stream4, encrypted_prod)
    resize!(data_stream4, size_encrypted_prod)

    print_line(@__LINE__)
    println("Ciphertext (secret-key): wrote ", size_encrypted_prod, " bytes")
  end

  # Client
  let
    parms = EncryptionParameters()
    load!(parms, parms_stream)
    context = SEALContext(parms)

    sk = SecretKey()
    load!(sk, context, sk_stream)
    decryptor = Decryptor(context, sk)
    encoder = CKKSEncoder(context)

    encrypted_result = Ciphertext()
    load!(encrypted_result, context, data_stream4)

    plain_result = Plaintext()
    decrypt!(plain_result, encrypted_result, decryptor)
    slot_count_ = slot_count(encoder)
    result = Vector{Float64}(undef, slot_count_)
    decode!(result, plain_result, encoder)

    print_line(@__LINE__)
    println("Result: ")
    print_vector(result, 3, 7)
  end

  pt = Plaintext("1x^2 + 3")
  stream = Vector{UInt8}(undef, save_size(pt))
  data_size = save!(stream, pt)
  resize!(stream, data_size)

  header = SEALHeader()
  load_header!(header, stream)

  print_line(@__LINE__)
  println("Size written to stream: ", data_size, " bytes")
  println("             ",
          "Size indicated in SEALHeader: ", header.size, " bytes")
  println()

  return
end


text:
using LibCImGui
using Test

@testset "Dock" begin
    @test IMGUI_HAS_DOCK == 1
end