text:
# an implementation of Uniform Manifold Approximation and Projection
# for Dimension Reduction, L. McInnes, J. Healy, J. Melville, 2018.

struct UMAP_{GraphType<:AbstractMatrix, EmbeddingType<:AbstractMatrix, IndexType}
    graph::GraphType
    embedding::EmbeddingType
    index::IndexType
    n_neighbors::Int
    knns::Matrix{Int32}
    dists::Matrix{Float32}
    a::Float32
    b::Float32

    function UMAP_(graph::GraphType, embedding::EmbeddingType, index::IndexType, n_neighbors, knns, dists, a::Real, b::Real) where {GraphType<:AbstractMatrix, EmbeddingType<:AbstractMatrix, IndexType}
        issymmetric(graph) || isapprox(graph, graph') || error("UMAP_ constructor expected graph to be a symmetric matrix")
        new{GraphType, EmbeddingType, IndexType}(graph, embedding, index, n_neighbors, knns, dists, a, b)
    end
end

const SMOOTH_K_TOLERANCE = 1e-5

"""
    umap(X::AbstractMatrix[, n_components=2]; <kwargs>) -> embedding

Embed the data `X` into a `n_components`-dimensional space. `n_neighbors` controls
how many neighbors to consider as locally connected.

See `UMAP_` for a description of keyword arguments.
"""
function umap(args...; kwargs...)
    # this is just a convenience function for now
    return UMAP_(args...; kwargs...).embedding
end

"""
    UMAP_(data_or_index, [, n_components=2]; <kwargs>) -> UMAP_ object

Create a model representing the embedding of data `X` into `n_components`-dimensional space.
The returned model has the following fields:

- `graph`: the graph representing the fuzzy simplicial set of the manifold of `X`.
- `embedding`: the `n-component`-dimensional embedding of the data `X`.
- `data_or_index`: a data matrix or an index to search knns (see SimilaritySearch package).
- `knns`: a matrix of indices of `X` representing each point's nearest neighbors according to `metric`.
          `knns[j, i]` is the index of point i's jth nearest neighbor.
- `dists`: the respective distances of the above neighbors.
           `dists[j, i]` is the distance of point i's jth nearest neighbor.

It uses all available threads for the projection.

# Keyword Arguments
- `n_neighbors::Integer = 15`: the number of neighbors to consider as locally connected. Larger values capture more global structure in the data, while small values capture more local structure.
- `metric::{SemiMetric, Symbol} = Euclidean()`: the metric to calculate distance in the input space. It is also possible to pass `metric = :precomputed` to treat `X` like a precomputed distance matrix.
- `n_epochs::Integer = 300`: the number of training epochs for embedding optimization
- `learning_rate::Real = 1`: the initial learning rate during optimization
- `layout::AbstractLayout = SpectralLayout()`: how to initialize the output embedding
- `min_dist::Real = 0.1`: the minimum spacing of points in the output embedding
- `spread::Real = 1`: the effective scale of embedded points. Determines how clustered embedded points are in combination with `min_dist`.
- `set_operation_ratio::Real = 1`: interpolates between fuzzy set union and fuzzy set intersection when constructing the UMAP graph (global fuzzy simplicial set). The value of this parameter should be between 1.0 and 0.0: 1.0 indicates pure fuzzy union, while 0.0 indicates pure fuzzy intersection.
- `local_connectivity::Integer = 1`: the number of nearest neighbors that should be assumed to be locally connected. The higher this value, the more connected the manifold becomes. This should not be set higher than the intrinsic dimension of the manifold.
- `repulsion_strength::Real = 1`: the weighting of negative samples during the optimization process.
- `neg_sample_rate::Integer = 5`: the number of negative samples to select for each positive sample. Higher values will increase computational cost but result in slightly more accuracy.
- `a = nothing`: this controls the embedding. By default, this is determined automatically by `min_dist` and `spread`.
- `b = nothing`: this controls the embedding. By default, this is determined automatically by `min_dist` and `spread`.

"""
function UMAP_(
    data_or_index,
    n_components::Integer = 2;
    n_neighbors::Integer = 15,
    n_epochs::Integer = 50,
    learning_rate::Real = 1f0,
    learning_rate_decay::Real = 0.9f0,
    layout::AbstractLayout = SpectralLayout(),
    min_dist::Real = 0.1f0,
    spread::Real = 1f0,
    set_operation_ratio::Real = 1f0,
    local_connectivity::Integer = 1,
    repulsion_strength::Float32 = 1f0,
    neg_sample_rate::Integer = 5,
    a = nothing,
    b = nothing
)
    min_dist = convert(Float32, min_dist)
    spread = convert(Float32, spread)
    set_operation_ratio = convert(Float32, set_operation_ratio)
    learning_rate = convert(Float32, learning_rate)
    learning_rate_decay = convert(Float32, learning_rate_decay)
    repulsion_strength = convert(Float32, repulsion_strength)
    index = data_or_index isa AbstractMatrix ? ExhaustiveSearch(; db=data_or_index, dist=SqEuclidean()) : data_or_index
    n = length(index)
    n > n_neighbors > 0 || throw(ArgumentError("the number of examples must be greater than n_neighbors and n_neighbors must be greater than 0"))
    n_components > 1 || throw(ArgumentError("n_components must be greater than 1"))

    @show (n, n_neighbors, n_components)
    # argument checking
    n_epochs > 0 || throw(ArgumentError("n_epochs must be greater than 0"))
    learning_rate > 0 || throw(ArgumentError("learning_rate must be greater than 0"))
    min_dist > 0 || throw(ArgumentError("min_dist must be greater than 0"))
    0 ≤ set_operation_ratio ≤ 1 || throw(ArgumentError("set_operation_ratio must lie in [0, 1]"))
    local_connectivity > 0 || throw(ArgumentError("local_connectivity must be greater than 0"))
    println(stderr, "*** computing allknn graph")
    timeallknn = @elapsed knns, dists = allknn(index, n_neighbors; parallel=Threads.nthreads() > 1)
    println(stderr, "*** computing graph")
    timegraph = @elapsed graph = fuzzy_simplicial_set(knns, dists, n, local_connectivity, set_operation_ratio)
    println(stderr, "*** layout embedding $(typeof(layout))")
    timeinit = @elapsed embedding = initialize_embedding(layout, graph, knns, dists, n_components)
    println(stderr, "*** fit ab / embedding")
    a, b = fit_ab(min_dist, spread, a, b)
    println(stderr, "*** opt embedding")
    timeopt = @elapsed embedding = optimize_embedding(graph, embedding, embedding, n_epochs, learning_rate, repulsion_strength, neg_sample_rate, a, b; learning_rate_decay)
    # TODO: if target variable y is passed, then construct target graph
    #       in the same manner and do a fuzzy simpl set intersection
    println(stderr,
    """
UMAP construction time cost report:
- allknn: $timeallknn
- fuzzy graph: $timegraph
- embedding init: $timeinit
- embedding opt: $timeopt
""")
    UMAP_(graph, embedding, index, n_neighbors, knns, dists, a, b)
end

"""
    UMAP_(umap_::UMAP_, n_components)

Reuses a previously computed model with a different number of components
"""
function UMAP_(U::UMAP_, n_components::Integer;
        n_epochs=50,
        learning_rate::Real = 1f0,
        learning_rate_decay::Real = 0.9f0,
        layout::AbstractLayout = SpectralLayout(),
        repulsion_strength::Float32 = 1f0,
        neg_sample_rate::Integer = 5,
        a = U.a,
        b = U.b
    )

    learning_rate = convert(Float32, learning_rate)
    learning_rate_decay = convert(Float32, learning_rate_decay)
    repulsion_strength = convert(Float32, repulsion_strength)

    graph = U.graph
    embedding = initialize_embedding(layout, graph, U.knns, U.dists, n_components)
    embedding = optimize_embedding(graph, embedding, embedding, n_epochs, learning_rate, repulsion_strength, neg_sample_rate, a, b; learning_rate_decay)
    # TODO: if target variable y is passed, then construct target graph
    #       in the same manner and do a fuzzy simplicial set intersection

    UMAP_(graph, embedding, U.index, U.n_neighbors, U.knns, U.dists, a, b)
end

"""
    transform(model::UMAP_, Q; <kwargs>) -> embedding

Use the given model to embed new points into an existing embedding. `Q` is an AbstractDatabase object (from `SimilaritySearch`)
The returned embedding is the embedding of these points in n-dimensional space, where
n is the dimensionality of `model.embedding`. This embedding is created by finding neighbors of `Q` in `model.embedding`
and optimizing cross entropy according to membership strengths according to these neighbors.

# Keyword Arguments
- `n_neighbors::Integer = 15`: the number of neighbors to consider as locally connected. Larger values capture more global structure in the data, while small values capture more local structure.
- `metric::{SemiMetric, Symbol} = Euclidean()`: the metric to calculate distance in the input space. It is also possible to pass `metric = :precomputed` to treat `X` like a precomputed distance matrix.
- `n_epochs::Integer = 30`: the number of training epochs for embedding optimization
- `learning_rate::Real = 1`: the initial learning rate during optimization
- `learning_rate_decay::Real = 0.8`: A decay factor for the `learning_rate` param (on each epoch)
- `layout::AbstractLayout = SpectralLayout()`: how to initialize the output embedding
- `min_dist::Real = 0.1`: the minimum spacing of points in the output embedding
- `spread::Real = 1`: the effective scale of embedded points. Determines how clustered embedded points are in combination with `min_dist`.
- `set_operation_ratio::Real = 1`: interpolates between fuzzy set union and fuzzy set intersection when constructing the UMAP graph (global fuzzy simplicial set). The value of this parameter should be between 1.0 and 0.0: 1.0 indicates pure fuzzy union, while 0.0 indicates pure fuzzy intersection.
- `local_connectivity::Integer = 1`: the number of nearest neighbors that should be assumed to be locally connected. The higher this value, the more connected the manifold becomes. This should not be set higher than the intrinsic dimension of the manifold.
- `repulsion_strength::Real = 1`: the weighting of negative samples during the optimization process.
- `neg_sample_rate::Integer = 5`: the number of negative samples to select for each positive sample. Higher values will increase computational cost but result in slightly more accuracy.
- `a = nothing`: this controls the embedding. By default, this is determined automatically by `min_dist` and `spread`.
- `b = nothing`: this controls the embedding. By default, this is determined automatically by `min_dist` and `spread`.
"""
function transform(model::UMAP_, Q;
                   n_neighbors::Integer = model.n_neighbors,
                   n_epochs::Integer = 30,
                   learning_rate::Real = 1.0,
                   learning_rate_decay::Real = 0.8,
                   set_operation_ratio::Real = 1.0,
                   local_connectivity::Integer = 1,
                   repulsion_strength::Real = 1.0,
                   neg_sample_rate::Integer = 5,
                   a = model.a,
                   b = model.b
    )

    set_operation_ratio = convert(Float32, set_operation_ratio)
    learning_rate = convert(Float32, learning_rate)
    repulsion_strength = convert(Float32, repulsion_strength)
    learning_rate_decay = convert(Float32, learning_rate_decay)
    Q = convert(AbstractDatabase, Q)

    # argument checking
    length(model.index) > n_neighbors > 0            || throw(ArgumentError("n_neighbors must be greater than 0"))
    learning_rate > 0                                || throw(ArgumentError("learning_rate must be greater than 0"))
    0 ≤ set_operation_ratio ≤ 1                      || throw(ArgumentError("set_operation_ratio must lie in [0, 1]"))
    local_connectivity > 0                           || throw(ArgumentError("local_connectivity must be greater than 0"))
    length(model.index) == size(model.embedding, 2)  || throw(ArgumentError("model.index must have same number of columns as model.embedding"))
    #size(model.data, 1) == size(Q, 1)                || throw(ArgumentError("size(model.data, 1) must equal size(Q, 1)"))

    n_epochs = max(1, n_epochs)
    # main algorithm
    n = length(model.index)
    println("===== inside transform")
    knns, dists = searchbatch(model.index, Q, n_neighbors; parallel=Threads.nthreads() > 1)
    graph = fuzzy_simplicial_set(knns, dists, n, local_connectivity, set_operation_ratio, false)
    E = initialize_embedding(graph, model.embedding)
    println("==== optimizing")
    optimize_embedding(graph, E, model.embedding, n_epochs, learning_rate, repulsion_strength, neg_sample_rate, a, b; learning_rate_decay)
end

"""
    fuzzy_simplicial_set(knns, dists,  n_points, local_connectivity, set_op_ratio, apply_fuzzy_combine=true) -> membership_graph::SparseMatrixCSC,

Construct the local fuzzy simplicial sets of each point represented by its distances
to its `n_neighbors` nearest neighbors, stored in `knns` and `dists`, normalizing the distances
on the manifolds, and converting the metric space to a simplicial set.
`n_points` indicates the total number of points of the original data, while `knns` contains
indices of some subset of those points (ie some subset of 1:`n_points`). If `knns` represents
neighbors of the elements of some set with itself, then `knns` should have `n_points` number of
columns. Otherwise, these two values may be inequivalent.
If `apply_fuzzy_combine` is true, use intersections and unions to combine
fuzzy sets of neighbors (default true).

The returned graph will have size (`n_points`, size(knns, 2)).
"""
function fuzzy_simplicial_set(knns::AbstractMatrix,
                              dists::AbstractMatrix,
                              n_points::Integer,
                              local_connectivity,
                              set_operation_ratio,
                              apply_fuzzy_combine=true)
    # @time σs, ρs = smooth_knn_dists(dists, n_neighbors, local_connectivity)
    # @time rows, cols, vals = compute_membership_strengths(knns, dists, σs, ρs)
    rows, cols, vals = compute_membership_strengths(knns, dists, local_connectivity)
    # transform uses n_points != size(knns, 2)
    fs_set = sparse(rows, cols, vals, n_points, size(knns, 2))

    if apply_fuzzy_combine
        dropzeros!(combine_fuzzy_sets(fs_set, convert(Float32, set_operation_ratio)))
    else
        dropzeros!(fs_set)
    end
end

"""
    smooth_knn_dists(dists, k, local_connectivity; <kwargs>) -> knn_dists, nn_dists

Compute the distances to the nearest neighbors for a continuous value `k`. Returns
the approximated distances to the kth nearest neighbor (`knn_dists`)
and the nearest neighbor (nn_dists) from each point.
"""
function smooth_knn_dists_vector(col::AbstractVector, k::Integer, local_connectivity::Integer; niter::Integer=64, bandwidth::Float32=1f0)
    local_connectivity = max(1, min(k, local_connectivity))
    ρ = _find_first_non_zero(col, local_connectivity) #col[local_connectivity]
    σ = smooth_knn_dist_opt_binsearch(col, ρ, k, bandwidth, niter)
    ρ, σ
end

# check if this is necessary with SimilaritySearch, maybe just adding eps for near duplicates, d(u,v) = 0, is enough.
# also check tests
function _find_first_non_zero(v, sp)
    @inbounds for i in sp:length(v)
        v[i] != 0 && return v[i]
    end

    v[1]
end

# calculate sigma for an individual point
function smooth_knn_dist_kernel(dists, ρ, mid)
    D::Float32 = 0.0
    invmid = -1f0/mid
    @fastmath @inbounds @simd for d in dists
        d = d - ρ
        D += d > 0 ? exp(d * invmid) : 1f0
    end

    D
end

@fastmath function smooth_knn_dist_opt_binsearch(dists::AbstractVector, ρ, k, bandwidth, niter)
    target = bandwidth * log2(k)
    lo::Float32 = 0
    mid::Float32 = 1
    hi::Float32 = Inf32

    for _ in 1:niter
        psum = smooth_knn_dist_kernel(dists, ρ, mid)
        abs(psum - target) < SMOOTH_K_TOLERANCE && break
        if psum > target
            hi = mid
            mid = (lo + hi) * 0.5f0
        else
            lo = mid
            if hi === Inf32
                mid += mid
            else
                mid = (lo + hi) * 0.5f0
            end
        end
    end

    # TODO: set according to min k dist scale
    mid
end

"""
    compute_membership_strengths(knns, dists, local_connectivity) -> rows, cols, vals

Compute the membership strengths for the 1-skeleton of each fuzzy simplicial set.
"""
function compute_membership_strengths(knns::AbstractMatrix, dists::AbstractMatrix, local_connectivity)
    n = length(knns)
    rows = Vector{Int32}(undef, n)
    cols = Vector{Int32}(undef, n)
    vals = Vector{Float32}(undef, n)
    n_neighbors, n = size(knns) # WARN n is now different

    Threads.@threads for i in 1:n
        D = @view dists[:, i]
        ρ, σ = smooth_knn_dists_vector(D, n_neighbors, local_connectivity)
        invσ = -1f0 / σ
        ii = (i-1) * n_neighbors
        @inbounds for k in 1:n_neighbors
            # if i == knns[j, i] # THIS CONDITION NEVER HAPPENS WITH SimilaritySearch
            d = D[k] - ρ
            d = d > 0 ? exp(d * invσ) : 1f0
            iii = ii + k
            cols[iii] = i
            rows[iii] = knns[k, i]
            vals[iii] = d
        end
    end

    rows, cols, vals
end


text:
module TestQueues

using Test
import TreeLas.Pwl.QueueUnion: Queues, Range, reset!
import TreeLas.TreeDP: Tree, TreeDPMem
import GraphIdx.Tree: hierarchy_string, ChildrenIndex, parenthesis_string, root_node
import GraphIdx.Tree: dfs_walk_rev

@testset "Queues: init 2 nodes           " begin
    tree = Tree([1, 1])
    n = length(tree)
    @test n == 2
    mem = TreeDPMem(n)
    cidx = ChildrenIndex(tree.parent)
    @test root_node(cidx) == 1
    @test collect(cidx[1]) == [2]
    @test collect(cidx[2]) == []
    vs = Int[]
    reset!(mem, tree)
    dfs_walk_rev(cidx) do v; push!(vs, v); end
    @test vs == [-1, -2, 2, 1]
    @test mem.queues.pq == Range[4:3, 3:2]
end

@testset "Queues: init 10 nodes          " begin
    """
    Compare to cxx/test/test_tree_dp.cpp: TEST(dptree, init_queues)
    """
    t = Tree(Int[1, 1, 2, 3, 4, 1, 8, 9, 4, 9])
    cidx = ChildrenIndex(t.parent, t.root)
    @test strip(hierarchy_string(cidx)) == strip("""
1
├─2
│ └─3
│   └─4
│     ├─5
│     └─9
│       ├─8
│       │ └─7
│       └─10
└─6""")
    @test parenthesis_string(cidx) == "1(2(3(4(5()9(8(7())10()))))6())"
    n = length(cidx)
    pq = Vector{Range}(undef, n)
    proc_order = reset!(pq, cidx, Int[])
    @test proc_order == [6, 10, 7, 8, 9, 5, 4, 3, 2, 1]
    @test n == 10
    rr = Range[20:19, 19:18, 18:17, 17:16, 16:15, 3:2, 12:11, 13:12, 14:13, 9:8]
    @test pq == rr
    """
      1
      ├─2
      │ └─3
      │   └─4
      │     ├─5
      │     └─9
      │       ├─8
      │       │ └─7
      │       └─10
      └─6
    """
end

end

ext: jl

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
module Config

    export add_config, get_config, set_config,
           get_config_dictionary

    if !(@isdefined GNICONFIG)
        const global GNICONFIG = Dict()
    end

    function add_config(name, value)
        if !haskey(GNICONFIG, name)
            GNICONFIG[name] = value
        end
    end

    function set_config(name, value)
        if haskey(GNICONFIG, name)
            GNICONFIG[name] = value
        else
            println("  WARNING: Unknown parameter $name.")
        end
    end

    function get_config(name)
        if haskey(GNICONFIG, name)
            return GNICONFIG[name]
        else
            println("  WARNING: Unknown parameter $name.")
            return nothing
        end
    end

    function get_config_dictionary()
        return Dict(GNICONFIG)
    end

end


text:
using Test
using LinearAlgebra
using SparseArrays
using BandedMatrices
using SummationByPartsOperators

@testset "Coupled Legendre operators" begin
  for T in (Float32, Float64)
    for degree in 1:5
      ymin = T(0)
      ymax = T(2)
      D = legendre_derivative_operator(ymin, ymax, degree+1)

      for N in 1:3
        xmin = T(-1)
        xmax = T(2)
        mesh = UniformMesh1D(xmin, xmax, N)
        cD_continuous = couple_continuously(D, mesh)
        cD_central    = couple_discontinuously(D, mesh)
        cD_plus       = couple_discontinuously(D, mesh, Val(:plus))
        cD_minus      = couple_discontinuously(D, mesh, Val(:minus))
        for cD in (cD_continuous, cD_central, cD_plus, cD_minus)
          print(cD)
          x = grid(cD)
          println(x)
          @test norm(cD * x.^0) < 100N * eps(float(T))
          for k in 1:degree
            @test cD * x.^k ≈ k .* x.^(k-1)
          end
          @test accuracy_order(cD) == accuracy_order(D)

          u = sinpi.(x)
          v = copy(u)
          SummationByPartsOperators.scale_by_mass_matrix!(u, cD)
          SummationByPartsOperators.scale_by_inverse_mass_matrix!(u, cD)
          @test u ≈ v
          @test integrate(u, cD) ≈ sum(mass_matrix(cD) * u)
          @test integrate(u->u^2, u, cD) ≈ sum(u' * mass_matrix(cD) * u)
          @test cD * u ≈ BandedMatrix(cD) * u atol=eps(float(T)) rtol=sqrt(eps(float(T)))
          @test cD * u ≈ Matrix(cD) * u atol=eps(float(T)) rtol=sqrt(eps(float(T)))
          @test cD * u ≈ sparse(cD) * u atol=eps(float(T)) rtol=sqrt(eps(float(T)))
        end
        M = mass_matrix(cD_central)
        @test M ≈ mass_matrix(cD_plus)
        @test M ≈ mass_matrix(cD_minus)
        @test sum(M) ≈ xmax - xmin
        Dp = Matrix(cD_plus)
        Dm = Matrix(cD_minus)
        diss = M * (Dp - Dm)
        @test diss ≈ diss'
        @test maximum(eigvals(Symmetric(diss))) < 100N * eps(float(T))
        Dc = Matrix(cD_central)
        res = M * Dc + Dc' * M
        res[1, 1] += 1
        res[end, end] -= 1
        @test norm(res) < 100N * eps(float(T))
        res = M * Dp + Dm' * M
        res[1, 1] += 1
        res[end, end] -= 1
        @test norm(res) < 100N * eps(float(T))

        cD1 = couple_discontinuously(cD_central, mesh)
        cD2 = couple_discontinuously(D, UniformMesh1D(xmin, xmax, N^2))
        @test grid(cD1) ≈ grid(cD2)
        @test mass_matrix(cD1) ≈ mass_matrix(cD2)
        @test Matrix(cD1) ≈ Matrix(cD2)

        Mcont = mass_matrix(cD_continuous)
        @test sum(Mcont) ≈ xmax - xmin
        Dcont = Matrix(cD_continuous)
        res = Mcont * Dcont + Dcont' * Mcont
        res[1, 1] += 1
        res[end, end] -= 1
        @test norm(res) < 100N * eps(float(T))
      end

      for N in 1:3
        xmin = T(-1)
        xmax = T(2)
        mesh = UniformPeriodicMesh1D(xmin, xmax, N)
        cD_continuous = couple_continuously(D, mesh)
        cD_central    = couple_discontinuously(D, mesh)
        cD_plus       = couple_discontinuously(D, mesh, Val(:plus))
        cD_minus      = couple_discontinuously(D, mesh, Val(:minus))
        for cD in (cD_continuous, cD_central, cD_plus, cD_minus)
          print(cD)
          x = grid(cD)
          println(x)
          @test norm(cD * x.^0) < 100N * eps(float(T))
          @test accuracy_order(cD) == accuracy_order(D)

          u = sinpi.(x)
          v = copy(u)
          SummationByPartsOperators.scale_by_mass_matrix!(u, cD)
          SummationByPartsOperators.scale_by_inverse_mass_matrix!(u, cD)
          @test u ≈ v
          @test integrate(u, cD) ≈ sum(mass_matrix(cD) * u)
          @test integrate(u->u^2, u, cD) ≈ sum(u' * mass_matrix(cD) * u)
          @test cD * u ≈ Matrix(cD) * u atol=eps(float(T)) rtol=sqrt(eps(float(T)))
          @test cD * u ≈ sparse(cD) * u atol=eps(float(T)) rtol=sqrt(eps(float(T)))
        end
        M = mass_matrix(cD_central)
        @test M ≈ mass_matrix(cD_plus)
        @test M ≈ mass_matrix(cD_minus)
        @test sum(M) ≈ xmax - xmin
        Dp = Matrix(cD_plus)
        Dm = Matrix(cD_minus)
        diss = M * (Dp - Dm)
        @test diss ≈ diss'
        @test maximum(eigvals(Symmetric(diss))) < 100N * eps(float(T))
        Dc = Matrix(cD_central)
        res = M * Dc + Dc' * M
        @test norm(res) < 100N * eps(float(T))
        res = M * Dp + Dm' * M
        @test norm(res) < 100N * eps(float(T))

        Mcont = mass_matrix(cD_continuous)
        @test sum(Mcont) ≈ xmax - xmin
        Dcont = Matrix(cD_continuous)
        res = Mcont * Dcont + Dcont' * Mcont
        @test norm(res) < 100N * eps(float(T))
      end
    end
  end
end

@testset "Coupled FD operators" begin
  for T in (Float32, Float64)
    for acc_order in 2:2:6
      for source in (MattssonNordström2004(),
                     MattssonAlmquistCarpenter2014Extended(),
                     MattssonAlmquistCarpenter2014Optimal())
        ymin = T(0)
        ymax = T(2)
        D = derivative_operator(source, 1, acc_order, ymin, ymax, 31)

        for N in 1:3
          xmin = T(-1)
          xmax = T(2)
          mesh = UniformMesh1D(xmin, xmax, N)
          cD_continuous = couple_continuously(D, mesh)
          cD_central    = couple_discontinuously(D, mesh)
          cD_plus       = couple_discontinuously(D, mesh, Val(:plus))
          cD_minus      = couple_discontinuously(D, mesh, Val(:minus))
          for cD in (cD_continuous, cD_central, cD_plus, cD_minus)
            print(cD)
            x = grid(cD)
            println(x)
            @test norm(cD * x.^0) < 100N * eps(float(T))
            @test accuracy_order(cD) == accuracy_order(D)

            u = sinpi.(x)
            v = copy(u)
            SummationByPartsOperators.scale_by_mass_matrix!(u, cD)
            SummationByPartsOperators.scale_by_inverse_mass_matrix!(u, cD)
            @test u ≈ v
            @test integrate(u, cD) ≈ sum(mass_matrix(cD) * u)
            @test integrate(u->u^2, u, cD) ≈ sum(u' * mass_matrix(cD) * u)
            @test cD * u ≈ BandedMatrix(cD) * u atol=eps(float(T)) rtol=sqrt(eps(float(T)))
            @test cD * u ≈ Matrix(cD) * u atol=eps(float(T)) rtol=sqrt(eps(float(T)))
            @test cD * u ≈ sparse(cD) * u atol=eps(float(T)) rtol=sqrt(eps(float(T)))
          end
          M = mass_matrix(cD_central)
          @test M ≈ mass_matrix(cD_plus)
          @test M ≈ mass_matrix(cD_minus)
          @test sum(M) ≈ xmax - xmin
          Dp = Matrix(cD_plus)
          Dm = Matrix(cD_minus)
          diss = M * (Dp - Dm)
          @test diss ≈ diss'
          @test maximum(eigvals(Symmetric(diss))) < 100N * eps(float(T))
          Dc = Matrix(cD_central)
          res = M * Dc + Dc' * M
          res[1, 1] += 1
          res[end, end] -= 1
          @test norm(res) < 100N * eps(float(T))
          res = M * Dp + Dm' * M
          res[1, 1] += 1
          res[end, end] -= 1
          @test norm(res) < 100N * eps(float(T))

          cD1 = couple_discontinuously(cD_central, mesh)
          cD2 = couple_discontinuously(D, UniformMesh1D(xmin, xmax, N^2))
          @test grid(cD1) ≈ grid(cD2)
          @test mass_matrix(cD1) ≈ mass_matrix(cD2)
          @test Matrix(cD1) ≈ Matrix(cD2)

          Mcont = mass_matrix(cD_continuous)
          @test sum(Mcont) ≈ xmax - xmin
          Dcont = Matrix(cD_continuous)
          res = Mcont * Dcont + Dcont' * Mcont
          res[1, 1] += 1
          res[end, end] -= 1
          @test norm(res) < 100N * eps(float(T))
        end

        for N in 1:3
          xmin = T(-1)
          xmax = T(2)
          mesh = UniformPeriodicMesh1D(xmin, xmax, N)
          cD_continuous = couple_continuously(D, mesh)
          cD_central    = couple_discontinuously(D, mesh)
          cD_plus       = couple_discontinuously(D, mesh, Val(:plus))
          cD_minus      = couple_discontinuously(D, mesh, Val(:minus))
          for cD in (cD_continuous, cD_central, cD_plus, cD_minus)
            print(devnull, cD)
            x = grid(cD)
            @test norm(cD * x.^0) < 100N * eps(float(T))
            @test accuracy_order(cD) == accuracy_order(D)

            u = sinpi.(x)
            v = copy(u)
            SummationByPartsOperators.scale_by_mass_matrix!(u, cD)
            SummationByPartsOperators.scale_by_inverse_mass_matrix!(u, cD)
            @test u ≈ v
            @test integrate(u, cD) ≈ sum(mass_matrix(cD) * u)
            @test integrate(u->u^2, u, cD) ≈ sum(u' * mass_matrix(cD) * u)
            @test cD * u ≈ Matrix(cD) * u atol=eps(float(T)) rtol=sqrt(eps(float(T)))
            @test cD * u ≈ sparse(cD) * u atol=eps(float(T)) rtol=sqrt(eps(float(T)))
          end
          M = mass_matrix(cD_central)
          @test M ≈ mass_matrix(cD_plus)
          @test M ≈ mass_matrix(cD_minus)
          @test sum(M) ≈ xmax - xmin
          Dp = Matrix(cD_plus)
          Dm = Matrix(cD_minus)
          diss = M * (Dp - Dm)
          @test diss ≈ diss'
          @test maximum(eigvals(Symmetric(diss))) < 100N * eps(float(T))
          Dc = Matrix(cD_central)
          res = M * Dc + Dc' * M
          @test norm(res) < 100N * eps(float(T))
          res = M * Dp + Dm' * M
          @test norm(res) < 100N * eps(float(T))

          Mcont = mass_matrix(cD_continuous)
          @test sum(Mcont) ≈ xmax - xmin
          Dcont = Matrix(cD_continuous)
          res = Mcont * Dcont + Dcont' * Mcont
          @test norm(res) < 100N * eps(float(T))
        end
      end
    end
  end
end

@testset "Coupled upwind operators" begin
  for T in (Float32, Float64, Rational{Int128})
    for acc_order in 2:7
      ymin = T(0)
      ymax = T(2)
      Dp = derivative_operator(Mattsson2017(:plus ), 1, acc_order, ymin, ymax, 31)
      Dm = derivative_operator(Mattsson2017(:minus), 1, acc_order, ymin, ymax, 31)

      for N in 1:3
        xmin = T(-1)
        xmax = T(2)
        mesh = UniformMesh1D(xmin, xmax, N)
        cDp_continuous = couple_continuously(Dp, mesh)
        cDm_continuous = couple_continuously(Dm, mesh)
        cDp_central    = couple_discontinuously(Dp, mesh)
        cDm_central    = couple_discontinuously(Dm, mesh)
        cDp_plus       = couple_discontinuously(Dp, mesh, Val(:plus))
        cDm_minus      = couple_discontinuously(Dm, mesh, Val(:minus))
        for cD in (cDp_continuous, cDm_continuous, cDp_central, cDm_central, cDp_plus, cDm_minus)
          print(cD)
          x = grid(cD)
          println(x)
          @test norm(cD * x.^0) < 100N * eps(float(T))
          @test accuracy_order(cD) == accuracy_order(Dp) == accuracy_order(Dm)

          u = sinpi.(x)
          v = copy(u)
          SummationByPartsOperators.scale_by_mass_matrix!(u, cD)
          SummationByPartsOperators.scale_by_inverse_mass_matrix!(u, cD)
          @test u ≈ v
          @test integrate(u, cD) ≈ sum(mass_matrix(cD) * u)
          @test integrate(u->u^2, u, cD) ≈ sum(u' * mass_matrix(cD) * u)
          @test cD * u ≈ BandedMatrix(cD) * u atol=eps(float(T)) rtol=sqrt(eps(float(T)))
          @test cD * u ≈ Matrix(cD) * u atol=eps(float(T)) rtol=sqrt(eps(float(T)))
          @test cD * u ≈ sparse(cD) * u atol=eps(float(T)) rtol=sqrt(eps(float(T)))
        end
        for (cDp,cDm) in ((cDp_continuous,cDm_continuous), (cDp_central,cDm_central), (cDp_plus,cDm_minus))
          cDp_dense = Matrix(cDp)
          cDm_dense = Matrix(cDm)
          M = mass_matrix(cDp)
          @test M ≈ mass_matrix(cDm)
          @test sum(M) ≈ xmax - xmin

          diss = M * (cDp_dense - cDm_dense)
          @test diss ≈ diss'
          @test maximum(eigvals(Symmetric(diss))) < 100N * eps(float(T))
          res = M * cDp_dense + cDm_dense' * M
          res[1, 1] += 1
          res[end, end] -= 1
          @test norm(res) < 100N * eps(float(T))
        end
      end

      for N in 1:3
        xmin = T(-1)
        xmax = T(2)
        mesh = UniformPeriodicMesh1D(xmin, xmax, N)
        cDp_continuous = couple_continuously(Dp, mesh)
        cDm_continuous = couple_continuously(Dm, mesh)
        cDp_central    = couple_discontinuously(Dp, mesh)
        cDm_central    = couple_discontinuously(Dm, mesh)
        cDp_plus       = couple_discontinuously(Dp, mesh, Val(:plus))
        cDm_minus      = couple_discontinuously(Dm, mesh, Val(:minus))
        for cD in (cDp_continuous, cDm_continuous, cDp_central, cDm_central, cDp_plus, cDm_minus)
          print(cD)
          x = grid(cD)
          println(x)
          @test norm(cD * x.^0) < 100N * eps(float(T))
          @test accuracy_order(cD) == accuracy_order(Dp) == accuracy_order(Dm)

          u = sinpi.(x)
          v = copy(u)
          SummationByPartsOperators.scale_by_mass_matrix!(u, cD)
          SummationByPartsOperators.scale_by_inverse_mass_matrix!(u, cD)
          @test u ≈ v
          @test integrate(u, cD) ≈ sum(mass_matrix(cD) * u)
          @test integrate(u->u^2, u, cD) ≈ sum(u' * mass_matrix(cD) * u)
          @test cD * u ≈ Matrix(cD) * u atol=eps(float(T)) rtol=sqrt(eps(float(T)))
          @test cD * u ≈ sparse(cD) * u atol=eps(float(T)) rtol=sqrt(eps(float(T)))
        end
        for (cDp,cDm) in ((cDp_continuous,cDm_continuous), (cDp_central,cDm_central), (cDp_plus,cDm_minus))
          cDp_dense = Matrix(cDp)
          cDm_dense = Matrix(cDm)
          M = mass_matrix(cDp)
          @test M ≈ mass_matrix(cDm)
          @test sum(M) ≈ xmax - xmin

          diss = M * (cDp_dense - cDm_dense)
          @test diss ≈ diss'
          @test maximum(eigvals(Symmetric(diss))) < 100N * eps(float(T))
          res = M * cDp_dense + cDm_dense' * M
          @test norm(res) < 100N * eps(float(T))
        end
      end
    end
  end
end

@testset "Coupled 2nd derivative FD operators" begin
  for T in (Float32, Float64)
    xmin = zero(T)
    xmax = one(T)
    D2op1 = derivative_operator(MattssonNordström2004(), 2, 2, xmin, xmax, 9)
    D2op2 = couple_continuously(D2op1, UniformMesh1D(xmin, xmax, 1))
    @test Matrix(D2op1) ≈ Matrix(D2op2) ≈ BandedMatrix(D2op1) ≈ BandedMatrix(D2op2)

    D2op2 = couple_continuously(D2op1, UniformPeriodicMesh1D(xmin, xmax, 1))

    D2op2 = couple_continuously(derivative_operator(MattssonNordström2004(), 2, 2, xmin, xmax, 5),
                               UniformMesh1D(xmin, xmax, 2))
    @test Matrix(D2op1) ≈ Matrix(D2op2) ≈ BandedMatrix(D2op1) ≈ BandedMatrix(D2op2)
  end
end

@testset "Coupled 2nd derivative Legendre operators" begin
  for T in (Float32, Float64)
    xmin = zero(T)
    xmax = one(T)
    D2op1 = legendre_second_derivative_operator(xmin, xmax, 5)
    D2op2 = couple_continuously(D2op1, UniformMesh1D(xmin, xmax, 1))
    @test Matrix(D2op1) ≈ Matrix(D2op2)

    D2op2 = couple_continuously(D2op1, UniformPeriodicMesh1D(xmin, xmax, 1))
  end
end


text:
using Base.Iterators, Combinatorics
export close!, close, up!, up, reset!, reset, zone_graph_fwd, get_clock_ceilings, get_diagonal_constraints, get_diagonal_guards, normalize_zone, normalize_zone_k, reachable_states, unreachable_states, prune, prune!, find_isomorphism

function close!(c::TAConstraint)
    n = size(c.D,1)
    for k ∈ 1:n
        for i ∈ 1:n
            for j ∈ 1:n
                c.D[i,j] = min(c.D[i,j], c.D[i,k]+c.D[k,j])
            end
        end
    end
    c
end
close(c::TAConstraint) = close!(deepcopy(c))

function up!(c::TAConstraint)
    c.D[2:end,1] .= typemax(eltype(c.D))
    close!(c)
    c
end
up(c::TAConstraint) = up!(deepcopy(c))

function reset!(c::TAConstraint, clocks::Set{Symbol})
    for clock ∈ clocks
        idx = findfirst(==(clock), c.clocks)+1
        c.D[idx,:] .= c.D[1,:]
        c.D[:,idx] .= c.D[:,1]
    end
    close!(c)
    c
end

reset(c::TAConstraint, clocks::Set{Symbol}) = reset!(deepcopy(c), clocks)

function get_diagonal_guards(ta::TA{State,T}) where {State,T}
    Gd = Vector{TADiagonalConstraint{T}}()
    if !isempty(ta.clocks)
        for arc in ta.arcs
            tmp = get_diagonal_constraints(arc.guard)[2:end]
            filter!(tmp) do arc
                !isnothing(arc.clock1) && !isnothing(arc.clock2) && arc.clock1 != arc.clock2
            end

            append!(Gd, tmp)
        end
    end
    unique!(Gd)
    Gd
end

function split(D::TAConstraint{T},Gd) where {T}
    Q = Set{TAConstraint{T}}([D])
    Q′= Set{TAConstraint{T}}()

    for g in Gd
        for D′ in Q
            D1 = D′ ∩ g
            D2 = D′ ∩ !g
            if !isempty(D1) && !isempty(D2)
                # if g cuts the region in two, keep both parts
                push!(Q′,D1)
                push!(Q′,D2)
            else
                # otherwise keep the region
                push!(Q′,deepcopy(D′))
            end
        end
        # swap buffers and empty Q′
        (Q,Q′) = (Q′,Q)
        empty!(Q′)
    end
    Q
end

function normalize_zone_k!(c::TAConstraint{T},k::Vector{T}) where T
    for i in 1:length(c.clocks)+1
        for j in 1:length(c.clocks)+1
            if !isinf(c.D[i,j]) && TABound(false,k[i]) < c.D[i,j]
                c.D[i,j] = typemax(TABound{T})
            elseif !isinf(c.D[i,j]) && c.D[i,j] < TABound(true,-k[j])
                c.D[i,j] = TABound(true,-k[j])
            end
        end
    end
    close!(c)
    c
end

function normalize_zone_k(c::TAConstraint,k)
    _c=deepcopy(c)
    normalize_zone_k!(_c,k)
    _c
end

function normalize_zone(D::TAConstraint{T},Gd,k) where T
    Q = Vector{TAConstraint{T}}()
    for D′ in split(D,Gd)
        G_unsat = Set{TADiagonalConstraint{T}}()
        for g in Gd
            ng = !g
            if isempty(D′ ∩ g)
                push!(G_unsat, ng)
            end
            if isempty(D′ ∩ ng)
                push!(G_unsat, g)
            end
        end

        normalize_zone_k!(D′,k)
        for g in G_unsat
            intersect!(D′,g)
        end
        push!(Q,D′)
    end
    Q
end

function get_clock_ceilings(c::TAConstraint{T}) where T
    max.(
        maximum(x->isinf(x) ? -1 : x.value, c.D, dims=2)[:],
        maximum(x->isinf(x) ? -1 : -x.value, c.D, dims=1)[:]
    )
end

function get_clock_ceilings(ta::TA{State,T}) where {State, T}
    k = zeros(T, length(ta.clocks)+1)
    for arc in ta.arcs
        k .= max.(k, get_clock_ceilings(arc.guard))
    end

    for inv in values(ta.invariants)
        k .= max.(k, get_clock_ceilings(inv))
    end
    k
end

function zone_graph_fwd(ta::TA{State, T}; initial_zone=ta.invariants[ta.initial_state], k_fun=kk->fill(maximum(kk),size(kk)), ignore_implicit_delays=true, max_iter = 1e4) where {State, T}
    lD_initial = (
        ta.initial_state,
        initial_zone
    )
    states = Tuple{State,TAConstraint{T}}[lD_initial]
    transitions = Pair{Tuple{State,TAConstraint{T}},Tuple{TAMessage,Tuple{State,TAConstraint{T}}}}[]

    Gd = get_diagonal_guards(ta)
    k  = k_fun(get_clock_ceilings(ta))

    function add_transition!((l,D),msg,(l′, D′))
        idx = findfirst(==((l′, D′)), states)
        if isnothing(idx)
            # if we just found a new state + clock region, save for later!
            push!(states, (l′, D′))
        end
        push!(transitions, (l,D) => (msg , (l′, D′)))
    end

    last = 0
    while last < length(states) && max_iter > 0
        max_iter -=1
        last += 1
        (l,D) = states[last]

        # Delay transition:
        # (l,D) ⇝ (l,D′), D′ = D↑ ∩ I(l)
        l′ = l
        D′ = up(D) ∩ ta.invariants[l]

        # we don't need loops within the same zone
        if D′!=D
            # split & normalize the zone
            for D′′ in normalize_zone(D′, Gd, k)
                # if this state is consistent, keep it
                if !isempty(D′′)
                    add_transition!((l,D),TAMessage(),(l′,D′′))
                end
            end
        end

        # Action transitions:
        # (l,D) ⇝ (l′, D′), D′ = r(D∧g) ∧ I(l′)  if  l →g,a,r→ l′
        for arc ∈ ta.arcs_by_state[l]
            # follow the arc
            l′ = arc.target

            # The arc constrains the zone
            # if the arc is compatible with the current time zone, keep going
            D′ = D ∩ arc.guard
            if isempty(D′)
                continue
            end

            # the transition resets some clocks and must satisfy the new state's guards
            if ignore_implicit_delays
                intersect!(up!(TimedAutomata.reset!(D′, arc.resets)), ta.invariants[l′])
            else
                intersect!(TimedAutomata.reset!(D′, arc.resets), ta.invariants[l′])
            end

            if isempty(D′)
                continue
            end

            # split & normalize the zone
            for D′′ in normalize_zone(D′, Gd, k)
                if !isempty(D′′)
                    add_transition!((l,D),arc.message,(l′,D′′))
                end
            end

        end
    end

    if max_iter==0
        println("Warning: reached maximum number of iterations!")
    end

    return (states=states, transitions=transitions)
end

#=
function language(ts::TTS{State,T}; ignore_silent=false) where {State,T}
    CombinedState = Tuple{State, TAConstraint{T}}
    sequences = NamedTuple{(:word,:history,:loops),Tuple{Vector{String},Vector{CombinedState}},Vector{Tuple{Int,Int}}}[]

    head = 0
    while head < length(sequences)
        head += 1
        sequence = sequences[head]
        state=last(sequence.history)

        split=0
        for trans in ts.transitions
            if trans[1]==state
                msg,target=trans[2]

                i = findfirst(==(target), sequence.history)

                if isnothing(i) ||
                    # not loop (yet), or the loop is inexpressible as a regex -> new state
                    split += 1
                    new_sequence = split==1 ? sequence : (word=copy(sequence.word), history=copy(sequence.history), loops=copy(sequence.history))
                    push!(new_sequence.history, target)
                    push!(new_sequence.word, ignore_silent && msg.direction == MSG_SILENT ? "" : repr(msg))

                    if split>1
                        push!(sequences, new_sequence)
                    end
                else
                    # loop detected
                    loop = (i,length(sequence.history)
                    if loop
                    push!(history.loops, (i,length(sequence.history)))
                end
            end
        end
    end

    # simplify
    for (word,history,loops) in sequences

    end
end

function language(ts::TTS{State,T}; ignore_silent=false) where {State,T}
    CombinedState = Tuple{State, TAConstraint{T}}

    function recurse(visited, state)
        if state in visited
            return Set{Tuple{String,Union{CombinedState,Nothing}}}([("",state)])
        end

        push!(visited, state)
        sequences = Set{Tuple{String,Union{CombinedState,Nothing}}}()
        for trans in ts.transitions
            if trans[1]==state
                msg,target=trans[2]
                m = ignore_silent && msg.direction == MSG_SILENT ? "" : repr(msg)
                paths = recurse(copy(visited), target)
                for (tail,loop) in paths
                    res=if loop==state
                        "($(m)$(tail))*",nothing
                    else
                        "$(m)$(tail)",loop
                    end
                    push!(sequences,res)
                end
            end
        end

        return sequences
    end

    res = recurse(Set{CombinedState}(), ts.states[1])
    return first.(res)
end
=#
function Base.:|(tas::(TA{S, T} where S)...; names=[], separate_clocks=true) where {T}
    N = length(tas)

    initial_state = |((ta.initial_state for ta in tas)...)
    TT = typeof(initial_state)

    clock_mappers = if separate_clocks
        [(s::Symbol)->Symbol("$(s)_$(i)") for i in 1:N]
    else
        fill(identity, N)
    end

    product_states = product((ta.states for ta in tas)...)
    all_states = collect(|(p...) for p in product_states)[:]
    all_edges = Vector{TAArc{TT,T}}()
    all_invariants = Dict{TT, TAConstraint{T}}()

    # compute clocks
    all_clocks = union((clock_mappers[i].(tas[i].clocks) for i in 1:N)...)

    # compute internal and external alphabet
    alphabet=Dict{Symbol, Symbol}()
    for tasᵢ in tas
        for arc in tasᵢ.arcs
            if haskey(alphabet, arc.message.symbol)
                # if we use the symbol as both in- and output, it must be an internal symbol (no mixed usage!)
                if alphabet[arc.message.symbol]==:input && arc.message.direction==MSG_OUT || alphabet[arc.message.symbol]==:output && arc.message.direction==MSG_IN
                    alphabet[arc.message.symbol] = :internal
                end
            else
                if arc.message.direction==MSG_IN
                    alphabet[arc.message.symbol] = :input
                elseif arc.message.direction == MSG_OUT
                    alphabet[arc.message.symbol] = :output
                end
            end
        end
    end

    # set the clock list for the combined automaton
    extend_clock_list(c::TAConstraint, automaton_id::Int) = set_clock_list!(deepcopy(c),all_clocks;mapping=clock_mappers[automaton_id])
    # compute the combined invariant for a given state from multiple automata
    compute_invariant(state) = intersect((extend_clock_list(tas[i].invariants[state[i]],i) for i in 1:N)...)

    for state in product_states
        combined_state = |(state...)
        invariant = compute_invariant(state)
        all_invariants[|(state...)] = invariant

        @assert !isempty(invariant) "Invariant empty!"

        # step 1: make list of all available edges for this symbol
        # step 2: take ONE sender out of edge list
        # step 3: build every possible subset of the remaining edges (receivers and senders alike, one edge per process each)
            # step 3.a: compute combined guard = union of guards
            # step 3.b: ignore edge set if combined invariant of current state, guard & reset is incompatible with invariant of target state
        # step 4: for every valid edge set, check every superset
            # step 4.a: if the guard of the superset includes the guard of this set, remove this set
            # step 4.b: else cut out the superset's guard region from this guard region <- this requires splitting this zone!!! -> multiple edges

        # step 1: make list of all available edges for this symbol
        options = DefaultDict{Symbol,NamedTuple{(:sender,:receiver,:silent),NTuple{3,Vector{Vector{Int}}}}}(
            ()->(sender=[Vector{Int}() for i in 1:N],receiver=[Vector{Int}() for i in 1:N],silent=[Vector{Int}() for i in 1:N]))
        for (i,(lᵢ,tasᵢ)) in enumerate(zip(state, tas))
            #check all processes' transitions starting in this state
            for (j,arc) in enumerate(tasᵢ.arcs)
                if arc.source==lᵢ
                    if arc.message.direction == MSG_IN
                        push!(options[arc.message.symbol].receiver[i],j)
                    elseif arc.message.direction == MSG_OUT
                        push!(options[arc.message.symbol].sender[i],j)
                    else
                        push!(options[arc.message.symbol].silent[i],j)
                    end
                end
            end
        end

        #see which transitions work
        for (symbol,transition_options) in pairs(options)
            potential_edges = TAArc{TT,T}[]

            edge_sets = if haskey(alphabet,symbol) && alphabet[symbol] in (:internal, :output)
                # step 2: fix exactly ONE process as sender
                tmp = []
                for (i,transition_ids) in enumerate(transition_options.sender)
                    if isempty(transition_ids)
                        continue
                    end
                    push!(tmp,product((j==i ? transition_ids : [nothing; transition_options.receiver[j]] for j=1:N)...))
                end
                flatten(tmp)
            else
                product(([nothing; transition_options.receiver[j]] for j=1:N)...)
            end

            # step 3: build every possible subset of the remaining edges (receivers and senders alike, one edge per process each)
            for edge_set in edge_sets
                # step 3.a: compute combined guard = union of guards
                guard = TAConstraint(TADiagonalConstraint{T}[],all_clocks)
                dir = MSG_SILENT
                invalid = false
                target_state = Any[]
                resets = Set{Symbol}()
                for (k,edge_id) in enumerate(edge_set)
                    if isnothing(edge_id)
                        push!(target_state, state[k])
                    else
                        edge = tas[k].arcs[edge_id]
                        push!(target_state, edge.target)
                        union!(resets, clock_mappers[k].(edge.resets))
                        if dir==MSG_SILENT || edge.message.direction==MSG_OUT
                            dir=edge.message.direction
                        end
                        intersect!(guard, extend_clock_list(edge.guard,k))
                        if isempty(guard)
                            invalid=true
                            break
                        end
                    end
                end
                if invalid
                    continue
                end

                # step 3.b: ignore edge if combined invariant of current state, guard & reset is incompatible with invariant of target state
                c = guard ∩ invariant
                if isempty(c)
                    continue
                end
                intersect!(TimedAutomata.reset!(c,resets),compute_invariant(target_state))
                if isempty(c)
                    continue
                end

                push!(potential_edges, TAArc(combined_state, |(target_state...), guard, TAMessage(dir,symbol), resets))
            end

            # step 4: for every valid edge set, check every superset
            to_remove = zeros(Bool,length(potential_edges))
            for i in 1:length(potential_edges)
                if to_remove[i]
                    continue
                end
                for j in i+1:length(potential_edges)
                    if to_remove[j]
                        continue
                    end
                    edgeᵢ = potential_edges[i]
                    edgeⱼ = potential_edges[j]

                    # step 4.a: if the guard of the superset includes the guard of this set, remove this set
                    if all((edgeᵢ.target .== combined_state) .| (edgeᵢ.target .== edgeⱼ.target)) &&
                        (edgeⱼ.guard ⊆ edgeᵢ.guard)
                        to_remove[i] = true
                    elseif all((edgeⱼ.target .== combined_state) .| (edgeᵢ.target .== edgeⱼ.target)) &&
                        (edgeᵢ.guard ⊆ edgeⱼ.guard)
                        to_remove[j] = true
                    elseif all((edgeᵢ.target .== combined_state) .| (edgeᵢ.target .== edgeⱼ.target)) ||
                        all((edgeⱼ.target .== combined_state) .| (edgeᵢ.target .== edgeⱼ.target))
                        # step 4.b: else cut out the superset's guard region from this guard region <- this requires splitting this zone!!! -> multiple edges
                        throw(ErrorException("Broadcast with non-trivially intersecting guards is not implemented yet"))
                    else
                        # sets don't have a super/subset relation at all
                    end
                end
            end
            deleteat!(potential_edges, to_remove)
            append!(all_edges, potential_edges)
        end
    end
    TA(all_states, initial_state, all_clocks, Set(keys(alphabet)), all_edges, all_invariants)
end

reachable_states(ta::TA) = first.(TTS(ta).states)
unreachable_states(ta::TA) = setdiff(ta.states, reachable_states(ta))

function prune!(ta::TA, prune_states=unreachable_states(ta))
    tts=TTS(ta)

    # prune states
    prune_at = ta.states .∈ Ref(prune_states)
    deleteat!(ta.states, prune_at)
    delete!.(Ref(ta.arcs_by_state), prune_states)
    delete!.(Ref(ta.invariants), prune_states)

    # prune edges
    prune_at = zeros(Bool, length(ta.arcs))
    for (i,arc) in enumerate(ta.arcs)
        if arc.source ∈ prune_states || arc.target ∈ prune_states
            prune_at[i] = true
        end
    end
    deleteat!(ta.arcs, prune_at)

    for (_, arcs) in pairs(ta.arcs_by_state)
        prune_at = zeros(Bool, length(arcs))
        for (i,arc) in enumerate(arcs)
            if arc.source ∈ prune_states || arc.target ∈ prune_states
                prune_at[i] = true
            end
        end
        deleteat!(arcs, prune_at)
    end
    ta
end
prune(ta::TA) = prune!(deepcopy(ta))

"""
Only compares the subgraphs connected to the respective initial states!
"""
function find_isomorphism(ta1::TA, ta2::TA)
    if length(ta1.states)!=length(ta2.states) || length(ta1.clocks) != length(ta2.clocks)
        return nothing
    end

    is_permutation(M::Matrix{Bool}) = all(sum(M,dims=1)' .== sum(M,dims=2) .== 1)

    function check_possible_assignment(M,clks,assignment)
        M′ = copy(M)

        # set assignments
        for (k,l) in assignment
            # check state
            i1 = copy(ta1.invariants[ta1.states[k]])
            set_clock_list!(i1,ta2.clocks;mapping=clks)

            if !M′[k,l] || i1 != ta2.invariants[ta2.states[l]]
                # assignment not possible!
                return nothing
            end
            M′[k,:] .= false
            M′[:,l] .= false
            M′[k,l] = true
        end

        if is_permutation(M′)
            # we did it!
            return M′
        end

        # check successor-states and recurse
        potential_assignments = []
        for (i,j) in assignment
            state1 = ta1.states[i]
            state2 = ta2.states[j]
            edges1 = ta1.arcs_by_state[state1]
            edges2 = ta2.arcs_by_state[state2]
            if length(edges1) != length(edges2)
                # this assignment didn't work out!
                return
            end

            # find all possible ways to identify the edges of the two states
            E = zeros(Bool, length(edges1), length(edges2))
            for (k,edge1) in enumerate(edges1)
                for (l,edge2) in enumerate(edges2)
                    id1 = findfirst(==(edge1.target), ta1.states)
                    id2 = findfirst(==(edge2.target), ta2.states)
                    e1 = deepcopy(edge1)
                    set_clock_list!(e1,ta2.clocks;mapping=clks)
                    if M′[id1,id2] && e1.message == edge2.message && e1.guard == edge2.guard && e1.resets == edge2.resets
                        E[k,l] = true
                    end
                end
            end

            # cannot map some edge
            if any(sum(E,dims=1).==0) || any(sum(E,dims=2).==0)
                # this assignment didn't work out!
                return
            end

            edge_combos = product((edges2[col] for col in eachcol(E))...)
            push!(potential_assignments,[[(findfirst(==(edges1[i].target), ta1.states),findfirst(==(e.target), ta2.states)) for (i,e) in enumerate(ee) ] for ee in (edge_combos)])
        end

        for potential_assignment in (product(potential_assignments...))
            res = check_possible_assignment(M′,clks, flatten(potential_assignment))
            if !isnothing(res)
                return res
            end
        end
        # we didn't find a working assignment
        return
    end

    id1 = findfirst(==(ta1.initial_state), ta1.states)
    id2 = findfirst(==(ta2.initial_state), ta2.states)

    M = ones(Bool, length(ta1.states), length(ta2.states))

    function clks_fact(d)
        (id)->d[id]
    end

    for permutation in permutations(ta2.clocks)
        d = Dict(zip(ta1.clocks,permutation))
        res = check_possible_assignment(M, clks_fact(d), [(id1,id2)])
        if !isnothing(res)
            s=Dict(ta1.states[ind[1]]=>ta2.states[ind[2]] for ind in findall(res))
            return (states=s,clocks=d)
        end
    end
    return
end


text:
# NDTensors.jl
@deprecate use_tblis() NDTensors.using_tblis()
@deprecate enable_tblis!() NDTensors.enable_tblis()
@deprecate disable_tblis!() NDTensors.disable_tblis()

@deprecate addblock!! insertblock!!
@deprecate addblock! insertblock!


text:
mutable struct EnumType <: NamedType
    type::String
    namespace::String
    name::String
    aliases::Union{Vector{String}, Nothing}
    doc::Union{String, Nothing}
    symbols::Vector{String}
    default::Union{String, Nothing}
end

==(x::EnumType, y::EnumType) = x.name == y.name && x.symbols == y.symbols
EnumType(name, symbols, namespace="", aliases=nothing, doc=nothing, default=nothing) =
    EnumType("enum", namespace, name, aliases, doc, symbols, default)
EnumType() = EnumType("enum", "", "", nothing, nothing, String[], nothing)
StructTypes.StructType(::Type{EnumType}) = StructTypes.Mutable()

struct EnumStructType <: StructTypes.StructType end

struct Enum{names}
    value::Int
end

function Base.show(io::IO, x::Enum{names}) where {names}
    print(io, names[x.value + 1], " = ", x.value)
end

StructTypes.StructType(::Type{<:Enum}) = EnumStructType()

juliatype(x::EnumType) = Enum{symtup(x.symbols)}

function sym(x)
    Base.isidentifier(x) || error("invalid avro enum value: `$x`")
    return String(x)
end

schematype(::EnumStructType, ::Type{Enum{names}}, name=Base.string("Enum_", join(names, "_"))) where {names} =
    EnumType(name, [sym(nm) for nm in names])

function writevalue(B::Binary, ::EnumType, x, buf, pos, len, opts)
    return writevalue(B, long, x.value, buf, pos, len, opts)
end

nbytes(::EnumType, x::Enum) = nbytes(x.value)

function readvalue(B::Binary, ::EnumType, ::Type{Enum{names}}, buf, pos, len, opts) where {names}
    x, pos = readvalue(B, long, Int, buf, pos, len, opts)
    return Enum{names}(x), pos
end

skipvalue(B::Binary, ::EnumType, ::Type{Enum{names}}, buf, pos, len, opts) where {names} =
    skipvalue(B, long, Int, buf, pos, len, opts)


text:
module RobotOS

using PyCall

#Empty imported modules for valid precompilation
const _py_sys = PyCall.PyNULL()
const _py_ros_callbacks = PyCall.PyNULL()
const __rospy__ = PyCall.PyNULL()

include("debug.jl")
include("time.jl")
include("gentypes.jl")
include("rospy.jl")
include("pubsub.jl")
include("services.jl")
include("callbacks.jl")

function __init__()
    #Put julia's ARGS into python's so remappings will work
    copy!(_py_sys, pyimport("sys"))
    if length(ARGS) > 0
        _py_sys.argv = ARGS
    end

    #Fill in empty PyObjects
    if ! (dirname(@__FILE__) in _py_sys."path")
        pushfirst!(_py_sys."path", dirname(@__FILE__))
    end
    copy!(_py_ros_callbacks, pyimport("ros_callbacks"))

    try
        copy!(__rospy__, pyimport("rospy"))
    catch ex
        if (isa(ex, PyCall.PyError) && ex.T.__name__ == "ModuleNotFoundError")
            @error """
                Unable to load the 'rospy' python package!
                Has an environment setup script been run?
                """
        else
            rethrow(ex)
        end
    end

    #Compile the callback notify function, see callbacks.jl
    CB_NOTIFY_PTR[] = @cfunction(_callback_notify, Cint, (Ptr{Cvoid},))
end

end


text:
function backwardpass!(prob::Problem,solver::iLQRSolver)
    if solver.opts.square_root
        return _backwardpass_sqrt!(prob,solver)
    else
        return _backwardpass!(prob,solver)
    end
end

function _backwardpass!(prob::Problem,solver::iLQRSolver)
    N = prob.N

    # Objective
    obj = prob.obj

    X = prob.X; U = prob.U; K = solver.K; d = solver.d

    S = solver.S
    Q = solver.Q

    # Terminal cost-to-go
    S[N].xx .= Q[N].xx
    S[N].x .= Q[N].x

    # Initialize expected change in cost-to-go
    ΔV = zeros(2)

    # Backward pass
    k = N-1
    while k >= 1
        fdx, fdu = solver.∇F[k].xx, solver.∇F[k].xu

        Q[k].x .+= fdx'*S[k+1].x
        Q[k].u .+= fdu'*S[k+1].x
        Q[k].xx .+= fdx'*S[k+1].xx*fdx
        Q[k].uu .+= fdu'*S[k+1].xx*fdu
        Q[k].ux .+= fdu'*S[k+1].xx*fdx

        if solver.opts.bp_reg_type == :state
            # Quu_reg = cholesky(Q[k].uu + solver.ρ[1]*fdu'*fdu,check=false)
            Quu_reg = Q[k].uu + solver.ρ[1]*fdu'*fdu
            Qux_reg = Q[k].ux + solver.ρ[1]*fdu'*fdx
        elseif solver.opts.bp_reg_type == :control
            # Quu_reg = cholesky(Q[k].uu + solver.ρ[1]*I,check=false)
            Quu_reg = Q[k].uu + solver.ρ[1]*Diagonal(ones(prob.model.m))
            Qux_reg = Q[k].ux
        end

        # Regularization
        # if Quu_reg.info == -1
        if !isposdef(Hermitian(Quu_reg))
            # increase regularization

            @logmsg InnerIters "Regularizing Quu "
            regularization_update!(solver,:increase)

            # reset backward pass
            k = N-1
            ΔV[1] = 0.
            ΔV[2] = 0.
            continue
        end

        # Compute gains
        K[k] = -1.0*(Quu_reg\Qux_reg)
        d[k] = -1.0*(Quu_reg\Q[k].u)

        # Calculate cost-to-go (using unregularized Quu and Qux)
        S[k].x .= Q[k].x + K[k]'*Q[k].uu*d[k] + K[k]'*Q[k].u + Q[k].ux'*d[k]
        S[k].xx .= Q[k].xx + K[k]'*Q[k].uu*K[k] + K[k]'*Q[k].ux + Q[k].ux'*K[k]
        S[k].xx .= 0.5*(S[k].xx + S[k].xx')

        # calculated change is cost-to-go over entire trajectory
        ΔV[1] += d[k]'*Q[k].u
        ΔV[2] += 0.5*d[k]'*Q[k].uu*d[k]

        k = k - 1;
    end

    # decrease regularization after backward pass
    regularization_update!(solver,:decrease)

    return ΔV
end

function _backwardpass_sqrt!(prob::Problem,solver::iLQRSolver)
    N = prob.N

    # Objective
    obj = prob.obj

    X = prob.X; U = prob.U; K = solver.K; d = solver.d

    S = solver.S
    Q = solver.Q # cost-to-go expansion

    # Terminal cost-to-go expansion
    S[N].xx .= Q[N].xx
    S[N].x .= Q[N].x

    # Initialize expected change in cost-to-go
    ΔV = zeros(2)

    tmp1 = []; tmp2 = []

    # Backward pass
    k = N-1
    while k >= 1
        fdx, fdu = solver.∇F[k].xx, solver.∇F[k].xu

        Q[k].x .+= fdx'*S[k+1].x
        Q[k].u .+= fdu'*S[k+1].x
        tmp_x = S[k+1].xx*fdx
        tmp_u = S[k+1].xx*fdu
        chol_plus!(Q[k].xx,tmp_x)
        chol_plus!(Q[k].uu,tmp_u)
        Q[k].ux .+= tmp_u'*tmp_x

        if solver.opts.bp_reg_type == :state
            Quu_reg = chol_plus(Q[k].uu, sqrt(solver.ρ[1])*fdu)
            Qux_reg = Q[k].ux + solver.ρ[1]*fdu'*fdx
        elseif solver.opts.bp_reg_type == :control
            Quu_reg = chol_plus(Q[k].uu, sqrt(solver.ρ[1])*Diagonal(ones(prob.model.m)))
            Qux_reg = Q[k].ux
        end

        #TODO regularization scheme
        if cond(Quu_reg) > 1e8
            regularization_update!(solver,:increase)

            # reset backward pass
            k = N-1
            ΔV[1] = 0.
            ΔV[2] = 0.
            continue
        end

        # Compute gains
        K[k] = -Quu_reg\(Quu_reg'\Qux_reg)
        d[k] = -Quu_reg\(Quu_reg'\Q[k].u)

        # Calculate cost-to-go (using unregularized Quu and Qux)
        S[k].x .= Q[k].x + (K[k]'*Q[k].uu')*(Q[k].uu*d[k]) + K[k]'*Q[k].u + Q[k].ux'*d[k]

        try
            tmp1 = (Q[k].xx')\Q[k].ux'
        catch SingularException
            tmp1 = pinv(Array(Q[k].xx'))*Q[k].ux'
        end

        tmp2 = chol_minus(Q[k].uu,tmp1)

        S[k].xx .= chol_plus(Q[k].xx + tmp1*K[k],tmp2*K[k])

        # calculated change is cost-to-go over entire trajectory
        ΔV[1] += d[k]'*Q[k].u

        tmp = Q[k].uu*d[k]
        ΔV[2] += 0.5*tmp'*tmp
        k = k - 1;
    end

    # decrease regularization after backward pass
    regularization_update!(solver,:decrease)

    return ΔV
end

#TODO neither of these are particularly fast, preallocating for the inplace version wasn't much faster either
function chol_plus(A::AbstractMatrix{T},B::AbstractMatrix{T}) where T
    n1,m = size(A)
    n2 = size(B,1)
    P = zeros(n1+n2,m)
    P[1:n1,:] = A
    P[n1+1:end,:] = B
    return qr(P).R
end

function chol_plus!(A::AbstractMatrix{T},B::AbstractMatrix{T}) where T
    A .= qr([A;B]).R
end

#TODO
function chol_minus(A,B)
    AmB = Cholesky(copy(A),:U,0)
    for i = 1:size(B,1)
        lowrankdowndate!(AmB,B[i,:])
    end
    return AmB.U
end

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
# Model Configuration
## obs likelihood model configuration
Random.seed!(1)
Nhs = [256,256,256,256] # hidden layer dims
dilations = [1,3,5,7]
Nks = 3 .*ones(Int,length(Nhs))         # conv kernel dims
strides = ones(Int,length(Nhs))
strides[end] = subsample                # output subsampling
dropout = zeros(length(Nhs))
fout = x -> logsoftmax(x,dims=2)        # function in last layer

# training opts
λ1 = 0f-5      # l1 output regularization
lr = 1f-4      # learning rate
Nb = 16        # batch size
epochs_cur = 5
epochs     = 15
opt = Flux.Optimise.Optimiser(WeightDecay(1e-5),ADAM(lr))
curriculum_training = true


text:
@testset "Tajima's D" begin
    n = (10, 77, 72)
    a1 = (2.828968, 4.914514, 4.846921)
    a2 = (1.539768, 1.631862, 1.630948)
    b1 = (0.407407, 0.342105, 0.342723)
    b2 = (0.279012, 0.228184, 0.228612)
    c1 = (0.053922, 0.138626, 0.136406)
    c2 = (0.047227, 0.086985, 0.085989)
    e1 = (0.019061, 0.028208, 0.028143)
    e2 = (0.004949, 0.003374, 0.003423)
    π = (3.888889, 8.438483, 15.339984)
    S = (16, 103, 88)
    D = (-1.446172, -2.021749, -0.525801)

    for i in eachindex(n)
        @test NaturalSelection.td_a1(n[i]) ≈ a1[i] atol=10e-5
        @test NaturalSelection.td_a2(n[i]) ≈ a2[i] atol=10e-5
        @test NaturalSelection.td_b1(n[i]) ≈ b1[i] atol=10e-5
        @test NaturalSelection.td_b2(n[i]) ≈ b2[i] atol=10e-5
        @test NaturalSelection.td_c1(a1[i], b1[i]) ≈ c1[i] atol=10e-5
        @test NaturalSelection.td_c2(n[i], a1[i], a2[i], b2[i]) ≈ c2[i] atol=10e-5
        @test NaturalSelection.td_e1(a1[i], c1[i]) ≈ e1[i] atol=10e-5
        @test NaturalSelection.td_e2(a1[i], a2[i], c2[i]) ≈ e2[i] atol=10e-5
        @test tajimad(π[i], S[i], n[i]) ≈ D[i] atol=10e-5
    end

    testpop = [dna"ATAATAAAAAAATAATAAAAAAATAAAAAAAATAAAAAAAA",
               dna"AAAAAAAATAAATAATAAAAAAATAAAAAAAAAAAAAAAAA",
               dna"AAAATAAAAATATAATAAAAAAATATAAAAAAAAAAAAAAA",
               dna"AAAAAAAAAAAATAATAAAAAAATAAATAAATAAAAAAAAA",
               dna"AAAATAAAAAAAATATAAAAAAATAAAAAAAAAAAAAAAAA",
               dna"AAAATAAAAAAAAAATAAAAAAAAAAAAAAAAAAATAAAAA",
               dna"AAAAAATAAAAATAATAAAAAAATAAAAAAAAAAAAAAAAA",
               dna"AAAAAAAAAAAAAAATAAAAAAATAAAAAAAAAAAAAAATA",
               dna"AAAAAAAAAAAAAAAAAAAAAAATAAAAAAAAAAAAAAAAA",
               dna"AAAAAAAAAAAAAAATAAAAAAATAATAAAAAAAAAAAAAA"]

    @test tajimad(testpop) ≈ -1.446172 atol=10e-5
end


text:
module Periodogram

export periodogram

""" Returns a maximum-likelihood estimate of the amplitude of a
sinusoidal signal in the series of `xs`, sampled at times `ts`, for
each of the frequencies `fs`. """
function periodogram(ts, xs, fs)
    amps = zeros(fs)

    for i in eachindex(amps)
        cs = cos(2*pi*fs[i]*ts)
        ss = sin(2*pi*fs[i]*ts)
        os = ones(ts)

        A = cat(2, cs, ss, os)
        a,b,m = A\xs

        amps[i] = sqrt(a*a + b*b)
    end

    amps
end

end

max_stars_repo_name: farr/RVChallenge

max_issues_repo_name: farr/RVChallenge

max_forks_repo_name: farr/RVChallenge


text:
include("DataStruct.jl")
include("Glouton.jl")
include("Exchange.jl")
include("Amelioration.jl")

function main()
	pbInit::Problem = Problem("Data/test.txt")
	solInit::Solution = glouton(pbInit)
	solInit = amelioration(solInit, pbInit)

	fnames = getfname("Data")

	for fname in fnames
		prob::Problem = Problem("Data/"*fname)

		# println("File : "*fname*" -> ")
		# print("			| Time : ")
		sol = glouton(prob)
		# println("		| z = "*string(sol.z))
		# print( "		| Time : ")
		newSol = amelioration(sol, prob)
		# println("		| z = "*string(newSol.z))
		println("Neeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeew")

	end
end

function getfname(target::String = "Data")
    # target : string := chemin + nom du repertoire ou se trouve les instances

    # recupere tous les fichiers se trouvant dans le repertoire data
    allfiles = readdir(target)

    # vecteur booleen qui marque les noms de fichiers valides
    flag = trues(size(allfiles))

    k=1
    for f in allfiles
        # traite chaque fichier du repertoire
        if f[1] == "."
            # fichier cache => supprimer
            flag[k] = false
        end
        k = k+1
    end

    # extrait les noms valides et retourne le vecteur correspondant
    finstances = allfiles[flag]
    return finstances
end

function evalZ(sol, prob)
	z = 0
	for ind = 1:length(sol.x)
		z+=sol.x[ind] * prob.C[ind]
	end
	return z
end


text:
###
# HTTP ROUTER

function make_router(settings::PlutoDeploySettings, server_session::ServerSession, notebook_sessions::AbstractVector{<:NotebookSession}; static_dir::Union{String,Nothing}=nothing)
    router = HTTP.Router()

    function get_sesh(request::HTTP.Request)
        uri = HTTP.URI(request.target)

        parts = HTTP.URIs.splitpath(uri.path)
        # parts[1] == "staterequest"
        notebook_hash = parts[2] |> HTTP.unescapeuri

        i = findfirst(notebook_sessions) do sesh
            sesh.current_hash == notebook_hash
        end

        if i === nothing
            #=
            ERROR HINT

            This means that the notebook file used by the web client does not precisely match any of the notebook files running in this server.

            If this is an automated setup, then this could happen inotebooketween deployments.

            If this is a manual setup, then running the .jl notebook file might have caused a small change (e.g. the version number or a whitespace change). Copy notebooks to a temporary directory before running them using the bind server. =#
            @info "Request hash not found. See errror hint in my source code." notebook_hash
            nothing
        else
            notebook_sessions[i]
        end
    end

    function get_bonds(request::HTTP.Request)
        request_body = if request.method == "POST"
            IOBuffer(HTTP.payload(request))
        elseif request.method == "GET"
            uri = HTTP.URI(request.target)

            parts = HTTP.URIs.splitpath(uri.path)
            # parts[1] == "staterequest"
            # notebook_hash = parts[2] |> HTTP.unescapeuri

            @assert length(parts) == 3

            base64decode(parts[3] |> HTTP.unescapeuri)
        end
        bonds_raw = Pluto.unpack(request_body)

        Dict{Symbol,Any}(Symbol(k) => v for (k, v) in bonds_raw)
    end

    "Happens whenever you move a slider"
    function serve_staterequest(request::HTTP.Request)
        sesh = get_sesh(request)

        response = if sesh isa RunningNotebookSession
            notebook = sesh.notebook

            bonds = try
                get_bonds(request)
            catch e
                @error "Failed to deserialize bond values" exception=(e, catch_backtrace())
                return HTTP.Response(500, "Failed to deserialize bond values") |> with_cors! |> with_not_cachable!
            end

            @debug "Deserialized bond values" bonds

            let lag = settings.SliderServer.simulated_lag
                lag > 0 && sleep(lag)
            end

            topological_order, new_state = withtoken(sesh.token) do
                try
                    notebook.bonds = bonds

                    names::Vector{Symbol} = Symbol.(keys(bonds))

                    topological_order = Pluto.set_bond_values_reactive(
                        session=server_session,
                        notebook=notebook,
                        bound_sym_names=names,
                        run_async=false,
                    )::Pluto.TopologicalOrder

                    new_state = Pluto.notebook_to_js(notebook)

                    topological_order, new_state
                catch e
                    @error "Failed to set bond values" exception=(e, catch_backtrace())
                    nothing, nothing
                end
            end
            topological_order === nothing && return (HTTP.Response(500, "Failed to set bond values") |> with_cors! |> with_not_cachable!)

            ids_of_cells_that_ran = [c.cell_id for c in topological_order.runnable]

            @debug "Finished running!" length(ids_of_cells_that_ran)

            # We only want to send state updates about...
            function only_relevant(state)
                new = copy(state)
                # ... the cells that just ran and ...
                new["cell_results"] = filter(state["cell_results"]) do (id, cell_state)
                    id ∈ ids_of_cells_that_ran
                end
                # ... nothing about bond values, because we don't want to synchronize among clients.
                new["bonds"] = Dict{String,Dict{String,Any}}()
                new
            end

            patches = Firebasey.diff(only_relevant(sesh.original_state), only_relevant(new_state))
            patches_as_dicts::Array{Dict} = patches

            HTTP.Response(200, Pluto.pack(Dict{String,Any}(
                "patches" => patches_as_dicts,
                "ids_of_cells_that_ran" => ids_of_cells_that_ran,
            ))) |> with_cachable! |> with_cors! |> with_msgpack!
        elseif sesh isa QueuedNotebookSession
            HTTP.Response(503, "Still loading the notebooks... check back later!") |> with_cors! |> with_not_cachable!
        else
            HTTP.Response(404, "Not found!") |> with_cors! |> with_not_cachable!
        end
    end

    function serve_bondconnections(request::HTTP.Request)
        sesh = get_sesh(request)

        response = if sesh isa RunningNotebookSession
            HTTP.Response(200, Pluto.pack(sesh.bond_connections)) |> with_cors! |> with_cachable! |> with_msgpack!
        elseif sesh isa QueuedNotebookSession
            HTTP.Response(503, "Still loading the notebooks... check back later!") |> with_cors! |> with_not_cachable!
        else
            HTTP.Response(404, "Not found!") |> with_cors! |> with_not_cachable!
        end
    end

    HTTP.@register(router, "GET", "/", r -> let
        done = count(x -> !(x isa QueuedNotebookSession), notebook_sessions)
        if static_dir !== nothing
            path = joinpath(static_dir, "index.html")
            if !isfile(path) || done < length(notebook_sessions)
                path = tempname() * ".html"
                write(path, temp_index(notebook_sessions))
            end
            Pluto.asset_response(path)
        else
            if done < length(notebook_sessions)
                HTTP.Response(503, "Still loading the notebooks... check back later! [$(done)/$(length(notebook_sessions)) ready]")
            else
                HTTP.Response(200, "Hi!")
            end
        end |> with_cors! |> with_not_cachable!
    end)

    # !!!! IDEAAAA also have a get endpoint with the same thing but the bond data is base64 encoded in the URL
    # only use it when the amount of data is not too much :o

    HTTP.@register(router, "POST", "/staterequest/*/", serve_staterequest)
    HTTP.@register(router, "GET", "/staterequest/*/*", serve_staterequest)
    HTTP.@register(router, "GET", "/bondconnections/*/", serve_bondconnections)

    if static_dir !== nothing
        function serve_pluto_asset(request::HTTP.Request)
            uri = HTTP.URI(request.target)

            filepath = Pluto.project_relative_path("frontend", relpath(HTTP.unescapeuri(uri.path), "/pluto_asset/"))
            Pluto.asset_response(filepath)
        end
        HTTP.@register(router, "GET", "/pluto_asset/*", serve_pluto_asset)
        function serve_asset(request::HTTP.Request)
            uri = HTTP.URI(request.target)

            filepath = joinpath(static_dir, relpath(HTTP.unescapeuri(uri.path), "/"))
            Pluto.asset_response(filepath)
        end
        HTTP.@register(router, "GET", "/*", serve_asset)
    end

    router
end

###
# HEADERS

function with_msgpack!(response::HTTP.Response)
    push!(response.headers, "Content-Type" => "application/msgpack")
    response
end

function with_cors!(response::HTTP.Response)
    push!(response.headers, "Access-Control-Allow-Origin" => "*")
    response
end

function with_cachable!(response::HTTP.Response)
    second = 1
    minute = 60second
    hour = 60minute
    day = 24hour
    year = 365day

    push!(response.headers, "Cache-Control" => "public, max-age=$(10year), immutable")
    response
end

function with_not_cachable!(response::HTTP.Response)
    push!(response.headers, "Cache-Control" => "no-store, no-cache, max-age=5")
    response
end

function temp_index(notebook_sessions::Vector{NotebookSession})
    default_index(temp_index.(notebook_sessions))
end
function temp_index(s::QueuedNotebookSession)
    without_pluto_file_extension(s.path) => nothing
end
function temp_index(s::Union{FinishedNotebookSession,RunningNotebookSession})
    without_pluto_file_extension(s.path) => without_pluto_file_extension(s.path)*".html"
end


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

mutable struct Drive <: SwaggerModel
    drive_id::Union{ Nothing, String } # drive_id
    path_on_host::Union{ Nothing, String } # path_on_host
    is_root_device::Union{ Nothing, Bool } # is_root_device
    partuuid::Union{ Nothing, String } # partuuid
    is_read_only::Union{ Nothing, Bool } # is_read_only
    rate_limiter::Union{ Nothing, RateLimiter } # rate_limiter

    function Drive(;drive_id=nothing, path_on_host=nothing, is_root_device=nothing, partuuid=nothing, is_read_only=nothing, rate_limiter=nothing)
        o = new()
        set_field!(o, :drive_id, drive_id)
        set_field!(o, :path_on_host, path_on_host)
        set_field!(o, :is_root_device, is_root_device)
        set_field!(o, :partuuid, partuuid)
        set_field!(o, :is_read_only, is_read_only)
        set_field!(o, :rate_limiter, rate_limiter)
        o
    end
end # type Drive

const _name_map_Drive = Dict{String,Symbol}(["drive_id"=>:drive_id, "path_on_host"=>:path_on_host, "is_root_device"=>:is_root_device, "partuuid"=>:partuuid, "is_read_only"=>:is_read_only, "rate_limiter"=>:rate_limiter])
const _field_map_Drive = Dict{Symbol,String}([:drive_id=>"drive_id", :path_on_host=>"path_on_host", :is_root_device=>"is_root_device", :partuuid=>"partuuid", :is_read_only=>"is_read_only", :rate_limiter=>"rate_limiter"])
Swagger.name_map(::Type{ Drive }) = _name_map_Drive
Swagger.field_map(::Type{ Drive }) = _field_map_Drive

function check_required(o::Drive)
    (o.drive_id === nothing) && (return false)
    (o.path_on_host === nothing) && (return false)
    (o.is_root_device === nothing) && (return false)
    (o.is_read_only === nothing) && (return false)
    true
end

function validate_field(o::Drive, name::Symbol, val)
end


text:
function mk_netsta(S::GphysData)
  netsta = Array{String,1}(undef, S.n)
  cha = Array{String,1}(undef, S.n)
  for i in 1:S.n
    id = split_id(S.id[i])
    netsta[i] = id[1]*"."*id[2]
    cha[i] = lowercase(id[4])
  end
  nsid = unique(netsta)
  return netsta, cha, nsid
end

function asdf_sxml(xbuf::IOBuffer, S::GphysData, chans::Array{Int64,1}, sta::HDF5Group)
  seekstart(xbuf)
  mk_xml!(xbuf, S, chans)
  sta["StationXML"] = take!(xbuf)
  return nothing
end

function asdf_mktrace(S::GphysData, xml_buf::IO, chan_numbers::Array{Int64,1}, wav::HDF5Group, ts::Array{Int64,1}, te::Array{Int64,1}, len::Int64, v::Int64, tag::String)
  nc = length(chan_numbers)
  netsta, cha, nsid = mk_netsta(S)
  trace_names = Array{Array{String,1},1}(undef, nc)
  (v>2) && println("traces to create = channels ", chan_numbers)

  # Build trace_names
  d0 = zeros(Int64, nc)
  for (i,j) in enumerate(chan_numbers)
    trace_names[i] = String[]
    id = S.id[j]
    cc = isempty(tag) ? cha[j] : tag

    # divide t0 and t1 by the time length using div to get d0, d1
    d0[i] = len*(div(ts[i], len))
    d1 = len*(div(te[i], len))
    Δ = round(Int64, 1.0e9/S.fs[j])

    # range of times is then d0:d1
    for d in d0[i]:len:d1
      # ...but we must end Δ ns *before* d1 to prevent a one-sample overlap
      s0 = string(u2d(d/1000000000))
      s1 = string(u2d((d+len-Δ)/1000000000))

      # create string like CI.SDD..HHZ__2019-07-07T00:00:00__2019-07-09T00:00:00__hhz_
      chan_str = join([id, s0, s1, cc], "__")
      push!(trace_names[i], chan_str)
    end
    (v > 2) && println("trace_names[", i, "] = ", trace_names[i])
  end

  # Check for trace names in each ns
  chans = falses(nc)
  for j in 1:length(nsid)
    id = nsid[j]
    fill!(chans, false)
    for (j,i) in enumerate(chan_numbers)
      if (netsta[i] == id) && (S.fs[i] != 0.0)
        chans[j] = true
      end
    end
    (maximum(chans) == false) && continue

    # does net.sta exist?
    if has(wav, id)
      sta = wav[id]

    else
      # Create Waveforms/net.sta
      (v > 0) && println("creating Waveforms/", id)
      sta = g_create(wav, id)

      # Write StationXML to sta
      asdf_sxml(xml_buf, S, chan_numbers[chans], sta)
    end

    # Create Waveforms/net.sta/chan_str
    for i = 1:length(chans)
      if chans[i]
        j = chan_numbers[i]
        (v > 2) && println("S.id[", j, "] = ", S.id[j], ", trace_names[", i, "] = ", trace_names[i])
        T = eltype(S.x[j])
        fs = S.fs[j]
        Δ = round(Int64, 1.0e9/fs)
        nx = div(len, Δ)
        for k = 1:length(trace_names[i])
          chan_str = trace_names[i][k]

          # is there a trace with the corresponding string?
          if has(sta, chan_str) == false
            (v > 1) && println("creating trace:  Waveforms/", id, "/", chan_str)
            sta[chan_str] = ones(T, nx).*T(NaN)

            # create a trace of all NaNs
            attrs(sta[chan_str])["sampling_rate"] = fs
            attrs(sta[chan_str])["starttime"] = d0[i] + (k-1)*len
          end
        end
      end
    end
  end

  return nothing
end

function asdf_write_chan(S::GphysData, sta::HDF5Group, i::Int64, tag::String, eid::String, v::Int64)
  fs = S.fs[i]
  tx = S.t[i]
  t = t_win(tx, fs)
  n_seg = size(tx,1)-1
  if n_seg == 0
    @warn(S.id[i] * ": malformed or empty :t field; can't write to file.")
    return
  end

  for k = 1:n_seg
    t0 = t[k,1]
    s0 = string(u2d(div(t0, 1000000)))
    s1 = string(u2d(div(t[k,2], 1000000)))

    si = S.t[i][k,1]
    ei = S.t[i][k+1,1] - (k == n_seg ? 0 : 1)

    # create string like CI.SDD..HHZ__2019-07-07T00:00:00__2019-07-09T00:00:00__hhz_
    chan_str = join([S.id[i], s0, s1, tag], "__")
    if has(sta, chan_str)
      (v > 0) && println("rewriting ", chan_str)
      o_delete(sta, chan_str)
    end
    sta[chan_str] = S.x[i][si:ei]

    # set dictionary attributes
    D = attrs(sta[chan_str])
    D["sampling_rate"] = fs
    D["starttime"] = t0*1000

    if !isempty(eid)
      D["event_id"] = eid
    end
  end

  return nothing
end

#AuxiliaryData/{...}/{...}/{TAG}
function asdf_getaux( hdf_out::String )
  if isfile(hdf_out)
    io = h5open(hdf_out, "r+")
    fmt = read(attrs(io)["file_format"])
    (fmt == "ASDF") || (close(io); error("invalid ASDF file!"))
    if has(io, "AuxiliaryData")
      aux = io["AuxiliaryData"]
    else
      aux = g_create(io, "AuxiliaryData")
    end
  else
    io = h5open(hdf_out, "cw")
    attrs(io)["file_format"] = "ASDF"
    attrs(io)["file_format_version"] = "1.0.2"
    aux = g_create(io, "AuxiliaryData")
  end
  return io, aux
end

"""
    asdf_waux(hdf_out, path, X)

Write `X` to AuxiliaryData/path in `hdf_out`. If an object already exists at
AuxiliaryData/path, it will be deleted and overwritten with `X`.
"""
function asdf_waux(hdf_out::String, path::String, X::Union{HDF5Type,HDF5Array})
  # Correct leading /
  while startswith(path, "/")
    path = path[nextind(path, 1):lastindex(path)]
  end

  # Correct paths that start with AuxiliaryData
  startswith(path, "AuxiliaryData") && (path = String(split(path, "/", limit=2, keepempty=true)[2]))

  io, aux = asdf_getaux(hdf_out)

  # Remove existing object if it exists
  has(aux, path) && o_delete(aux, path)
  aux[path] = X
  close(io)
  return nothing
end


text:
using ClassicalOrthogonalPolynomials, ContinuumArrays, FillArrays, Test
import ClassicalOrthogonalPolynomials: jacobimatrix, oneto, OrthonormalWeighted
import DomainSets: ℝ

@testset "Hermite" begin
    @testset "Basics" begin
        H = Hermite()
        w = HermiteWeight()
        @test axes(H) == (Inclusion(ℝ), oneto(∞))
        x = axes(H,1)
        @test H[0.1,1:4] ≈ hermiteh.(0:3,0.1) ≈ [1,2*0.1,4*0.1^2-2,8*0.1^3-12*0.1]
        @test w[0.1] ≈ exp(-0.1^2)

        X = jacobimatrix(H)
        @test 0.1 * H[0.1,1:10]' ≈ H[0.1,1:11]'*X[1:11,1:10]

        @test (H'*(w .* H))[1,1] ≈ sqrt(π)
        @test (H'*(w .* H))[2,2] ≈ 2sqrt(π)
        @test (H'*(w .* H))[3,3] ≈ 8sqrt(π)

        D = Derivative(x)
        @test (D*H)[0.1,1:4] ≈ [0,2,8*0.1,24*0.1^2-12]
    end

    @testset "OrthonormalWeighted" begin
        H = Hermite()
        Q = OrthonormalWeighted(H)
        @testset "evaluation" begin
            x = 0.1
            @test Q[x,1] ≈ exp(-x^2/2)/π^(1/4)
            @test Q[x,2] ≈ 2x*exp(-x^2/2)/(sqrt(2)π^(1/4))
            # Trap test of L^2 orthonormality
            x = range(-20,20; length=1_000)
            @test Q[x,1:5]'Q[x,1:5] * step(x) ≈ I
        end

        @testset "Differentiation" begin
            x = axes(Q,1)
            D = Derivative(x)
            D¹ = Q \ (D * Q)
            @test D¹[1:10,1:10] ≈ -D¹[1:10,1:10]'
            D² = Q \ (D^2 * Q)
            X = Q \ (x .* Q)
            @test (D² - X^2)[1:10,1:10] ≈ -Diagonal(1:2:19)
        end
    end
end


text:
"""
$(SIGNATURES)

Insert an nGraph `move` node between `producer` at output `index` and all `consumers`.
Return the newly created node.

If optional argument `concurrent` is provided, the generated `move` node will be happen
across `concurrent`.
"""
function insert_move_node!(
        producer::NodeDescriptor,
        index,
        consumers::Vector{NodeDescriptor},
        consumer_inputs,
        concurrent = nothing
    )

    if isnothing(concurrent)
        move_node = nGraph.move(nGraph.Node(producer), index)
    else
        move_node = nGraph.moveasync(nGraph.Node(producer), index, nGraph.Node(concurrent))
    end
    for (consumer, input) in zip(consumers, consumer_inputs)
        nGraph.splice(nGraph.Node(producer), index, nGraph.Node(consumer), input, move_node)
    end

    return NodeDescriptor(move_node)
end

#####
##### General helpers
#####

const NON_PROFILED_NODES = (
    "Parameter",
    "Constant",
    "Move",
    "MoveAsync",
    "Result",
)

"""
$(SIGNATURES)

Return `true` if nodelike `x` should be profiled.
"""
hasprofile(description::String) = !in(description, NON_PROFILED_NODES)
hasprofile(x::nGraph.NodeLike) = hasprofile(nGraph.description(x))

"""
$(SIGNATURES)

Return `true` if nodelike `x` is a `Move` node.
"""
ismove(description::String) = startswith(description, "Move")
ismove(x::nGraph.NodeLike) = ismove(nGraph.description(x))

"""
$(SIGNATURES)

Return `true` if nodelike `x` is a `MoveAsync` node.
"""
ismoveasync(description::String) = startswith(description, "MoveAsync")
ismoveasync(x::nGraph.NodeLike) = ismoveasync(nGraph.description(x))

"""
$(SIGNATURES)

Return `true` if nodelike `x` is a synchronous `Move` node.
"""
ismovesync(description::String) = description == "Move"
ismovesync(x::nGraph.NodeLike) = ismovesync(nGraph.description(x))

"""
$(SIGNATURES)

Return `true` if nodelike `x` is a constant.
"""
isconstant(description::String) = startswith(description, "Constant")
isconstant(x::nGraph.NodeLike) = isconstant(nGraph.description(x))

"""
$(SIGNATURES)

Return `true` if nodelike `x` is a parameter.
"""
isparam(str::String) = startswith(str, "Parameter")
isparam(t) = isparam(nGraph.description(t))

"""
$(SIGNATURES)

Return `true` is nodelike `x` is a result.
"""
isresult(str::String) = startswith(str, "Result")
isresult(t) = isresult(nGraph.description(t))

"""
$(SIGNATURES)

Return all of input tensor descriptors for `fex`.
"""
input_tensors(fex::nGraph.FluxExecutable) = input_tensors(fex.ex.ngraph_function)
function input_tensors(f::nGraph.NFunction)
    params = NodeDescriptor.(nGraph.get_parameters(f))
    return Iterators.flatten(outputs.(params))
end

"""
$(SIGNATURES)

Return all of the output tensor descriptors for `fex`.
"""
output_tensors(fex::nGraph.FluxExecutable) = output_tensors(fex.ex.ngraph_function)
function output_tensors(f::nGraph.NFunction)
    params = NodeDescriptor.(nGraph.get_results(f))
    return Iterators.flatten(outputs.(params))
end

"""
$(SIGNATURES)

Mark that `tensor` should be allocated into PMM.
"""
make_persistent(tensor::TensorDescriptor) = nGraph.make_persistent(tensor)

#####
##### For managing callbacks
#####

"""
Mutation of the ngraph graph happens via callbacks that are called by ngraph during the
function compilation process.

At the time of writing this docstring, there are two callback sites in the modified ngraph:
one right before memory assignment and one right after memory assignment.

Both of these callbacks are passed an `nGraph.NFunction`.

The `CallbackChain` is a way of setting up these callbacks easily to seemlessly allow 0, 1,
or 2 callbacks to be passed to ngraph without trouble.

Functions will be called in the order they are added via [`AutoTM.Utils.callback!`](@ref)

$(METHODLIST)
"""
mutable struct CallbackChain
    # Vector of callback functions
    fns::Vector
    # Number of times this has been invoked.
    num_invocations::Int64
end
CallbackChain() = CallbackChain([], 1)

"""
$(SIGNATURES)

Add `f` to the [`CallbackChain`](@ref) `chain`.
"""
callback!(chain::CallbackChain, f) = push!(chain.fns, f)

function (G::CallbackChain)(args...)
    if G.num_invocations > length(G.fns)
        return nothing
    else
        f = G.fns[G.num_invocations]
        G.num_invocations += 1
        return f(args...)
    end
end

# In case we need to gracefully exit from a GPU compilation callback function
"""
Custom exception for gracefully exiting out of ngraph callbacks.
"""
struct CompilerExit <: Exception end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
try
    run(`lualatex -v`)
catch ex
    @warn("Unable to run `lualatex -v`. Lualatex is required for visualizing LaserTag.")
    if Sys.islinux()
        @info("Attempting to install lualatex with apt-get.")
        # cmds = [`sudo apt-get install texlive-latex-base`, `sudo apt-get install texlive-binaries`]
        cmds = [`sudo apt-get install texlive-latex-extra`]
        for c in cmds
            @info("Running $c")
            run(c)
        end
    end
end

try
    run(`which pdf2svg`)
catch
    @warn("Unable to run `which pdf2svg`. pdf2svg is required for visualizing LaserTag.")
    if Sys.islinux()
        @info("Attempting to install pdf2svg with apt-get")
        cmd = `sudo apt-get install pdf2svg`
        @info("Running $cmd")
        run(cmd)
    end
end


text:
using Documenter
using Inequality

DocMeta.setdocmeta!(Inequality, :DocTestSetup, :(using Inequality); recursive=true)

makedocs(
    sitename = "Inequality.jl",
    format = Documenter.HTML(),
    modules = [Inequality],
    pages = [
		"Introduction" => "index.md",
		"Getting Started" => "gettingstarted.md",
		"Library" => "library.md",
        "Acknowledgements" => "acknowledgements.md"]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
#=deploydocs(
    repo = "<repository url>"
)=#
deploydocs(
    repo = "github.com/JosepER/Inequality.jl.git"
)


text:
sigmoid(x::Real) = inv(one(x) + exp(-x))
binary_step(x::Real) = x < 0 ? 0 : 1

ActivationType = Union{typeof(sigmoid), typeof(binary_step)}


text:
using JuliaZH
using Test


text:
using ForwardDiff
import Rotations: jacobian, ErrorMap
import Rotations: CayleyMap, ExponentialMap, MRPMap, IdentityMap, QuatVecMap

@testset "Quaternion Maps" begin
    ϕ = @SVector rand(3)
    v = 0.1 * @SVector rand(3)
    g = @SVector rand(3)
    p = params(rand(MRP))

    @testset "Forward Maps" begin
        check_forward_jacobian(emap::ErrorMap, ϕ) =
            ForwardDiff.jacobian(x -> params(emap(x)), ϕ) ≈ jacobian(emap, ϕ)

        # Exponential
        @test check_forward_jacobian(ExponentialMap(), ϕ)

        ϕ = 1e-6 * @SVector rand(3)
        @test check_forward_jacobian(ExponentialMap(), ϕ)

        # Vector Part
        @test check_forward_jacobian(QuatVecMap(), v)

        # Gibbs Vectors
        @test check_forward_jacobian(CayleyMap(), g)

        # MRPs
        @test check_forward_jacobian(MRPMap(), p)

        μ0 = 1 / Rotations.scaling(QuatVecMap)
        jac_eye = [@SMatrix zeros(1, 3); μ0 * Diagonal(@SVector ones(3))]
        @test jacobian(ExponentialMap(), p * 1e-10) ≈ jac_eye
        @test jacobian(MRPMap(), p * 1e-10) ≈ jac_eye
        @test jacobian(CayleyMap(), p * 1e-10) ≈ jac_eye
        @test jacobian(QuatVecMap(), p * 1e-10) ≈ jac_eye
    end

    ############################################################################################
    #                                 INVERSE RETRACTION MAPS
    ############################################################################################
    @testset "Inverse Maps" begin

        # Exponential Map
        Random.seed!(1)
        function test_inverse_map(emap, invmap, q, ϕ)
            μ = Rotations.scaling(emap)
            @test inv(emap)(q) ≈ μ * invmap(params(q))
            @test ForwardDiff.jacobian(invmap, params(q)) * μ ≈ jacobian(inv(emap), q)
            @test emap(inv(emap)(q)) ≈ q
            @test inv(emap)(emap(ϕ)) ≈ ϕ
            b = @SVector rand(3)
            @test ForwardDiff.jacobian(
                q -> jacobian(inv(emap), QuatRotation(q, false))'b,
                params(q),
            ) ≈ Rotations.∇jacobian(inv(emap), q, b)
        end

        q = rand(QuatRotation)

        function invmap(q)
            v = @SVector [q[2], q[3], q[4]]
            s = q[1]
            θ = norm(v)
            M = 2 * atan(θ, s) / θ
            return M * v
        end
        test_inverse_map(ExponentialMap(), invmap, q, ϕ)

        qI = ExponentialMap()(v * 1e-5)
        @test ForwardDiff.jacobian(invmap, params(q)) * Rotations.scaling(ExponentialMap) ≈
              jacobian(inv(ExponentialMap()), q)
        @test ForwardDiff.jacobian(invmap, params(qI)) *
              Rotations.scaling(ExponentialMap) ≈ jacobian(inv(ExponentialMap()), qI)

        # Vector Part
        invmap(q) = sign(q[1]) * @SVector [q[2], q[3], q[4]]
        test_inverse_map(QuatVecMap(), invmap, q, v)

        # Cayley
        invmap(q) = 1 / q[1] * @SVector [q[2], q[3], q[4]]
        test_inverse_map(CayleyMap(), invmap, q, g)

        # MRP
        invmap(q) = 1 / (1 + q[1]) * @SVector [q[2], q[3], q[4]]
        test_inverse_map(MRPMap(), invmap, q, p)

        # Test near origin
        μ0 = Rotations.scaling(CayleyMap)
        jacT_eye = [@SMatrix zeros(1, 3); μ0 * Diagonal(@SVector ones(3))]'
        @test isapprox(jacobian(inv(ExponentialMap()), qI), jacT_eye, atol = 1e-5)
        @test isapprox(jacobian(inv(QuatVecMap()), qI), jacT_eye, atol = 1e-5)
        @test isapprox(jacobian(inv(CayleyMap()), qI), jacT_eye, atol = 1e-5)
        @test isapprox(jacobian(inv(MRPMap()), qI), jacT_eye, atol = 1e-5)

    end
end


text:
module MPCHelper

DIR = @__DIR__
@show DIR

using JuMP
using LinearAlgebra
using LazySets
using Polyhedra
import LazySets.Approximations.overapproximate
import LazySets: Interval
# using ControlSystems
using Plots

include("JuMPHelper.jl")
include("approximate_mRPI_set.jl")
include("set_calculations_for_MPC.jl")
include("analytic_solution_lmax.jl")
include("PlotHelper.jl")
include("MatlabHelper.jl")

# copied functionalities
# include("ControlSystem.jl")

# approximate_mRPI_set
export approx_mRPI
# tightend_sets
export compute_tightened_sets
# JuMPHelper
export constrain_variable_by_set!
# PlotHelper
export plot_lin_contour!, plot_vec!
# analytic_solution_lmax
export lmax_mask, ω_opt, hΩ_opt, get_Ellipsoid
# MatlabHelper
export eye

end # module


text:
function comparesorts(tosort)
    a = shuffle(["i", "m", "q"])
    iavg = mavg = qavg = 0.0
    for c in a
        if c == "i"
            iavg = sum(i -> @elapsed(sort(tosort, alg=InsertionSort)), 1:100) / 100.0
        elseif c == "m"
            mavg = sum(i -> @elapsed(sort(tosort, alg=MergeSort)), 1:100) / 100.0
        elseif c == "q"
            qavg = sum(i -> @elapsed(sort(tosort, alg=QuickSort)), 1:100) / 100.0
        end
    end
    iavg, mavg, qavg
end

allones = fill(1, 40000)
sequential = collect(1:40000)
randomized = collect(shuffle(1:40000))

comparesorts(allones)
comparesorts(allones)
iavg, mavg, qavg = comparesorts(allones)
println("Average sort times for 40000 ones:")
println("\tinsertion sort:\t$iavg\n\tmerge sort:\t$mavg\n\tquick sort\t$qavg")

comparesorts(sequential)
comparesorts(sequential)
iavg, mavg, qavg = comparesorts(sequential)
println("Average sort times for 40000 presorted:")
println("\tinsertion sort:\t$iavg\n\tmerge sort:\t$mavg\n\tquick sort\t$qavg")

comparesorts(randomized)
comparesorts(randomized)
iavg, mavg, qavg = comparesorts(randomized)
println("Average sort times for 40000 randomized:")
println("\tinsertion sort:\t$iavg\n\tmerge sort:\t$mavg\n\tquick sort\t$qavg")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module GIDsTest

using Comm, GIDs
using Base.Test

function test()
    t = "hello"
    gid = makegid(t)
    t = nothing
    @test islocal(gid)
    @test !isnull(gid)
    @test hasobj(gid)
    @test getobj(gid) == "hello"
    free(gid)
    sleep(1)
    @test !hasobj(gid)

    gid = nullgid()
    @test isnull(gid)
    @test islocal(gid)

    gid = makegid(t)
    @test gid.proc == 1
    @test islocal(gid)
    @test !isnull(gid)
    global TERMINATE = false
    rexec(mod1(2, Comm.nprocs())) do
        test2(gid)
    end
    while !TERMINATE yield() end
end

function test2(gid::GID)
    @test gid.proc == 1
    @test islocal(gid) == (Comm.myproc()==1)
    @test !isnull(gid)
    free(gid)
    rexec(1) do
        test3(gid)
    end
end

function test3(gid::GID)
    @test gid.proc == 1
    @test islocal(gid)
    @test !isnull(gid)
    sleep(1)
    @test !hasobj(gid)
    global TERMINATE = true
end

function main()
    test()
end

end