text:
using Random

"""
    function information_set(rng::AbstractRNG, j::Int, n_I::Int)

Generates random information sets for chance and decision nodes.
"""
function information_set(rng::AbstractRNG, j::Int, n_I::Int)
    m = min(rand(rng, 0:n_I), j-1)
    return shuffle(rng, 1:(j-1))[1:m]
end

"""
    function information_set(rng::AbstractRNG, leaf_nodes::Vector{Node}, n::Int)

Generates random information sets for value nodes.
"""
function information_set(rng::AbstractRNG, leaf_nodes::Vector{Node}, n::Int)
    @assert n ≥ 1
    non_leaf_nodes = shuffle(rng, setdiff(1:n, leaf_nodes))
    l = length(non_leaf_nodes)
    if isempty(leaf_nodes)
        m = rand(rng, 1:l)
    else
        m = rand(rng, 0:l)
    end
    return [leaf_nodes; non_leaf_nodes[1:m]]
end

"""
    function random_diagram(rng::AbstractRNG, n_C::Int, n_D::Int, n_V::Int, m_C::Int, m_D::Int)

Generate random decision diagram with `n_C` chance nodes, `n_D` decision nodes, and `n_V` value nodes.
Parameter `m_C` and `m_D` are the upper bounds for the size of the information set.

# Examples
```julia
rng = MersenneTwister(3)
random_diagram(rng, 5, 2, 3, 2)
```
"""
function random_diagram(rng::AbstractRNG, n_C::Int, n_D::Int, n_V::Int, m_C::Int, m_D::Int)
    n = n_C + n_D
    n_C ≥ 0 || throw(DomainError("There should be `n_C ≥ 0` chance nodes."))
    n_D ≥ 0 || throw(DomainError("There should be `n_D ≥ 0` decision nodes"))
    n ≥ 1 || throw(DomainError("There should be at least one chance or decision node `n_C+n_D≥1`."))
    n_V ≥ 1 || throw(DomainError("There should be `n_V ≥ 1` value nodes."))
    m_C ≥ 1 || throw(DomainError("Maximum size of information set should be `m_C ≥ 1`."))
    m_D ≥ 1 || throw(DomainError("Maximum size of information set should be `m_D ≥ 1`."))

    # Create node indices
    U = shuffle(rng, 1:n)
    C_j = sort(U[1:n_C])
    D_j = sort(U[(n_C+1):n])
    V_j = collect((n+1):(n+n_V))

    # Create chance and decision nodes
    C = [ChanceNode(j, information_set(rng, j, m_C)) for j in C_j]
    D = [DecisionNode(j, information_set(rng, j, m_D)) for j in D_j]

    # Assign each leaf node to a random value node
    leaf_nodes = setdiff(1:n, (c.I_j for c in C)..., (d.I_j for d in D)...)
    leaf_nodes_j = Dict(j=>Node[] for j in V_j)
    for i in leaf_nodes
        k = rand(rng, V_j)
        push!(leaf_nodes_j[k], i)
    end

    # Create values nodes
    V = [ValueNode(j, information_set(rng, leaf_nodes_j[j], n)) for j in V_j]

    return C, D, V
end

"""
    function States(rng::AbstractRNG, states::Vector{State}, n::Int)

Generate `n` random states from `states`.

# Examples
```julia
rng = MersenneTwister(3)
S = States(rng, [2, 3], 10)
```
"""
function States(rng::AbstractRNG, states::Vector{State}, n::Int)
    States(rand(rng, states, n))
end

"""
    function Probabilities(rng::AbstractRNG, c::ChanceNode, S::States; n_inactive::Int=0)

Generate random probabilities for chance node `c` with `S` states.

# Examples
```julia
rng = MersenneTwister(3)
c = ChanceNode(2, [1])
S = States([2, 2])
Probabilities(rng, c, S)
```
"""
function Probabilities(rng::AbstractRNG, c::ChanceNode, S::States; n_inactive::Int=0)
    states = S[c.I_j]
    state = S[c.j]
    if !(0 ≤ n_inactive ≤ prod([states...; (state - 1)]))
        throw(DomainError("Number of inactive states must be < prod([S[I_j]...;, S[j]-1])"))
    end

    # Create the random probabilities
    data = rand(rng, states..., state)

    # Create inactive chance states
    if n_inactive > 0
        # Count of inactive states per chance stage.
        a = zeros(Int, states...)
        # There can be maximum of (state - 1) inactive chance states per stage.
        b = repeat(vec(CartesianIndices(a)), state - 1)
        # Uniform random sample of n_inactive states.
        r = shuffle(rng, b)[1:n_inactive]
        for s in r
            a[s] += 1
        end
        for s in CartesianIndices(a)
            indices = CartesianIndices(data[s, :])
            i = shuffle(rng, indices)[1:a[s]]
            data[s, i] .= 0.0
        end
    end

    # Normalize the probabilities
    for s in CartesianIndices((states...,))
        data[s, :] /= sum(data[s, :])
    end

    Probabilities(c.j, data)
end

scale(x::Float64, low::Float64, high::Float64) = x * (high - low) + low

"""
    function Consequences(rng::AbstractRNG, v::ValueNode, S::States; low::Float64=-1.0, high::Float64=1.0)

Generate random consequences between `low` and `high` for value node `v` with `S` states.

# Examples
```julia
rng = MersenneTwister(3)
v = ValueNode(3, [1])
S = States([2, 2])
Consequences(rng, v, S; low=-1.0, high=1.0)
```
"""
function Consequences(rng::AbstractRNG, v::ValueNode, S::States; low::Float64=-1.0, high::Float64=1.0)
    if !(high > low)
        throw(DomainError("high should be greater than low"))
    end
    data = rand(rng, S[v.I_j]...)
    data = scale.(data, low, high)
    Consequences(v.j, data)
end

"""
    function LocalDecisionStrategy(rng::AbstractRNG, d::DecisionNode, S::States)

Generate random decision strategy for decision node `d` with `S` states.

# Examples
```julia
rng = MersenneTwister(3)
d = DecisionNode(2, [1])
S = States([2, 2])
LocalDecisionStrategy(rng, d, S)
```
"""
function LocalDecisionStrategy(rng::AbstractRNG, d::DecisionNode, S::States)
    states = S[d.I_j]
    state = S[d.j]
    data = zeros(Int, states..., state)
    for s in CartesianIndices((states...,))
        s_j = rand(rng, 1:state)
        data[s, s_j] = 1
    end
    LocalDecisionStrategy(d.j, data)
end


text:
using FastRefine
using LazySets
import FastRefine: forward_network, forward_affine_map, ishull

nnet = read_nnet("nnet/mnist-20.nnet")

delta = 0.1

solver5 = SpeGuid(delta)
solver6 = HullTree(delta)
solver7 = DimTree(delta)
solver8 = FastTree(delta)

in_hyper = Hyperrectangle(fill(1.0, 400), fill(1.0, 400))
out_hyper = Hyperrectangle(fill(0.0, 10), fill(10.0, 10))
problem = Problem(nnet, in_hyper, out_hyper)

file = open("results/group2.txt", "a")
print(file, "Test Result of Group 2:\n\n")

#solver5

time5 = 0

solve(solver5, problem)
for i = 1:10
    timed_result =@timed solve(solver5, problem)
    print(file, "SpeGuid - test " * string(i) * " - Time: " * string(timed_result.time) * " s")
    print(file, " - Output: " * string(timed_result.value) * "\n")
    global time5 += timed_result.time
end

print(file, "Average time: " * string(time5/10) * " s\n\n")

#solver6

time6 = 0

solve(solver6, problem)
for i = 1:10
    timed_result =@timed solve(solver6, problem)
    print(file, "HullTree - test " * string(i) * " - Time: " * string(timed_result.time) * " s")
    print(file, " - Output: " * string(timed_result.value) * "\n")
    global time6 += timed_result.time
end

print(file, "Average time: " * string(time6/10) * " s\n\n")

#solver7

time7 = 0

solve(solver7, problem)
for i = 1:10
    timed_result =@timed solve(solver7, problem)
    print(file, "DimTree - test " * string(i) * " - Time: " * string(timed_result.time) * " s")
    print(file, " - Output: " * string(timed_result.value) * "\n")
    global time7 += timed_result.time
end

print(file, "Average time: " * string(time7/10) * " s\n\n")

#solver8

time8 = 0

solve(solver8, problem)
for i = 1:10
    timed_result =@timed solve(solver8, problem)
    print(file, "FastTree - test " * string(i) * " - Time: " * string(timed_result.time) * " s")
    print(file, " - Output: " * string(timed_result.value) * "\n")
    global time8 += timed_result.time
end

print(file, "Average time: " * string(time8/10) * " s\n\n")

close(file)

ext: jl


text:
using Defintions

module Results

   include( "results.jl" )
end


text:
using Printf

function output_result(stepnum, Qbase, cellxmax, cellymax, cellzmax, out_file_front, out_ext, out_dir, Rd, yplus)

    stepnum = string(stepnum)
    while length(stepnum) < 6
        stepnum = "0"*stepnum
    end

    fff = out_dir*"/" * out_file_front*stepnum*out_ext
    open(fff,"w") do f
        write(f,"result:rho[kg/m^3], u[m/s], v[m/s], w[m/s], p[Pa], T[K], yplus[-]\n")
        for i in 2:cellxmax-1
            for j in 2:cellymax-1
                for k in 2:cellzmax-1
                    a1 = @sprintf("%8.8f", Qbase[i,j,k,1])
                    a2 = @sprintf("%8.8f", Qbase[i,j,k,2])
                    a3 = @sprintf("%8.8f", Qbase[i,j,k,3])
                    a4 = @sprintf("%8.8f", Qbase[i,j,k,4])
                    a5 = @sprintf("%8.8f", Qbase[i,j,k,5])
                    T = Qbase[i,j,k,5]/(Qbase[i,j,k,1]*Rd)
                    a6 = @sprintf("%8.8f", T)
                    a7 = @sprintf("%8.8f", yplus[i,j,k])
                    write(f, a1*" "*a2*" "*a3*" "*a4*" "*a5*" "*a6*" "*a7*"\n")
                end
            end
        end
    end
    println("\nwrite "*fff)
end

max_stars_repo_name: hide-dog/test

max_issues_repo_name: hide-dog/test

max_forks_repo_name: hide-dog/test


text:
using DesignStructureMatrix
using Test

@testset "Clustering" begin
 DSM =[ 0 1 0 0 1 1 0;
    0 0 0 1 0 0 1;
    0 1 0 1 0 0 1;
    0 1 1 0 1 0 1;
    0 0 0 1 0 1 0;
    1 0 0 0 1 0 0;
    0 1 1 1 0 0 0];

    label = ["A","B","C","D","E","F","G"];

    A,labelA = Clustering(DSM,label)

        @test A[1,1] == 0
        @test A[7,7] == 0

end

@testset "Sequencing" begin

         DSM = [0 0 0 0 0 0 0 0 0 0 0 0 0;
                0 0 0 1 0 0 0 0 0 0 0 0 0;
                1 0 0 1 1 0 0 1 0 1 0 0 1;
                1 0 0 0 0 0 0 0 0 0 0 0 0;
                1 1 1 1 0 0 1 1 0 1 1 0 1;
                0 1 0 0 0 0 0 0 1 0 0 0 0;
                0 0 1 1 1 0 0 1 0 1 0 0 1;
                0 1 0 1 0 0 0 0 0 1 0 0 0;
                0 0 1 1 0 0 0 0 0 1 0 0 1;
                0 1 0 1 0 0 0 1 0 0 0 0 0;
                0 0 1 1 1 0 0 1 0 1 0 0 1;
                0 0 0 0 0 0 0 0 0 0 0 0 1;
                0 1 1 1 1 0 1 1 1 1 1 1 0];

        label = ["1", "2", "3", "4", "5", "6", "7", "8",
                "9", "10", "11", "12", "13"];

        A,labelA = Sequencing(DSM,label)

    @test labelA == ["1","4","2","8","10","3","5","7","9","11","12","13","6"]

end

@testset "Plot" begin

        using Luxor

        DSM = [0 0 0 0 0 0 0 0 0 0 0 0 0;
               0 0 0 1 0 0 0 0 0 0 0 0 0;
               1 0 0 1 1 0 0 1 0 1 0 0 1;
               1 0 0 0 0 0 0 0 0 0 0 0 0;
               1 1 1 1 0 0 1 1 0 1 1 0 1;
               0 1 0 0 0 0 0 0 1 0 0 0 0;
               0 0 1 1 1 0 0 1 0 1 0 0 1;
               0 1 0 1 0 0 0 0 0 1 0 0 0;
               0 0 1 1 0 0 0 0 0 1 0 0 1;
               0 1 0 1 0 0 0 1 0 0 0 0 0;
               0 0 1 1 1 0 0 1 0 1 0 0 1;
               0 0 0 0 0 0 0 0 0 0 0 0 1;
               0 1 1 1 1 0 1 1 1 1 1 1 0];

       label = ["1", "2", "3", "4", "5", "6", "7", "8",
               "9", "10", "11", "12", "13"];

       p = plotDSM(DSM,label)

   @test  typeof(p) == Luxor.Drawing

end


text:
using Pkg
Pkg.activate("scripts")
#Pkg.add(["Plots", "CSV", "DataFrames", "SimilaritySearch", "TextSearch", "Formatting", "Latexify", "UMAP", "Glob"])
using Plots, JLD2, LinearAlgebra, Glob

files = glob("data/SpanishLang/embeddings/umap-embeddings-*.jld2")

for filename in files
    @info filename
    X, C = load(filename, "e2", "e3")

    for i in (1, 2, 3)
        min_, max_ = extrema(C[i, :])
        # @info i => (min_, max_)
        for j in 1:size(C, 2)
            C[i, j] = (C[i, j] - min_) / (max_ - min_)
        end
    end

    C = [RGB(c...) for c in eachcol(C)]

    title = last(split(replace(basename(filename), ".jld2" => ""), '-'))
    display(scatter(X[1, :], X[2, :], c=C, label="", ms=6, ma=0.7, title=title, fmt=:png))
end

max_stars_repo_name: sadit/RegionalSpanish

max_issues_repo_name: sadit/RegionalSpanish

max_forks_repo_name: sadit/RegionalSpanish


text:
export
    Incidency,
    Hypersurface,
    Contact,
    O1_i,
    O1,
    R1,
    Psi,
    Jet,
    Euler_inv

using LightGraphs: SimpleEdge
using Nemo

"""
    Incidency(g, c, w, s, r)

Equivariant class of the cycle parameterizing curves meeting a linear subspace of codimension `r`.
# Arguments
- `g::SimpleGraph`: the graph.
- `c::Vector{UInt8}`: the coloration.
- `w::Vector{Int64}`: the weights.
- `s::Rational{BigInt}`: the scalars.
- `r::Int64`: the codimension of the subvariety. Alternatively, it can be an array of integers, meaning the multiplication of the equivariant class defined by each element of the array.

# Example
The following Gromov-Witten invariants
```math
\\begin{aligned}
\\int_{\\overline{M}_{0,0}(\\mathbb{P}^{3},1)}\\delta_{*}(\\mathrm{ev}^{*}\\mathcal{O}_{\\mathbb{P}^{3}}(1)^{3})^{2} &= 1 \\\\
\\int_{\\overline{M}_{0,0}(\\mathbb{P}^{3},1)}\\delta_{*}(\\mathrm{ev}^{*}\\mathcal{O}_{\\mathbb{P}^{3}}(1)^{2})^{2}\\cdot \\delta_{*}(\\mathrm{ev}^{*}\\mathcal{O}_{\\mathbb{P}^{3}}(1)^{3}) &= 1 \\\\
\\int_{\\overline{M}_{0,0}(\\mathbb{P}^{3},3)}\\delta_{*}(\\mathrm{ev}^{*}\\mathcal{O}_{\\mathbb{P}^{3}}(1)^{2})^{2}\\cdot \\mathrm{c_{top}}(\\delta_{*}(\\mathrm{ev}^{*}\\mathcal{O}_{\\mathbb{P}^{3}}(3))) &= 756 \\\\
\\end{aligned}
```
can be computed as
```julia-repl
julia> P = (g,c,w,s,m) -> Incidency(g,c,w,s,3)^2;
julia> AtiyahBottFormula(3,1,0,P);
Result: 1//1
julia> P = (g,c,w,s,m) -> Incidency(g,c,w,s,[2,2,3]);
julia> AtiyahBottFormula(3,1,0,P);
Result: 1//1
julia> P = (g,c,w,s,m) -> Incidency(g,c,w,s,[2,2])*Hypersurface(g,c,w,s,3);
julia> AtiyahBottFormula(3,3,0,P);
Result: 756//1
```
!!! warning "Attention!"

    The program will stop if `r` is not positive.

"""
function Incidency(g::SimpleGraph, col::Vector{UInt8}, weights::Vector{Int64}, scalars::Vector{fmpq}, r::Int64)::fmpq

    local p1 = fmpq(0); #the final result
    r -= 1

    # col = Dict(vertices(g).=> coloration) #assign colors to vertices
    d = Dict(edges(g).=> weights) #assign weights to edges

    for e in edges(g)
        for t in (0:r)
            p1 += d[e]*(scalars[col[src(e)]]^(t))*(scalars[col[dst(e)]]^(r-t))
        end
    end

    return p1

end

function Incidency(g::SimpleGraph, coloration::Vector{UInt8}, weights::Vector{Int64}, scalars::Vector{fmpq}, r::Vector{Int64})::fmpq

    local p1 = fmpq(1);

    for j in unique(r)
        p1 *= Incidency(g, coloration, weights, scalars, j)^count(x->x==j, r)
    end
    return p1

end

"""
    Hypersurface(g, c, w, s, b)

Equivariant class of the Euler class of the bundle equal to the direct image under the forgetful map of ``\\mathrm{ev}^*\\mathcal{O}_{\\mathbb{P}^n}(b)``. It parameterizes curves contained in a hypersurface of degree `b`.
# Arguments
- `g::SimpleGraph`: the graph.
- `c::Vector{UInt8}`: the coloration.
- `w::Vector{Int64}`: the weights.
- `s::Rational{BigInt}`: the scalars.
- `b::Int64`: the degrees of the hypersurface. Alternatively, it can be an array of integers, meaning the multiplication of the equivariant class defined by each element of the array.

# Example
The following Gromov-Witten invariants of Calabi-Yau threefolds
```math
\\begin{aligned}
\\int_{\\overline{M}_{0,0}(\\mathbb{P}^{4},1)}\\mathrm{c_{top}}(\\delta_{*}(\\mathrm{ev}^{*}\\mathcal{O}_{\\mathbb{P}^{4}}(5))) &= 2875 \\\\
\\int_{\\overline{M}_{0,0}(\\mathbb{P}^{5},2)}\\mathrm{c_{top}}(\\delta_{*}(\\mathrm{ev}^{*}\\mathcal{O}_{\\mathbb{P}^{5}}(3))^{\\oplus 2}) &= \\frac{423549}{8} \\\\
\\int_{\\overline{M}_{0,0}(\\mathbb{P}^{5},3)}\\mathrm{c_{top}}(\\delta_{*}(\\mathrm{ev}^{*}\\mathcal{O}_{\\mathbb{P}^{5}}(4)))\\cdot\\mathrm{c_{top}}(\\delta_{*}(\\mathrm{ev}^{*}\\mathcal{O}_{\\mathbb{P}^{5}}(2))) &= \\frac{422690816}{27} \\\\
\\int_{\\overline{M}_{0,0}(\\mathbb{P}^{7},4)}\\mathrm{c_{top}}(\\delta_{*}(\\mathrm{ev}^{*}\\mathcal{O}_{\\mathbb{P}^{7}}(2)))^4 &= 25705160 \\\\
\\end{aligned}
```
can be computed as
```julia-repl
julia> P = (g,c,w,s,m) -> Hypersurface(g,c,w,s,5);
julia> AtiyahBottFormula(4,1,0,P);
Result: 2875//1
julia> P = (g,c,w,s,m) -> Hypersurface(g,c,w,s,[3,3]);
julia> AtiyahBottFormula(5,2,0,P);
Result: 423549//8
julia> P = (g,c,w,s,m) -> Hypersurface(g,c,w,s,4)*Hypersurface(g,c,w,s,2);
julia> AtiyahBottFormula(5,3,0,P);
Result: 422690816//27
julia> P = (g,c,w,s,m) -> Hypersurface(g,c,w,s,2)^4;
julia> AtiyahBottFormula(7,4,0,P);
Result: 25705160//1
```
!!! warning "Attention!"

    The program will stop if `b` is not positive.

"""
function Hypersurface(g::SimpleGraph, col::Vector{UInt8}, weights::Vector{Int64}, scalars::Vector{fmpq}, b::Int64)::fmpq

    local p1 = fmpq(1)
    local q1 = fmpq(1)

    # col = Dict(vertices(g).=> coloration) #assign colors to vertices
    d = Dict(edges(g).=> weights) #assign weights to edges

    for e in edges(g)
        for alph in 0:(b*d[e])
            a = fmpq( 1 ); c = fmpq( 1 )
            mul!( a, scalars[col[src(e)]], fmpq(alph) )
            mul!( c, scalars[col[dst(e)]], fmpq(b*d[e]-alph))
            a = addeq!( a, c )
            mul!( a, a, inv(fmpq(d[e])))
            p1 = mul!( p1, p1, a )
        end
    end

    for v in vertices(g)
        a = fmpq( 1 )
        mul!( a, scalars[col[v]], fmpq( b ))
        a = a^(1-length(all_neighbors(g, v)))
        mul!( q1, q1, a )
    end

    return p1*q1

end

function Hypersurface(g::SimpleGraph, coloration::Vector{UInt8}, weights::Vector{Int64}, scalars::Vector{fmpq}, b::Vector{Int64})::fmpq

    local p1 = fmpq(1)

    for j in unique(b)
        p1 *= Hypersurface(g, coloration, weights, scalars, j)^count(x->x==j, b)
    end
    return p1
end

"""
    Contact(g, c, w, s)

Equivariant class of the Euler class of the bundle equal to the direct image under the forgetful map of ``\\mathrm{ev}^*\\mathcal{O}_{\\mathbb{P}^n}(2)`` tensor the dualizing sheaf of the forgetful map. It parameterizes contact curves in an odd dimensional projective space.
# Arguments
- `g::SimpleGraph`: the graph.
- `c::Vector{UInt8}`: the coloration.
- `w::Vector{Int64}`: the weights.
- `s::Rational{BigInt}`: the scalars.

# Example
```math
\\begin{aligned}
\\int_{\\overline{M}_{0,2}(\\mathbb{P}^{3},1)}\\mathrm{ev}_{1}^{*}\\mathcal{O}_{\\mathbb{P}^{3}}(1)^{2}\\cdot\\mathrm{ev}_{2}^{*}\\mathcal{O}_{\\mathbb{P}^{3}}(1)^{3}\\cdot\\mathrm{c_{top}}(\\delta_{*}(\\omega_{\\delta}\\otimes\\mathrm{ev}^{*}\\mathcal{O}_{\\mathbb{P}^{3}}(2))) &= 1 \\\\
\\end{aligned}
```
can be computed as
```julia-repl
julia> P = (g,c,w,s,m) -> O1_i(g,c,w,s,m,1)^2*O1_i(g,c,w,s,m,2)^3*Contact(g,c,w,s);
julia> AtiyahBottFormula(3,1,2,P);
Result: 1//1
```
"""
function Contact(g::SimpleGraph, col::Vector{UInt8}, weights::Vector{Int64}, scalars::Vector{fmpq})::fmpq

    local p1 = fmpq(1)
    local q1 = fmpq(1)

    # col = Dict(vertices(g).=> coloration) #assign colors to vertices
    d = Dict(edges(g).=> weights) #assign weights to edges

    for e in edges(g)
        for alph in (1:2*d[e]-1)
            p1 *= (alph*scalars[col[src(e)]]+(2*d[e]-alph)*scalars[col[dst(e)]])//d[e]
        end
    end

    for v in vertices(g)
        q1 *= (2*scalars[col[v]])^(length(all_neighbors(g, v))-1)
    end

    return p1*q1

end

"""
    O1_i(g, c, w, s, m, i)

Equivariant class of the pull-back of ``\\mathcal{O}_{\\mathbb{P}^n}(1)`` with respect to the i-th evaluation map.
# Arguments
- `g::SimpleGraph`: the graph.
- `c::Vector{UInt8}`: the coloration.
- `w::Vector{Int64}`: the weights.
- `s::Rational{BigInt}`: the scalars.
- `m::Marks`: the marks.
- `i::Int64`: the evaluation map.

# Example
The following Gromov-Witten invariants
```math
\\begin{aligned}
\\int_{\\overline{M}_{0,2}(\\mathbb{P}^{2},1)}\\mathrm{ev}_{1}^{*}\\mathcal{O}_{\\mathbb{P}^{2}}(1)\\cdot\\mathrm{ev}_{2}^{*}\\mathcal{O}_{\\mathbb{P}^{2}}(1) &= 1 \\\\
\\int_{\\overline{M}_{0,1}(\\mathbb{P}^{3},1)}\\mathrm{ev}_{1}^{*}\\mathcal{O}_{\\mathbb{P}^{2}}(1)^{2}\\cdot\\mathrm{c_{top}}(\\delta_{*}(\\mathrm{ev}^{*}\\mathcal{O}_{\\mathbb{P}^{3}}(2))) &= 4
\\end{aligned}
```
can be computed as
```julia-repl
julia> P = (g,c,w,s,m) -> O1_i(g,c,w,s,m,1)*O1_i(g,c,w,s,m,2);
julia> AtiyahBottFormula(2,1,2,P);
Result: 1//1
julia> P = (g,c,w,s,m) -> O1_i(g,c,w,s,m,1)^2*Hypersurface(g,c,w,s,2);
julia> AtiyahBottFormula(3,1,1,P);
Result: 4//1
```
!!! warning "Attention!"

    The program will stop if `i` is not between 1 and the number of marks.

"""
function O1_i(g::SimpleGraph, col::Vector{UInt8}, weights::Vector{Int64}, scalars::Vector{fmpq}, mark::Marks, i::Int64)::fmpq

    # col = Dict(vertices(g).=> coloration) #assing colors to vertices

    return scalars[col[mark.get_vertex[i]]]
end

"""
    O1(g, c, w, s, m)

Equivariant class of the pull-back of ``\\mathcal{O}_{\\mathbb{P}^n}(1)`` with respect to the product of all evaluation maps.
# Arguments
- `g::SimpleGraph`: the graph.
- `c::Vector{UInt8}`: the coloration.
- `w::Vector{Int64}`: the weights.
- `s::Rational{BigInt}`: the scalars.
- `m::Marks`: the marks.

This function is equivalent to the product of the function `O1_i(g,c,w,s,m,i)` where `i` runs from 1 to the number of marks.

# Example
The following Gromov-Witten invariants
```math
\\begin{aligned}
\\int_{\\overline{M}_{0,8}(\\mathbb{P}^{2},3)}\\prod_{i=1}^{8}\\mathrm{ev}_{i}^{*}\\mathcal{O}_{\\mathbb{P}^{2}}(1) &= 12 \\\\
\\int_{\\overline{M}_{0,1}(\\mathbb{P}^{3},2)}\\mathrm{ev}_{1}^{*}\\mathcal{O}_{\\mathbb{P}^{3}}(1)^2\\cdot\\mathrm{c_{top}}(\\delta_{*}(\\mathrm{ev}^{*}\\mathcal{O}_{\\mathbb{P}^{3}}(3))) &= 81 \\\\
\\end{aligned}
```
can be computed as
```julia-repl
julia> P = (g,c,w,s,m) -> O1(g,c,w,s,m);
julia> AtiyahBottFormula(2,3,8,P);
Result: 12//1
julia> P = (g,c,w,s,m) -> O1(g,c,w,s,m)^2*Hypersurface(g,c,w,s,3);
julia> AtiyahBottFormula(3,2,1,P);
Result: 81//1
```

In order to remove `O1_i(g,c,w,s,m,i)` for some `i`, it is enough to divide by that function.

# Example
```julia-repl
julia> P = (g,c,w,s,m) -> O1(g,c,w,s,m)//O1_i(g,c,w,s,m,1);
```
Here `P` is the product of all `O1_i(g,c,w,s,m,i)` where `i` runs from 2 to `m`.
"""
function O1(g::SimpleGraph, col::Vector{UInt8}, weights::Vector{Int64}, scalars::Vector{fmpq}, mark::Marks)::fmpq

    local p1 = fmpq(1)
    # col = Dict(vertices(g).=> coloration)

    for t in 1:mark.m
        p1 *= scalars[col[mark.get_vertex[t]]]
    end

    return p1
end

"""
    R1(g, c, w, s, k)

The equivariant class of the first derived functor of the pull-back of ``\\mathcal{O}_{\\mathbb{P}^n}(-k)``.
# Arguments
- `g::SimpleGraph`: the graph.
- `c::Vector{UInt8}`: the coloration.
- `w::Vector{Int64}`: the weights.
- `s::Rational{BigInt}`: the scalars.
- `k::Int64`: a positive integer.

# Example
```math
\\begin{aligned}
\\int_{\\overline{M}_{0,0}(\\mathbb{P}^{1},d)}\\mathrm{c_{top}}(R^{1}\\delta_{*}(\\mathrm{ev}^{*}\\mathcal{O}_{\\mathbb{P}^{3}}(-1)))^2 &= \\frac{1}{d^3} \\\\
\\end{aligned}
```
can be computed as
```julia-repl
julia> d = 1; #for other values of d, change this line
julia> P = (g,c,w,s,m) -> R1(g,c,w,s,1)^2;
julia> AtiyahBottFormula(1,d,0,P);
Result: 1//1
```
!!! warning "Attention!"

    The program will stop if `k` is not positive.

"""
function R1(g::SimpleGraph, col::Vector{UInt8}, weights::Vector{Int64}, scalars::Vector{fmpq}, k::Int64)::fmpq

    local p1 = fmpq(1)
    local q1 = fmpq(1)

    # col = Dict(vertices(g).=> coloration) #assign colors to vertices
    d = Dict(edges(g).=> weights) #assign weights to edges

    for e in edges(g)
        for alph in 1:(k*d[e]-1)
            p1 *= (-1)*(alph*scalars[col[src(e)]]+(k*d[e]-alph)*scalars[col[dst(e)]])//d[e]
        end
    end

    for v in vertices(g)
        q1 *= (-k*scalars[col[v]])^(length(all_neighbors(g, v))-1)
    end

    return p1*q1

end

"""
    Psi(g, c, w, s, m, a)

Equivariant class of the cycle of ``\\psi``-classes.
# Arguments
- `g::SimpleGraph`: the graph.
- `c::Vector{UInt8}`: the coloration.
- `w::Vector{Int64}`: the weights.
- `s::Rational{BigInt}`: the scalars.
- `m::Marks`: the marks.
- `a::Vector{Int64}`: the vector of the exponents of the ``\\psi`` classes. It is ordered, meaning that the first element is the exponent of ``\\psi_1``, the second is the exponent of ``\\psi_2``, and so on.

!!! note

    The size of `a` must be at most `m`. If it is smaller, missing exponents will be considered as zeros.
    If `a` is a number, it will be considered as the exponent of ``\\psi_1``.

!!! warning "Attention!"

    The program will stop if we have one of the following conditions:

    * the size of `a` is bigger than `m`,
    * `a` contains a negative number.

# Example
The following Gromov-Witten invariants
```math
\\begin{aligned}
\\int_{\\overline{M}_{0,2}(\\mathbb{P}^{6},2)}\\mathrm{ev}_{1}^{*}\\mathcal{O}_{\\mathbb{P}^{6}}(1)^{5}\\cdot\\mathrm{ev}_{2}^{*}\\mathcal{O}_{\\mathbb{P}^{6}}(1)^{2}\\cdot\\mathrm{c_{top}}(\\delta_{*}(\\mathrm{ev}^{*}\\mathcal{O}_{\\mathbb{P}^{6}}(5)))\\cdot\\psi_{1}\\psi_{2}^{0} &= 495000 \\\\
\\int_{\\overline{M}_{0,2}(\\mathbb{P}^{10},2)}\\mathrm{ev}_{1}^{*}\\mathcal{O}_{\\mathbb{P}^{10}}(1)^{8}\\cdot\\mathrm{ev}_{2}^{*}\\mathcal{O}_{\\mathbb{P}^{10}}(1)^{6}\\cdot\\mathrm{c_{top}}(\\delta_{*}(\\mathrm{ev}^{*}\\mathcal{O}_{\\mathbb{P}^{10}}(7)))\\cdot\\psi_{1}^{2} &= 71804533752 \\\\
\\int_{\\overline{M}_{0,1}(\\mathbb{P}^{2},2)}\\mathrm{ev}_{1}^{*}\\mathcal{O}_{\\mathbb{P}^{2}}(1)^{2}\\cdot\\psi_{1}^{4} &= \\frac{1}{8} \\\\
\\int_{\\overline{M}_{0,1}(\\mathbb{P}^{2},2)}\\delta_{*}(\\mathrm{ev}^{*}\\mathcal{O}_{\\mathbb{P}^{2}}(1)^{2})^{4}\\cdot\\mathrm{ev}_{1}^{*}\\mathcal{O}_{\\mathbb{P}^{2}}(1)\\cdot(\\mathrm{ev}_{1}^{*}\\mathcal{O}_{\\mathbb{P}^{2}}(1)+\\psi_{1}) &= 2 \\\\
\\int_{\\overline{M}_{0,1}(\\mathbb{P}^{2},2)}\\mathrm{ev}_{1}^{*}\\mathcal{O}_{\\mathbb{P}^{2}}(1)^{2}\\cdot(\\psi_{1}^{3}\\cdot\\mathrm{ev}_{1}^{*}\\mathcal{O}_{\\mathbb{P}^{2}}(1)+\\psi_{1}^{2}\\cdot\\mathrm{ev}_{1}^{*}\\mathcal{O}_{\\mathbb{P}^{2}}(1)^{2}) &= \\frac{1}{8} \\\\
\\int_{\\overline{M}_{0,1}(\\mathbb{P}^{3},2)}\\mathrm{ev}_{1}^{*}\\mathcal{O}_{\\mathbb{P}^{2}}(1)\\cdot(\\psi_{1}^{7}\\cdot\\mathrm{ev}_{1}^{*}\\mathcal{O}_{\\mathbb{P}^{2}}(1)+\\psi_{1}^{6}\\cdot\\mathrm{ev}_{1}^{*}\\mathcal{O}_{\\mathbb{P}^{2}}(1)^{2}) &= -\\frac{5}{16} \\\\
\\end{aligned}
```
can be computed as
```julia-repl
julia> P = (g,c,w,s,m) -> O1_i(g,c,w,s,m,1)^5*O1_i(g,c,w,s,m,2)^2*Hypersurface(g,c,w,s,5)*Psi(g,c,w,s,m,[1,0]);
julia> AtiyahBottFormula(6,2,2,P);
Result: 495000//1
julia> P = (g,c,w,s,m) -> O1_i(g,c,w,s,m,1)^8*O1_i(g,c,w,s,m,2)^6*Hypersurface(g,c,w,s,7)*Psi(g,c,w,s,m,2);
julia> AtiyahBottFormula(10,2,2,P);
Result: 71804533752//1
julia> P = (g,c,w,s,m) -> O1(g,c,w,s,m)^2*Psi(g,c,w,s,m,4);
julia> AtiyahBottFormula(2,2,1,P);
Result: 1//8
julia> P = (g,c,w,s,m) -> Incidency(g,c,w,s,2)^4*O1_i(g,c,w,s,m,1)*(O1_i(g,c,w,s,m,1) + Psi(g,c,w,s,m,1))
julia> AtiyahBottFormula(2,2,1,P); #number of plane conics through four points and tangent to a line
Result: 2
julia> P = (g,c,w,s,m) -> O1(g,c,w,s,m)^2*(Psi(g,c,w,s,m,3)*O1(g,c,w,s,m)+Psi(g,c,w,s,m,2)*O1(g,c,w,s,m)^2);
julia> AtiyahBottFormula(2,2,1,P);
Result: 1//8
julia> P = (g,c,w,s,m) -> O1(g,c,w,s,m)*(Psi(g,c,w,s,m,7)*O1(g,c,w,s,m)+Psi(g,c,w,s,m,6)*O1(g,c,w,s,m)^2);
julia> AtiyahBottFormula(3,2,1,P);
Result: -5//16
```
!!! warning "Psi is singleton!"

    `Psi` cannot be multiplied by itself.
    ```julia-repl
    julia> P = (g,c,w,s,m) -> O1(g,c,w,s,m)^2*Psi(g,c,w,s,m,1)^4;                  #this is **wrong**
    julia> AtiyahBottFormula(2,2,1,P);
    Warning: more instances of Psi has been found. Type:
    julia> ?Psi
    for support.
    julia> P = (g,c,w,s,m) -> O1(g,c,w,s,m)^2*Psi(g,c,w,s,m,3)*Psi(g,c,w,s,m,1);   #this is **wrong**
    julia> AtiyahBottFormula(2,2,1,P);
    Warning: more instances of Psi has been found. Type:
    julia> ?Psi
    for support.
    julia> P = (g,c,w,s,m) -> O1(g,c,w,s,m)^2*Psi(g,c,w,s,m,4);
    julia> AtiyahBottFormula(2,2,1,P);
    Result: 1//8
    ```
"""
function Psi(g::SimpleGraph, col::Vector{UInt8}, weights::Vector{Int64}, scalars::Vector{fmpq}, mark::Marks, a::Int64)::fmpq

    return Psi(g, col, weights, scalars, mark, [a])
end
function Psi(g::SimpleGraph, col::Vector{UInt8}, weights::Vector{Int64}, scalars::Vector{fmpq}, mark::Marks, a::Vector{Int64})::fmpq

    if findfirst(x -> x>0, a) === nothing #if all of them are zero or a is empty
        return fmpq(1)
    end

    local q1 = fmpq(1)

    local inv_marks::Dict{Int64,Vector{Int64}} = invert_marks(mark)

    for v in vertices(g)
        S_v = num_marks(mark, v)
        S_v == 0 && continue

        nghbrs = all_neighbors(g, v)
        N = S_v + length(nghbrs)
        d = Dict(edges(g).=> weights) #assign weights to edges

        if N == 2
            i = inv_marks[v][1]
            (i > length(a) || a[i] == 0) && continue
            w = nghbrs[1]
            e = SimpleEdge(v,w)
            d_e = haskey(d,e) ? d[e] : d[reverse(e)]
            q1 *= ((scalars[col[w]]-scalars[col[v]])//d_e)^a[i]
            continue
        end
        #If no previous condition holds, then N>2
        if length(a) < mark.m #correct the dimension of a, if it is too low
            a = vcat(a,zeros(Int64,mark.m-length(a)))
        end
        a_v = [a[i] for i in inv_marks[v]]
        Sum_ai = sum(a_v)
        if Sum_ai > N-3
            return fmpq(0)
        end

        local s1 = fmpq(0)

        for w in nghbrs
            e = (SimpleEdge(v,w) in edges(g)) ? SimpleEdge(v,w) : SimpleEdge(w,v)
            s1 += d[e]//(scalars[col[v]]-scalars[col[w]])
        end
        s1 ^= -Sum_ai
        q1 *= multinomial(N-3-Sum_ai,a_v...)*s1

    end

    return q1
end
"""
    Jet(g, c, w, s, m, p, q)

Equivariant class of the jet bundle ``J^p`` of the pull back of ``\\mathcal{O}_{\\mathbb{P}^n}(q)`` with respect to the first ``\\psi``-class.
# Arguments
- `g::SimpleGraph`: the graph.
- `c::Vector{UInt8}`: the coloration.
- `w::Vector{Int64}`: the weights.
- `s::Rational{BigInt}`: the scalars.
- `m::Marks`: the marks.
- `p::Int64`: the exponent of the Jet bundle. In particular, it is a bundle of rank p+1.
- `q::Int64`: the degree of the line bundle that is pulled back.

!!! note

    In order to define this bundle, the number of marks must be at least 1.
    You cannot multiply this bundle by the class `Psi(g,c,w,s,m,a)`.

# Example
```math
\\begin{aligned}
\\int_{\\overline{M}_{0,1}(\\mathbb{P}^{2},2)}\\delta_{*}(\\mathrm{ev}^{*}\\mathcal{O}_{\\mathbb{P}^{2}}(1)^{2})^{4}\\cdot\\mathrm{c_{top}}(J^{1}(\\mathrm{ev}_{1}^{*}\\mathcal{O}_{\\mathbb{P}^{2}}(1))) &= 2 \\\\
\\int_{\\overline{M}_{0,1}(\\mathbb{P}^{2},2)}\\delta_{*}(\\mathrm{ev}^{*}\\mathcal{O}_{\\mathbb{P}^{2}}(1)^{2})^{4}\\cdot(\\mathrm{c_{top}}(J^{1}(\\mathrm{ev}_{1}^{*}\\mathcal{O}_{\\mathbb{P}^{2}}(1)))+\\mathrm{ev}_{1}^{*}\\mathcal{O}_{\\mathbb{P}^{2}}(1)^{2}) &= 3 \\\\
\\int_{\\overline{M}_{0,1}(\\mathbb{P}^{3},d)}\\frac{\\mathrm{ev}_{1}^{*}\\mathcal{O}_{\\mathbb{P}^{3}}(1)^{2}}{k}\\cdot\\mathrm{c_{top}}(J^{4d-2}(\\mathrm{ev}_{1}^{*}\\mathcal{O}_{\\mathbb{P}^{3}}(k))) &= \\frac{(4d-2)!}{(d!)^{4}} \\\\
\\end{aligned}
```
can be computed as
```julia-repl
julia> P = (g,c,w,s,m) -> Incidency(g,c,w,s,2)^4*Jet(g,c,w,s,m,1,1);
julia> AtiyahBottFormula(2,2,1,P);
Result: 2//1
julia> P = (g,c,w,s,m) -> Incidency(g,c,w,s,2)^4*(Jet(g,c,w,s,m,1,1)+O1(g,c,w,s,m)^2);
julia> AtiyahBottFormula(2,2,1,P);
Result: 3//1
julia> P = (g,c,w,s,m) -> (O1(g,c,w,s,m)^2)//k*Jet(g,c,w,s,m,4*d-2,k);
julia> d=1;k=1;AtiyahBottFormula(3,d,1,P);   #The value of this integral does not depend on k, only on d
```
"""
function Jet(g::SimpleGraph, col::Vector{UInt8}, weights::Vector{Int64}, scalars::Vector{fmpq}, mark::Marks, p::Int64, q::Int64)::fmpq

    local s1 = fmpq(0)
    for h in 0:p
        s1 += stirlings1(p+1, p+1-h)*(q*O1_i(g, col, weights, scalars, mark, 1))^(p+1-h)*Psi(g, col, weights, scalars, mark, [h])
    end

    return s1
end

"""
    Euler_inv(g, c, w, s, m)

The inverse of the (equivariant) Euler class of the normal bundle. This function is invoked automatically.
# Arguments
- `g::SimpleGraph`: the graph.
- `c::Vector{UInt8}`: the coloration.
- `w::Vector{Int64}`: the weights.
- `s::Rational{BigInt}`: the scalars.
- `m::Marks`: the marks.

"""
function Euler_inv(g::SimpleGraph, col::Vector{UInt8}, weights::Vector{Int64}, scalars::Vector{fmpq}, mark::Marks)::fmpq

    local V = fmpq(1)
    local E = fmpq(1)

    # col = Dict(vertices(g).=> coloration) #assing colors to vertices
    d = Dict(edges(g).=> weights) #assign weights to edges

    omega_inv = Dict(edges(g).=> [d[e]//(scalars[col[src(e)]]-scalars[col[dst(e)]]) for e in edges(g)])
    merge!(omega_inv,Dict(reverse.(edges(g)).=> [d[e]//(scalars[col[dst(e)]]-scalars[col[src(e)]]) for e in edges(g)]))

    max_col = length(scalars)

    for e in edges(g)
        local q1 = fmpq(1)
        for j in 1:max_col
            if j != col[src(e)] && j != col[dst(e)]
                for alph in 0:d[e]
                    a = fmpq( 1 )
                    mul!( a, scalars[col[src(e)]], fmpq( alph ))
                    b = fmpq( 1 )
                    mul!( b, scalars[col[dst(e)]], fmpq( d[e]-alph ))
                    addeq!( a, b )
                    mul!( a, a, inv( fmpq( d[e] )))
                    addeq!( a, -fmpq(scalars[j]))
                    mul!( q1, q1, a )
                end
            end
        end
        mul!( E, E, omega_inv[e]^(2*d[e])*((-1)^d[e]))
        mul!( E, E, inv(fmpq(factorial(d[e])^2)))
        mul!( E, E, inv( fmpq( q1 )))
    end

    for v in vertices(g)
        nghbrs = all_neighbors(g, v)
        local p1 = fmpq(1)
        for j in 1:max_col
            if j != col[v]
                a = deepcopy( scalars[col[v]] )
                addeq!( a, -scalars[j] )
                mul!( p1, p1, a )
            end
        end
        p1 ^= length(nghbrs)-1

        local s1 = fmpq(0)
        local r1 = fmpq(1)

        for w in nghbrs
            e = SimpleEdge(v,w)
            s1 += omega_inv[e]
            r1 *= omega_inv[e]
        end
        s1 ^= length(nghbrs) + num_marks(mark,v) - 3
        mul!( V, V, p1 )
        mul!( V, V, s1 )
        mul!( V, V, r1 )
        #V *= p1*s1*r1

    end

    return V*E
end

#####################################################
###List of the class useful for computing the rank###
#####################################################
function Incidency(n::Int64, deg::Int64, c::Int64, s::Int64, r::Int64)::Cycle

    try
        r < 0 && error(string("Incidency requires a positive integer, correct ", r))
    catch e
        printstyled(stderr,"ERROR: ", bold=true, color=:red)
        printstyled(stderr,sprint(showerror,e), color=:light_red)
        println(stderr)
        return error_cycle()
    end

    return Cycle(r-1, 0)
end

function Incidency(n::Int64, deg::Int64, c::Int64, s::Int64, r::Vector{Int64})::Cycle

    try
        if findfirst(x -> x<1, r) !== nothing #if some of them is non positive
            error(string("Incidency requires a positive integer, correct ", r))
        end
    catch e
        printstyled(stderr,"ERROR: ", bold=true, color=:red)
        printstyled(stderr,sprint(showerror,e), color=:light_red)
        println(stderr)
        return error_cycle()
    end

    return Cycle(sum(r)-length(r), 0)
end

function Hypersurface(n::Int64, deg::Int64, c::Int64, s::Int64, b::Int64)::Cycle

    try
        b < 1 && error(string("Hypersurface requires a positive integer, correct ", b))
    catch e
        printstyled(stderr,"ERROR: ", bold=true, color=:red)
        printstyled(stderr,sprint(showerror,e), color=:light_red)
        println(stderr)
        return error_cycle()
    end

    return Cycle(b*deg+1, 0)
end

function Hypersurface(n::Int64, deg::Int64, c::Int64, s::Int64, b::Vector{Int64})::Cycle

    try
        if findfirst(x -> x<1, b) !== nothing #if some of them is negative
            error(string("Hypersurface requires a positive integer, correct ", b))
        end
    catch e
        printstyled(stderr,"ERROR: ", bold=true, color=:red)
        printstyled(stderr,sprint(showerror,e), color=:light_red)
        println(stderr)
        return error_cycle()
    end

    return Cycle(sum(b)*deg+length(b), 0)
end

function Contact(n::Int64, deg::Int64, c::Int64, s::Int64)::Cycle

    return Cycle(2*deg-1, 0)
end

function O1_i(n::Int64, deg::Int64, n_marks::Int64, c::Int64, s::Int64, i::Int64)::Cycle

    try
        if (i < 1 || i > n_marks)
            error(string("O1_i requires a positive integer between 1 and ", n_marks, ", correct ",i))
        end
    catch e
        printstyled(stderr,"ERROR: ", bold=true, color=:red)
        printstyled(stderr,sprint(showerror,e), color=:light_red)
        println(stderr)
        return error_cycle()
    end

    return Cycle(1, 0)
end

function O1(n::Int64, deg::Int64, n_marks::Int64, c::Int64, s::Int64)::Cycle

    return Cycle(n_marks, 0)
end

function R1(n::Int64, deg::Int64, c::Int64, s::Int64, k::Int64)::Cycle

    try
        k < 1 && error(string("R1 requires a positive integer, correct ", k))
    catch e
        printstyled(stderr,"ERROR: ", bold=true, color=:red)
        printstyled(stderr,sprint(showerror,e), color=:light_red)
        println(stderr)
        return error_cycle()
    end

    return Cycle(k*deg-1, 0)
end

function Psi(n::Int64, deg::Int64, n_marks::Int64, c::Int64, s::Int64, a::Int64)::Cycle

    try
        a < 0 && error(string("exponents of psi classes must be nonnegative, correct ", a))
    catch e
        printstyled(stderr,"ERROR: ", bold=true, color=:red)
        printstyled(stderr,sprint(showerror,e), color=:light_red)
        println(stderr)
        return error_cycle()
    end
    psi_deg = a==0 ? 0 : 1

    return Cycle(a, psi_deg)
end

function Psi(n::Int64, deg::Int64, n_marks::Int64, c::Int64, s::Int64, a::Vector{Int64})::Cycle

    try
        length(a) > n_marks && error(string("size of ",a," is greater than ",n_marks))

        if findfirst(x -> x<0, a) !== nothing #if some of them is negative
            error(string("exponents of psi classes must be nonnegative, correct ", a))
        end
    catch e
        printstyled(stderr,"ERROR: ", bold=true, color=:red)
        printstyled(stderr,sprint(showerror,e), color=:light_red)
        println(stderr)
        return error_cycle()
    end
    psi_deg = sum(a)==0 ? 0 : 1
    return Cycle(sum(a), psi_deg)
end

function Jet(n::Int64, deg::Int64, n_marks::Int64, c::Int64, s::Int64, p::Int64, q::Int64)::Cycle

    try
        1 > n_marks && error("Jet is defined when there is at least one mark")
        0 > p       && error(string("exponents of Jet must be nonnegative, correct ", p))
    catch e
        printstyled(stderr,"ERROR: ", bold=true, color=:red)
        printstyled(stderr,sprint(showerror,e), color=:light_red)
        println(stderr)
        return error_cycle()
    end
    psi_deg = p==0 ? 0 : 1
    return Cycle(p+1, psi_deg)
end


text:
function LogKernelAsymptotics(S::Fourier,k)
    d=domain(S)
    r=d.radius
    if k==-1  # coefficient in front of logabs term
        FiniteOperator([2r].',S,ConstantSpace())
    else
        error("Not implemented")
    end
end

LogKernelAsymptotics(S::Laurent,k)=LogKernelAsymptotics(Fourier(domain(S)),k)*Conversion(S,Fourier(domain(S)))

function LogKernelAsymptotics{CD<:ChebyshevDirichlet{1,1}}(S::JacobiWeight{CD},k)
    @assert S.α==S.β==-0.5
    d=domain(S)
    if k==-1  # coefficient in front of logabs term
        FiniteOperator([arclength(d)/2,0,-arclength(d)/2].',S)
    else
        error("Not implemented")
    end
end

function LogKernelAsymptotics{C<:Chebyshev}(S::JacobiWeight{C},k)
    @assert S.α==S.β==-0.5
    d=domain(S)
    if k==-1  # coefficient in front of logabs term
        FiniteOperator([arclength(d)/2].',S,ConstantSpace())
    else
        error("Not implemented")
    end
end


text:
module Destruct
export destruct

using Base.Cartesian

unpack_broadcast(w::Array{<:Tuple}) = Tuple((v->v[i]).(w) for i=1:length(w[1]))

"""
    destruct(v::Array{<:Tuple,N})

Destructure an array of tuples to a tuple of arrays. Works for tuples with
elements of varying types.
## Examples
```julia-repl
julia> f(a, b) = a+b, a*b, a-b;
julia> v = f.(rand(3,1), rand(1,4));
julia> x, y, z = destruct(v);
julia> x
3×4 Array{Float64,2}:
 0.301013  0.888299  1.03866  1.0867
 0.853248  1.44053   1.5909   1.63894
 0.687546  1.27483   1.4252   1.47324
julia> v = f.(rand(100,1,1), rand(1,100,100));
julia> @btime destruct(v);
  7.138 ms (7 allocations: 22.89 MiB)
julia> x, y, z = f.(rand(100,1,1), rand(1,100,100)) |> destruct;
```
"""
@generated function destruct(v::Array{T,N}) where {T <: Tuple, N}
    TT = T.types
    M = length(TT)
    quote
        @nexprs $M i -> out_i = similar(v,$TT[i])
        @inbounds for j in eachindex(v)
            @nexprs $M i -> (out_i[j] = v[j][i])
        end
        return @ntuple $M out
    end
end

@generated function destruct_cartesian(v::Array{T,N}) where {T <: Tuple, N}
    TT = T.types
    M = length(TT)
    quote
        @nexprs $M i -> out_i = similar(v,$TT[i])
        @inbounds @nloops $N j v begin
            @nexprs $M i -> ((@nref $N out_i j) = (@nref $N v j)[i])
        end
        return @ntuple $M out
    end
end

end # module Destruct


text:
# merge multiple JSON objects garnered via `--meta-json`
function merge_json_objects(objs::Vector)
    merged = Dict()
    for obj in objs
        for k in keys(obj)
            if !haskey(merged, k)
                merged[k] = obj[k]
            else
                if merged[k] != obj[k]
                    if !isa(merged[k], Array)
                        merged[k] = [merged[k]]
                    end

                    if isa(obj[k], Array)
                        append!(merged[k], obj[k])
                    else
                        push!(merged[k], obj[k])
                    end
                    merged[k] = unique(merged[k])
                end
            end
        end
    end
    return merged
end

# Fix various things that go wrong with the JSON translation
function cleanup_merged_object!(meta::Dict)
    # Turn `source` values back into AbstractSource's
    meta["sources"] = sourcify.(meta["sources"])

    # Turn `dependencies` back into AbstractDependency's
    meta["dependencies"] = dependencify.(meta["dependencies"])

    # Turn `version` back into a VersionNumber (not a string)
    meta["version"] = VersionNumber(meta["version"])

    # Reconstruct our `Product`s
    function reconstruct_product(p::Dict)
        if p["type"] == "exe"
            return ExecutableProduct(p)
        elseif p["type"] == "lib"
            return LibraryProduct(p)
        elseif p["type"] == "framework"
            return FrameworkProduct(p)
        elseif p["type"] == "file"
            return FileProduct(p)
        else
            error("Unknown product type $(p["type"])")
        end
    end
    meta["products"] = Product[reconstruct_product(p) for p in meta["products"]]

    if haskey(meta, "platforms")
        # Convert platforms back to actual Platform objects
        meta["platforms"] = [platform_key_abi(p) for p in meta["platforms"]]
    else
        # If the key isn't there it's because this is a platform-independent
        # build, so use `AnyPlatform()`.
        meta["platforms"] = [AnyPlatform()]
    end

    # Return the cleaned-up meta for fun
    return meta
end


text:
using BenchmarkTools
using NiGraphEmbedding
using ForwardDiff, NiLang, Zygote, ReverseDiff

function get_hessian_reversediff(params0::AbstractArray{T}) where T
    N = length(params0)
    params = Dual.(params0, zero(T))
    hes = zeros(T, N, N)
    for i=1:N
        @inbounds i !== 1 && (params[i-1] = Dual(params0[i-1], zero(T)))
        @inbounds params[i] = Dual(params0[i], one(T))
        res = ReverseDiff.gradient(embedding_loss, params)
        hes[:,i] .= vec(ForwardDiff.partials.(res, 1))
    end
    hes
end

suite = BenchmarkGroup()
suite["NiLang"] = BenchmarkGroup(["Term"])
suite["ForwardDiff"] = BenchmarkGroup(["Term"])
suite["Zygote"] = BenchmarkGroup(["Term"])
suite["ReverseDiff"] = BenchmarkGroup(["Term"])

cases = [("NiLang", "Call"), ("NiLang", "Uncall"),
         ("NiLang", "Gradient"), ("NiLang", "Hessian"),
         ("ForwardDiff", "Call"), ("ForwardDiff", "Gradient"),
         ("ForwardDiff", "Hessian"),
         ("Zygote", "Gradient"),
         ("ReverseDiff", "Gradient"), ("ReverseDiff", "Hessian")
         ]

for (lang, term) in cases
    suite[lang][term] = BenchmarkGroup(["dimension"])
end

for k=1:10
    suite["NiLang"]["Call"][k] = @benchmarkable embedding_loss(0.0, $(randn(k, 10)))
    suite["NiLang"]["Uncall"][k] = @benchmarkable (~embedding_loss)(0.0, $(randn(k, 10)))
    suite["NiLang"]["Gradient"][k] = @benchmarkable NiLang.AD.gradient(Val(1), embedding_loss, (0.0, $(randn(k, 10))))
    suite["NiLang"]["Hessian"][k] = @benchmarkable get_hessian($(randn(k, 10)))
    suite["ForwardDiff"]["Call"][k] = @benchmarkable embedding_loss($(randn(k, 10)))
    suite["ForwardDiff"]["Gradient"][k] = @benchmarkable ForwardDiff.gradient(embedding_loss, $(randn(k, 10)))
    suite["ForwardDiff"]["Hessian"][k] = @benchmarkable ForwardDiff.hessian(embedding_loss, $(randn(k, 10)))
    suite["Zygote"]["Gradient"][k] = @benchmarkable Zygote.gradient(embedding_loss, $(randn(k, 10)))
    suite["ReverseDiff"]["Gradient"][k] = @benchmarkable ReverseDiff.gradient!(ctape, $((randn(k, 10),))) setup=(ctape = ReverseDiff.compile(ReverseDiff.GradientTape(embedding_loss, (randn(k, 10),))))
    suite["ReverseDiff"]["Hessian"][k] = @benchmarkable get_hessian_reversediff($(randn(k, 10)))
    #suite["ReverseDiff"]["Hessian"][k] = @benchmarkable ReverseDiff.hessian(embedding_loss, $(randn(k, 10)))
end

tune!(suite)
res = run(suite)#; seconds=100, samples=1000)

function analyze_res(res)
    times = zeros(10, length(cases))
    for (k, (lang, term)) in enumerate(cases)
        for i=1:10
            @show lang, term
            times[i,k] = minimum(res[lang][term][i].times)
        end
    end
    return times
end

times = analyze_res(res)
using DelimitedFiles
writedlm(joinpath(dirname(@__FILE__), "bench_graphembedding.dat"), times)


text:
function AstronomicalAmplitudeReader(arg0::JString, arg1::jint, arg2::jint, arg3::jint, arg4::jdouble)
    return AstronomicalAmplitudeReader((JString, jint, jint, jint, jdouble), arg0, arg1, arg2, arg3, arg4)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_astronomical_amplitudes_map(obj::AstronomicalAmplitudeReader)
    return jcall(obj, "getAstronomicalAmplitudesMap", Map, ())
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_supported_names(obj::AstronomicalAmplitudeReader)
    return jcall(obj, "getSupportedNames", JString, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function load_data(obj::AstronomicalAmplitudeReader, arg0::InputStream, arg1::JString)
    return jcall(obj, "loadData", void, (InputStream, JString), arg0, arg1)
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function still_accepts_data(obj::AstronomicalAmplitudeReader)
    return jcall(obj, "stillAcceptsData", jboolean, ())
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
csconv = Dict{Type, AbstractVector{Color}}(
    Lab{Float64} => Lab{Float64}[
Lab{Float64}(61.40554408859231,-61.26004899789989,57.97749162785381),
Lab{Float64}(55.54652888796126,30.809120864024266,53.44054155482845),
Lab{Float64}(44.15905874913364,-27.770886035217714,34.72162815967184),
Lab{Float64}(62.25658951557699,-14.029208420051065,52.59616361156139),
Lab{Float64}(52.41763760855498,78.65495176055043,28.797415938527948),
Lab{Float64}(37.36754271570359,12.494669958826199,4.788494682968725),
Lab{Float64}(75.82115093728481,-24.72761078970176,64.60074310217375),
Lab{Float64}(71.8426650650214,-7.604216749949655,57.77684960484335),
Lab{Float64}(37.754274173620324,70.3180945761324,-83.16037630169984),
Lab{Float64}(42.15029927306029,68.36754602428618,-47.55058470357214),
Lab{Float64}(19.332795440287164,42.919076249245904,-62.726071239338566),
Lab{Float64}(50.02974879505459,76.5885392906267,45.398569298282474),
Lab{Float64}(40.561202169048386,-31.962174166736716,24.242059437675913),
Lab{Float64}(60.757760020025074,-52.783977694143616,58.95267093429221),
Lab{Float64}(77.9139415958467,16.03655624480843,-21.775158356104264),
Lab{Float64}(57.14306499534186,6.570139697934629,-24.826155752006727),
Lab{Float64}(44.97666070400284,70.24363046832543,-20.532082159114996),
Lab{Float64}(60.664366637100485,-35.19163329720709,35.983642777064205),
Lab{Float64}(17.93658797983788,19.195613517492216,1.4523602503233035),
Lab{Float64}(22.70092488166712,42.09789731817473,-67.66437344882776)],
    RGB{Float64} => RGB{Float64}[
RGB{Float64}(0.11228396215602665,0.6732949141086532,0.09859409271099773),
RGB{Float64}(0.7844278907439127,0.42956445925185305,0.14202724843192838),
RGB{Float64}(0.2915703179069309,0.45173923490473555,0.16825129849152776),
RGB{Float64}(0.6035334387508868,0.609857169290216,0.19624892225695534),
RGB{Float64}(0.9644226672159013,0.059403296858277166,0.3174797751494354),
RGB{Float64}(0.432484993970361,0.31581143833886605,0.31642899914629097),
RGB{Float64}(0.6947521946179401,0.7727221308629347,0.21756330942282065),
RGB{Float64}(0.7580992752894541,0.6964329890060948,0.2510418460157361),
RGB{Float64}(0.4326423261471005,0.12533218110917999,0.892922405285977),
RGB{Float64}(0.6554980398035374,0.13015145038942355,0.703223224556291),
RGB{Float64}(0.10975505072305158,0.09235233871920157,0.5557379427193866),
RGB{Float64}(0.9337598483853324,0.00782029356933489,0.18443366775765324),
RGB{Float64}(0.1874608064216865,0.4231093851641665,0.21203084253232074),
RGB{Float64}(0.2661788389076388,0.6555731749379143,0.0773468081126768),
RGB{Float64}(0.7978302602115971,0.7229225246920243,0.9138004107795679),
RGB{Float64}(0.48760377543392364,0.5312092935824387,0.7067152176969306),
RGB{Float64}(0.7689582640588691,0.10881793827304542,0.5577889667548762),
RGB{Float64}(0.39600674521787516,0.6317663735284889,0.3134289899365913),
RGB{Float64}(0.2729387941236906,0.1264998653293029,0.16620356290215066),
RGB{Float64}(0.03723463407032801,0.13430330431357462,0.6224972592798952)],
    HSI{Float64} => HSI{Float64}[
HSI{Float64}(118.8039444580707,0.6654701189568601,0.2947243229918925),
HSI{Float64}(26.537330142953778,0.6857849652100698,0.4520065328092315),
HSI{Float64}(94.29190735403657,0.4462751501965313,0.30385361710106473),
HSI{Float64}(60.76444178490796,0.5823423264482519,0.4698798434326861),
HSI{Float64}(343.93235717256283,0.8671370103299607,0.44710191307453795),
HSI{Float64}(359.7366651138475,0.110160904373796,0.3549084771518394),
HSI{Float64}(67.45289969116865,0.6126555782798495,0.5616792116345651),
HSI{Float64}(53.60149020968084,0.5584328271084302,0.5685247034370949),
HSI{Float64}(263.43653587702016,0.740852337559448,0.4836323041807525),
HSI{Float64}(295.6962678920551,0.7377516913969178,0.4962909049164173),
HSI{Float64}(241.89845939343923,0.6344148279341627,0.2526151107205466),
HSI{Float64}(349.6518042121588,0.9791646598774854,0.3753379365707736),
HSI{Float64}(125.4417442482248,0.31633635755402123,0.2742003447060579),
HSI{Float64}(101.32416888416884,0.7677502773116757,0.3330329406527433),
HSI{Float64}(262.9200122451308,0.10917223828930589,0.8115177318943965),
HSI{Float64}(229.16502448773292,0.1522530592134964,0.5751760955710976),
HSI{Float64}(318.2512056299271,0.7725956147105435,0.47852172302893026),
HSI{Float64}(105.5262070776537,0.2989222401885988,0.4470673695609852),
HSI{Float64}(344.80359875018354,0.32908191611333415,0.18854740745171472),
HSI{Float64}(231.09848635556082,0.8593212208480958,0.26467839922126596)],
    YIQ{Float64} => YIQ{Float64}[
YIQ{Float64}(0.4400367458354851,-0.14957309026873586,-0.29743867192165413),
YIQ{Float64}(0.5028893832345075,0.30377289091436527,-0.0144248883188723),
YIQ{Float64}(0.37153110397128625,-0.004341001642347811,-0.1220716976086112),
YIQ{Float64}(0.5608150336971648,0.1291098788046794,-0.13002519270561796),
YIQ{Float64}(0.3594248071203988,0.4562240955867101,0.2716684010514613),
YIQ{Float64}(0.3507672334047295,0.0693059044288959,0.024863468161431534),
YIQ{Float64}(0.6861210142815083,0.13190404991219223,-0.18921655074739152),
YIQ{Float64}(0.6640966183039183,0.17982328817008927,-0.12553706704197126),
YIQ{Float64}(0.3047232000316731,-0.06352876782830047,0.3038067584283862),
YIQ{Float64}(0.35256026287926645,0.12885061147524787,0.28939037486655916),
YIQ{Float64}(0.15038170846437382,-0.13850157531382595,0.14785538777008367),
YIQ{Float64}(0.30481014511678645,0.49485766540490517,0.25074607668125404),
YIQ{Float64}(0.32858750626013455,-0.07256790289262055,-0.1155032282283292),
YIQ{Float64}(0.47322646264678486,-0.04620569899723531,-0.26224622936178826),
YIQ{Float64}(0.7670800166263565,-0.016698265745364704,0.0752284812098746),
YIQ{Float64}(0.5381789190050847,-0.08236006454815478,0.04538538725375002),
YIQ{Float64}(0.3573825929299354,0.2490183747926694,0.27928123367606156),
YIQ{Float64}(0.524983782934139,-0.03817575987380177,-0.1488986898079142),
YIQ{Float64}(0.17481132656212944,0.07448068391233167,0.04331860007147402),
YIQ{Float64}(0.16093388277700443,-0.21466401451697095,0.1313684734434902)],
    LCHab{Float64} => LCHab{Float64}[
LCHab{Float64}(61.40554408859231,84.34561718716007,136.57692947148365),
LCHab{Float64}(55.54652888796126,61.68543920640755,60.03602910879027),
LCHab{Float64}(44.15905874913364,44.46137169768345,128.65335148070065),
LCHab{Float64}(62.25658951557699,54.43505410622245,104.93504261794598),
LCHab{Float64}(52.41763760855498,83.76092526465487,20.108859052334054),
LCHab{Float64}(37.36754271570359,13.380824291082131,20.96898020598853),
LCHab{Float64}(75.82115093728481,69.17160360089987,110.94564405157224),
LCHab{Float64}(71.8426650650214,58.27511014696499,97.49780692927679),
LCHab{Float64}(37.754274173620324,108.90492464281962,310.2169286804031),
LCHab{Float64}(42.15029927306029,83.27772484304839,325.1808465885591),
LCHab{Float64}(19.332795440287164,76.00399410038371,304.38111497800685),
LCHab{Float64}(50.02974879505459,89.03277174727756,30.65773681622481),
LCHab{Float64}(40.561202169048386,40.11555836885015,142.82110494046856),
LCHab{Float64}(60.757760020025074,79.13005567736441,131.8400549984751),
LCHab{Float64}(77.9139415958467,27.04308890689674,306.3701994798936),
LCHab{Float64}(57.14306499534186,25.680824462490943,284.8232668592222),
LCHab{Float64}(44.97666070400284,73.18288064267016,343.7064570495635),
LCHab{Float64}(60.664366637100485,50.33163619069481,134.3624622173942),
LCHab{Float64}(17.93658797983788,19.25047865923466,4.326814913733189),
LCHab{Float64}(22.70092488166712,79.69128178686806,301.8881183495275)],
    LMS{Float64} => LMS{Float64}[
LMS{Float64}(0.2307455933908811,0.3966129226546093,0.06205539564220773),
LMS{Float64}(0.3107019890330763,0.1897841714619035,0.0498207460521362),
LMS{Float64}(0.1218228377056502,0.17064451545425868,0.04614786486182357),
LMS{Float64}(0.30782939233464807,0.3408975716081407,0.0796305359309192),
LMS{Float64}(0.36293517359571514,0.06999350791647939,0.09885173109323558),
LMS{Float64}(0.10665192537186992,0.08953048173500038,0.09046458812156158),
LMS{Float64}(0.47953245310048775,0.5692794259137218,0.11805172223345846),
LMS{Float64}(0.4529818748418708,0.46457163956456615,0.11708679460289854),
LMS{Float64}(0.07608222895850802,0.02604919712027899,0.7296345053537655),
LMS{Float64}(0.1636469929148561,0.04282906750024994,0.43464808179574577),
LMS{Float64}(0.011814257293131639,0.00973850643965921,0.253539933890513),
LMS{Float64}(0.3347188452051175,0.06123620787020029,0.04649844014616222),
LMS{Float64}(0.09407290872609109,0.1462548602213991,0.054453038775816975),
LMS{Float64}(0.23553355846222834,0.3771163500455532,0.05733013093653634),
LMS{Float64}(0.5062615568672968,0.5071983572055632,0.838514346078231),
LMS{Float64}(0.2174546897679178,0.2500580133905722,0.4643617224024501),
LMS{Float64}(0.22442774528787113,0.05052631522063127,0.2684043405961488),
LMS{Float64}(0.2477609356957181,0.35305569878781684,0.12364967459867356),
LMS{Float64}(0.031897534250467936,0.01841899562603634,0.025306259760502797),
LMS{Float64}(0.013112766066299388,0.016272875996157916,0.32555480144136284)],
    HSL{Float64} => HSL{Float64}[
HSL{Float64}(118.57074822913228,0.7445381606943033,0.38594450340982545),
HSL{Float64}(26.855877022639817,0.6933963827794761,0.46322756958792055),
HSL{Float64}(93.89962599981922,0.4572455886709776,0.30999526669813166),
HSL{Float64}(60.917350742112355,0.5130940596657895,0.4030530457735857),
HSL{Float64}(342.89032345092164,0.9271086272022048,0.5119129820370892),
HSL{Float64}(359.6824160518213,0.1559188987062825,0.3741482161546135),
HSL{Float64}(68.42677085192335,0.5606048507387105,0.4951427201428777),
HSL{Float64}(52.703041187462325,0.5117352633925326,0.5045705606525951),
HSL{Float64}(264.0214220055365,0.7818628063239984,0.5091272931975785),
HSL{Float64}(295.0032243529561,0.6876520146285906,0.4166873374728573),
HSL{Float64}(242.25333439627224,0.7150016244212027,0.3240451407192941),
HSL{Float64}(348.5556218047031,0.9833890006365249,0.47079007097733366),
HSL{Float64}(126.25593489468521,0.3859483839694542,0.3052850957929265),
HSL{Float64}(100.40573294866358,0.7889351910129655,0.36645999152529557),
HSL{Float64}(263.54627989285797,0.5254333530120703,0.8183614677357961),
HSL{Float64}(228.05935891850672,0.27195805833187997,0.5971594965654271),
HSL{Float64}(319.1931281022064,0.7520599488025355,0.4388881011659573),
HSL{Float64}(104.43580436274233,0.3367953291950933,0.4725976817325401),
HSL{Float64}(343.732318489466,0.3666118071669045,0.19971932972649675),
HSL{Float64}(230.04870640337006,0.8871219219636491,0.3298659466751116)],
    HSV{Float64} => HSV{Float64}[
HSV{Float64}(118.57074822913228,0.8535647742987599,0.6732949141086532),
HSV{Float64}(26.855877022639817,0.8189416132345362,0.7844278907439127),
HSV{Float64}(93.89962599981922,0.627547741061258,0.45173923490473555),
HSV{Float64}(60.917350742112355,0.678205107459243,0.609857169290216),
HSV{Float64}(342.89032345092164,0.9384053290350767,0.9644226672159013),
HSV{Float64}(359.6824160518213,0.26977480665951337,0.432484993970361),
HSV{Float64}(68.42677085192335,0.7184456084105452,0.7727221308629347),
HSV{Float64}(52.703041187462325,0.6688535997875945,0.7580992752894541),
HSV{Float64}(264.0214220055365,0.8596382167507156,0.892922405285977),
HSV{Float64}(295.0032243529561,0.8149215699303098,0.703223224556291),
HSV{Float64}(242.25333439627224,0.8338203465696532,0.5557379427193866),
HSV{Float64}(348.5556218047031,0.9916249412706513,0.9337598483853324),
HSV{Float64}(126.25593489468521,0.5569448161757232,0.4231093851641665),
HSV{Float64}(100.40573294866358,0.8820165146018949,0.6555731749379143),
HSV{Float64}(263.54627989285797,0.20888356345200657,0.9138004107795679),
HSV{Float64}(228.05935891850672,0.31004206047388555,0.7067152176969306),
HSV{Float64}(319.1931281022064,0.8584865481532628,0.7689582640588691),
HSV{Float64}(104.43580436274233,0.503884658839859,0.6317663735284889),
HSV{Float64}(343.732318489466,0.5365266204262058,0.2729387941236906),
HSV{Float64}(230.04870640337006,0.940185063443651,0.6224972592798952)],
    Luv{Float64} => Luv{Float64}[
Luv{Float64}(61.40554408859231,-55.412996764588925,72.15094916974675),
Luv{Float64}(55.54652888796126,73.81632939424085,47.28648729212383),
Luv{Float64}(44.15905874913364,-20.1916300181127,41.666371386767686),
Luv{Float64}(62.25658951557699,3.2697611242111053,60.31219107239103),
Luv{Float64}(52.41763760855498,152.259323372186,15.212166995935663),
Luv{Float64}(37.36754271570359,18.436507756805966,3.8578775297180368),
Luv{Float64}(75.82115093728481,-7.409794283217813,77.29226601764863),
Luv{Float64}(71.8426650650214,15.684580802243243,66.90904075772004),
Luv{Float64}(37.754274173620324,7.754363614644651,-117.7874863449133),
Luv{Float64}(42.15029927306029,48.20313623793672,-76.81447399397281),
Luv{Float64}(19.332795440287164,-4.309127046305885,-66.66113003627298),
Luv{Float64}(50.02974879505459,157.6413223995605,26.989718351905513),
Luv{Float64}(40.561202169048386,-27.00395794142001,32.04972428697958),
Luv{Float64}(60.757760020025074,-45.680120676009444,71.02883623327125),
Luv{Float64}(77.9139415958467,7.2620670193642765,-36.88967721685302),
Luv{Float64}(57.14306499534186,-7.94128077035882,-38.54634624531412),
Luv{Float64}(44.97666070400284,85.27844273261951,-40.087521394485464),
Luv{Float64}(60.664366637100485,-29.74711054177352,49.71477206711368),
Luv{Float64}(17.93658797983788,20.000790452137903,-1.0590760936318002),
Luv{Float64}(22.70092488166712,-7.783710802693105,-77.2290167510072)],
    XYZ{Float64} => XYZ{Float64}[
XYZ{Float64}(0.153663016662531,0.29712708024737594,0.05852499318777285),
XYZ{Float64}(0.29674743038157375,0.23463442881668406,0.04651156755042157),
XYZ{Float64}(0.09437874875578477,0.13948552143012788,0.044709910030614676),
XYZ{Float64}(0.25690546129316694,0.3070356352157358,0.07594481890187685),
XYZ{Float64}(0.39636754695382753,0.20517774641604547,0.09647367409102692),
XYZ{Float64}(0.10846840739164045,0.0973772826964292,0.0903140653393484),
XYZ{Float64}(0.3884459435137322,0.49596865344466834,0.11200041764904389),
XYZ{Float64}(0.38836664204882565,0.4342527151240371,0.11187294869947646),
XYZ{Float64}(0.2094684185109957,0.09950961331508182,0.7399357020105221),
XYZ{Float64}(0.24686350958413994,0.1259742829096898,0.4394897712217018),
XYZ{Float64}(0.0565673219895524,0.028259210828662944,0.25725636227097265),
XYZ{Float64}(0.358313785231492,0.18443569110836341,0.04363959059527557),
XYZ{Float64}(0.0722806398419704,0.11592628720860365,0.05354850452537528),
XYZ{Float64}(0.1634846895954648,0.2897296286168442,0.05379230620150586),
XYZ{Float64}(0.5667179078583392,0.5306598592674642,0.8436009947819997),
XYZ{Float64}(0.2534837034751407,0.25069489921914667,0.46795995593110057),
XYZ{Float64}(0.28095997402550754,0.1452502685062495,0.27006920583932),
XYZ{Float64}(0.1957167557114015,0.28867334765309877,0.12114761724980394),
XYZ{Float64}(0.03445175720203697,0.025039782925762675,0.02528204539465763),
XYZ{Float64}(0.06932878517091455,0.03713556742882635,0.3303251691771589)],
    YCbCr{Float64} => YCbCr{Float64}[
YCbCr{Float64}(112.36804733797123,84.838038954416,75.63437414224273),
YCbCr{Float64}(126.13277492835714,82.38305926806305,172.9819070849848),
YCbCr{Float64}(97.36531176971168,102.30325567748676,115.22453057007604),
YCbCr{Float64}(138.81849237967913,81.91489437546983,134.82520279105893),
YCbCr{Float64}(94.71403275936733,122.6975484272026,224.66156450445874),
YCbCr{Float64}(92.81802411563575,123.65925642822798,141.056189978181),
YCbCr{Float64}(166.2605021276503,68.76924167895929,129.37902991427086),
YCbCr{Float64}(161.4371594085581,75.7853913624277,143.01897834216263),
YCbCr{Float64}(82.7343808069364,202.35470355580298,148.43784790109092),
YCbCr{Float64}(93.21069757055935,172.3275136666039,176.40088871970542),
YCbCr{Float64}(48.93359415369787,179.2414173424112,121.50899835317183),
YCbCr{Float64}(82.75342178057623,112.7829605557114,228.4883941419257),
YCbCr{Float64}(87.96066387096947,113.26601255596279,105.45194375633868),
YCbCr{Float64}(119.63659531964588,77.95658463450972,94.91964940996402),
YCbCr{Float64}(183.99052364117208,146.5470355623716,132.91301656099364),
YCbCr{Float64}(133.86118326211357,149.3048212702825,119.91951706554495),
YCbCr{Float64}(94.26678785165586,153.33343129623827,193.75815817524418),
YCbCr{Float64}(130.97144846257646,101.25721970896373,107.39311873395408),
YCbCr{Float64}(54.28368051710635,126.91186193651748,143.67799687737957),
YCbCr{Float64}(51.24452032816397,186.3466274854119,108.23634423699981)],
)


text:
#!/usr/bin/env julia

module Rasim

if !isinteractive()
    print("Loading modules... ")
end

using Params
using Agents.BaseAgent
using Agents
using Agents.ContextAware
using Traffic.Simple
using Channel.Gilbert
using HDF5
using JLD
using Movement

if verbose
    using Plots
end

using Base.Collections: heapify!, heappush!, heappop!
using DataStructures: Queue, Deque, enqueue!, dequeue!, front, back, DefaultDict
using Util

type Event
    t :: Float64
    a :: Action
end

# I'd define Base.<, but it was giving me all sorts of errors
Base.isless(e1 :: Event, e2 :: Event) = e1.t < e2.t # generated code for < is shorter for Floats

if !isinteractive()
    println("DONE!")
end

const agent_types = [ContextQ, OptHighestSNR, CooperativeQ, RandomChannel]

# generate ith channel
function genchan(i, P)
    goodness = shuffle!(append!(2 * ones(Int, P.n_good_channel), ones(Int, n_channel - P.n_good_channel)))
    GilbertChannel(base_freq + chan_bw * i, noise[goodness[i], :], chan_trans_prob)
end

function run_simulation{AgentT <: Agent}(:: Type{AgentT}, at_no :: Int, P :: ParamT, init_positions :: Array{Point{Float64}, 2})
    const output_dir = joinpath("data/", P.prefix)
    const n_agent = int(P.n_agent)
    avg_energies = zeros(Float64, n_agent, t_total)
    avg_bits = zeros(Float64, n_agent, t_total)
    en_idle = zeros(n_agent, t_total)
    en_sense = zeros(n_agent, t_total)
    en_sw = zeros(n_agent, t_total)
    en_tx = zeros(n_agent, t_total)
    buf_overflow = falses(n_runs, n_agent, t_total)
    buf_levels = zeros(Int16, n_runs, n_agent, t_total)
    trajectories = Array(Point{Float64}, t_total, n_agent) # trajectories for a single run
    transmission_channels = zeros(Int8, n_runs, n_agent, t_total) # channels tried for transmission
    # Q = zeros(n_agent, int(Params.n_channel), P.buf_levels + 1, Params.n_channel * length(Params.P_levels) + 1)
    expertness = zeros(n_agent)
    pkt_sent = zeros(Int, n_agent)
    interference = zeros(n_agent, n_channel)
    release_time = zeros(n_channel)
    generated_packets = zeros(Int64, n_agent, n_runs)
    tried_packets = zeros(Int64, n_agent, n_runs)
    sent_packets = zeros(Int64, n_agent, n_runs)
    init_distances = zeros(Float64, n_agent)
    final_distances = zeros(Float64, n_agent)
    latencyhist = Array(DefaultDict{Int, Int, Int}, n_agent, n_runs)
    lostinchan = 0
    for i=1:n_agent, j=1:n_runs
        latencyhist[i, j] = DefaultDict(Int, Int, 0)
    end
    energies = zeros(Float64, n_agent, t_total)

    channels = [genchan(i, P) for i=1:n_channel]
    traffics = [SimpleTraffic() for i=1:n_channel]
    env = Environment(channels, traffics)
    # coordinator for cooperative Q learning type of agents
    coordinator = initcoordinator(AgentT, P)

    println("Agent type: ", AgentT, " - ", P.iteration)

    for n_run=1:n_runs
        agents = [AgentT(i, P, init_positions[n_run, i]) for i::Int8=1:n_agent]
        fill!(energies, 0.0)
        bits_transmitted = zeros(n_agent, t_total)
        if verbose
            @printf("Run %d of %d (agent), %d of %d (total)\n", n_run, n_runs, n_run + (at_no - 1) * n_runs, length(agent_types) * n_runs)
        end

        rates = zeros(5)

        for a=agents
            init_distances[a.s.id] += (a.s.pos.x .^ 2 + a.s.pos.y .^ 2)
        end

        for t=1:t_total
            # iterate environment
            fill!(release_time, 0.)

            for traf in traffics
                iterate(traf)
            end

            for i=1:endof(channels)
                iterate!(channels[i])
                if traffics[i].traffic
                    # add interference from a stationary PU, we may change this in the future
                    interfere!(channels[i], Params.P_tx, Point(Params.r_init / 2, 0.))
                end
            end

            # get initial actions
            actions = [Event(Params.t_slot - a.s.t_remaining, initial_action(a, env, t, P)) for a in agents]
            heapify!(actions)
            # save initial statistics
            @inbounds trajectories[t, :] = [a.s.pos for a in agents]
            @simd for i=1:n_agent
                @inbounds buf_overflow[n_run, i, t] = agents[i].s.buf_overflow
                @inbounds buf_levels[n_run, i, t] = B - agents[i].s.B_empty
            end
            @inbounds generated_packets[:, n_run] += buf_levels[n_run, :, t]'
            if t > 1
                @inbounds generated_packets[:, n_run] -= buf_levels[n_run, :, t - 1]'
            end

            cooperate(agents, P, coordinator, t)

            # interpret actions
            last_actions = Array(Action, n_agent)
            # whether agents[i] has collided
            collided = zeros(Bool, n_agent)
            # resolve all actions
            fill!(pkt_sent, 0)
            @inbounds begin
                while ! isempty(actions)
                    const event = heappop!(actions)
                    const action = event.a
                    const tt = event.t
                    const agentid = int(action.i)
                    a = agents[agentid]
                    if isa(action, Switch)
                        switch!(agents[agentid], action.chan)
                        heappush!(actions, Event(Params.t_slot - a.s.t_remaining, act(a, env, t)))
                        last_actions[agentid] = action
                    elseif isa(action, Sense)
                        heappush!(actions, Event(Params.t_slot - a.s.t_remaining, act(a, env, t)))
                        last_actions[agentid] = action
                    elseif isa(action, Transmit)
                        i = int(action.chan)
                        traffics[i].occupancy = min(tt, traffics[i].occupancy)
                        # add interference
                        interfere!(channels[i], action.power, a.s.pos)
                        interference[:, i] = max(interference[:, i], channels[i].interference)
                        # keep current interference
                        interference[agentid, i] = channels[i].interference
                        # use previous occupier data to mark as collided
                        if traffics[i].occupier > -1 && traffics[i].occupier != agentid
                            if traffics[i].occupier > 0
                                collided[traffics[i].occupier] = true # also mark the occupier
                                if release_time[i] > tt
                                    collided[agentid] = true
                                end
                            else
                                collided[agentid] = true
                            end
                        end
                        if traffics[i].occupier != 0
                            # mark self as occupier
                            traffics[i].occupier = agentid
                        end
                        release_time[i] = max(release_time[i], tt + Params.pkt_size / action.bitrate)
                        # Enqueue end of this transmission
                        heappush!(actions, Event(Params.t_slot - a.s.t_remaining, EndTransmission(action)))
                        last_actions[agentid] = action
                    elseif isa(action, EndTransmission)
                        #=
                        Remove interference. Order of this is important, o/w we'll interpret
                        our signal as part of interference.
                        =#
                        prev_interference = channels[action.chan].interference
                        interfere!(channels[action.chan], - action.power, a.s.pos)
                        tmp = channels[action.chan].interference
                        # change channel interference for us, then bring it back
                        channels[action.chan].interference = interference[agentid, action.chan] - (prev_interference - tmp)
                        # if no PU collision
                        if ! traffics[i].traffic
                            # Resolve transmission result
                            n_success = transmission_successes(channels[action.chan], action.power, action.bitrate, a.s.pos.x, a.s.pos.y, a.s.indoor)
                            # Compute latency figures
                            for _ = 1:n_success
                                latencyhist[agentid, n_run][t - dequeue!(a.s.pktqueue)] += 1
                            end
                            # compute transmission counts
                            pkt_sent[agentid] += n_success
                        end
                        # bring interference back
                        channels[action.chan].interference = tmp
                        # Enqueue next action
                        if action.n_pkt > 1
                            heappush!(actions, Event(Params.t_slot - a.s.t_remaining, Transmit(action)))
                        else
                            heappush!(actions, Event(Params.t_slot - a.s.t_remaining, act(a, env, t)))
                        end
                        # collect statistics
                        tried_packets[agentid, n_run] += 1
                    else # Idle case
                        if !(isdefined(last_actions, int(agentid)) && isa(last_actions[agentid], Transmit))
                            last_actions[agentid] = action
                        end
                    end
                end
            end

            # collect packet statistics
            sent_packets[:, n_run] += pkt_sent

            # resolve collisions
            for i=1:n_agent
                if isa(last_actions[i], Transmit)
                    transmission_channels[n_run, i, t] = int8(last_actions[i].chan)
                else
                    transmission_channels[n_run, i, t] = int8(0)
                end
            end

            @simd for i=1:n_agent
                @inbounds s = agents[i].s
                @inbounds energies[i, t] = s.E_slot
                @inbounds en_idle[i, t] += s.E_idle
                @inbounds en_sense[i, t] += s.E_sense
                @inbounds en_tx[i, t] += s.E_tx
                @inbounds en_sw[i, t] += s.E_sw
            end

            for i=1:n_agent
                a = agents[i]
                act = last_actions[i]

                if isa(act, Idle)
                    feedback(agents[i], Success, true)
                    rates[4] += 1
                    continue
                elseif isa(act, Transmit)
                    ch = channels[act.chan]

                    if pkt_sent[i] == 0
                        feedback(agents[i], ifelse(collided[i], Collision, LostInChannel))
                    else
                        feedback(agents[i], Success, false, pkt_sent[i])
                    end
                    if collided[i]
                        # interference-based collision scheme for SUs
                        # collision with PU
                        if traffics[act.chan].traffic
                            rates[5] += 1
                            rates[1] += 1
                        # SU Collision
                        else
                            # Complete SU Collision
                            if pkt_sent[i] == 0
                                rates[1] += 1
                            else
                                rates[2] += pkt_sent[i] / a.s.n_pkt_slot
                                rates[3] += (a.s.n_pkt_slot - pkt_sent[i]) / a.s.n_pkt_slot
                                # collect bit transmission statistics
                                bits_transmitted[i, t] = pkt_sent[i] * Params.pkt_size
                            end
                        end
                    else
                        rates[2] += pkt_sent[i] / a.s.n_pkt_slot
                        rates[3] += (a.s.n_pkt_slot - pkt_sent[i]) / a.s.n_pkt_slot
                        # collect bit transmission statistics
                        bits_transmitted[i, t] = pkt_sent[i] * Params.pkt_size
                    end
                end
            end
        end

        avg_energies += energies
        avg_bits += bits_transmitted

        rates[5] /= t_total * n_channel / 100
        lostinchan += rates[3]
        if verbose
            println("Collisions: ", rates[1])
            println("Successes: ", rates[2])
            println("Lost in Channel: ", rates[3])
            println("Idle: ", rates[4])
            println("%PU Collisions: ", rates[5])
            println("Success: %", rates[2]/(t_total * n_agent - rates[4]) * 100)
            println("Collided Channels: %", rates[1]/(t_total * n_channel) * 100, "\n")
        end

        for a=agents
            final_distances[a.s.id] += (a.s.pos.x .^ 2 + a.s.pos.y .^ 2)
        end
    end
    avg_energies ./= float(n_runs)
    avg_bits ./= float(n_runs)
    en_idle ./= float(n_runs)
    en_sense ./= float(n_runs)
    en_tx ./= float(n_runs)
    en_sw ./= float(n_runs)
    latencies = [sum([k*latencyhist[i, j][k] for k=keys(latencyhist[i, j])]) for i=1:n_agent, j=1:n_runs]
    latencies ./= sent_packets
    if verbose
        println("Statistics (mean of agents and runs)")
        println("Throughput:\t", mean(avg_bits))
        println("Efficiency:\t", sum(avg_bits)/sum(avg_energies))
        println("Latency:\t", mean(latencies))
        println("Buffer Overflows:\t", mean(sum(buf_overflow, 3)))
        print("\nMin / Max Latency:\t")
        print(minimum([minimum(keys(latencyhist[i, j])) for i=1:n_agent, j=1:n_runs]), "\t")
        println(maximum([maximum(keys(latencyhist[i, j])) for i=1:n_agent, j=1:n_runs]), "\t")
        if !isinteractive()
            plot_ee(avg_energies, avg_bits, string(AgentT), at_no)
            plot_buf(reshape(mean(buf_levels, 1), size(buf_levels)[2:end]), string(AgentT), at_no)
        end

        if Params.debug
            println("Energies")
            println(cumsum(sum(avg_energies, 1), 2))
            println("bits_transmitted")
            println(cumsum(sum(avg_bits, 1), 2))
            println("Efficiency")
            println(cumsum(sum(avg_energies, 1), 2) ./ cumsum(sum(avg_bits, 1), 2))
            readline()
        end
    end
    avg_buf_levels = mean(buf_levels)
    buf_overflows = sum(buf_overflow)
    buf_levels = reshape(mean(buf_levels, 1), size(buf_levels)[2:end])
    latencyhist = convert(Matrix{Dict{Int, Int}}, latencyhist)
    @save joinpath(output_dir, string(AgentT, ".jld")) buf_levels avg_energies avg_bits avg_buf_levels buf_overflows generated_packets tried_packets sent_packets latencies lostinchan
    @save joinpath(output_dir, string(AgentT, "-extra.jld")) init_positions init_distances final_distances latencyhist
    @save joinpath(output_dir, string(AgentT, "-energies.jld")) en_idle en_tx en_sw en_sense
    #@save "trajectories.jld" trajectories
end

function run_whole_simulation(P :: ParamT)
    output_dir = joinpath("data/", P.prefix)

    if !isdir(output_dir)
        mkpath(output_dir, 0o755)
    end

    if verbose
        initplots()
    end

    # Generate same initial positions among all agent types
    positions = Array(Point{Float64}, Params.n_runs, P.n_agent)
    for i=1:Params.n_runs, j=1:P.n_agent
        r = sqrt(rand()) * Params.r_init
        theta = rand()*2*pi
        positions[i, j] = Point(r * cos(theta), r * sin(theta))
    end

    for i=1:endof(agent_types)
        run_simulation(agent_types[i], i, P, positions)
    end
    if verbose
        displayplots()
        readline()
    end
    println("DONE ", P.iteration)
    nothing
end

end # module


text:
module Examples_ent
using LinearAlgebra
using Random

srcDir = dirname(@__FILE__)*"\\"
include(srcDir *"Utils_states.jl")
using .Utils_states
const us = Utils_states

export get_ent_examples

## Entangled States
"""
Entangled state examples:
"""
function get_ent_example()
    ρ = Dict()
    # Eq22HHH96 1a
    ρ["E1i"] = get_ρ_Eq22HHH96(0.25,0.5,0.25)
    ρ["E1ii"] = get_ρ_Eq22HHH96(0.5,0.25,0.25)
    ρ["E1iii"] = get_ρ_Eq22HHH96(0.75,0.5,0.85)

    # Eq25HHH96
    ρ["E2"]  = get_ρ_Eq25HHH96(0.3)

    # Eq32Hor97
    ρ["E3i"]   = get_ρ_Eq32Hor97(0.25)
    ρ["E3ii"]  = get_ρ_Eq32Hor97(0.5)
    ρ["E3iii"] = get_ρ_Eq32Hor97(0.75)

    # Eq14Hor97
    ρ["E4i"]   = get_ρ_Eq14Hor97(0.25)
    ρ["E4ii"]  = get_ρ_Eq14Hor97(0.5)
    ρ["E4iii"] = get_ρ_Eq14Hor97(0.75)

    #Eq.(6)CHRUSCINSKI20112793
    ρ["E5i"]   = get_ρ_Eq6CR11(0.25,0.5,0.25)
    ρ["E5ii"]  = get_ρ_Eq6CR11(0.5,0.25,0.25)
    ρ["E5iii"] = get_ρ_Eq6CR11(0.75,0.5,0.85)

    # Page12Ha_2014
    ρ["E6i"]   = get_ρ_p12HK14(0.5)
    ρ["E6ii"]   = get_ρ_p12HK14(2)
    ρ["E6ii"]   = get_ρ_p12HK14(3)

    # Tura
    ρ["E7i"]   = get_ρ_Eq5TAH12(0.25)
    ρ["E7ii"]  = get_ρ_Eq5TAH12(0.5)
    ρ["E7iii"] = get_ρ_Eq5TAH12(0.75)

    # Eq2-3BP00
    Random.seed!(343)
    ρ["E8i"]   = get_ρ_Eq2_3BP00(randn(8) + im *randn(8)...)
    ρ["E8ii"]  = get_ρ_Eq2_3BP00(randn(8) + im *randn(8)...)
    ρ["E8iii"] = get_ρ_Eq2_3BP00(randn(8) + im *randn(8)...)
    ρ["E8iv"]  = get_ρ_Eq2_3BP00(randn(8) + im *randn(8)...)
    ρ["E8v"]   = get_ρ_Eq2_3BP00(randn(8) + im *randn(8)...)

    for key in keys(ρ)
        ρ[key] = Utils_states.maketraceone(ρ[key])
    end
    return ρ
end

## Examples

"""
##    Separability of Mixed States: Necessary and Sufficient Conditions
##    Michal Horodecki, Pawel Horodecki, Ryszard Horodecki
## http://arxiv.org/abs/quant-ph/9605038v2  page 9 eq. (22)
"""
function get_ρ_Eq22HHH96(a,b,p)
    @assert 0 < a
    @assert 0 < b
    @assert a != b
    @assert p != 0.5
    ψ₁ = a*us.ψ(2,2,2) + b*us.ψ(1,1,2)
    ψ₂ = a*us.ψ(2,1,2) + b*us.ψ(1,2,2)
    ρ_temp = p*us.sq(ψ₁) + (1 - p)*us.sq(ψ₂)
    return ρ_temp
end

"""
##   page 10 eq. (25) fails PPT
"""
function get_ρ_Eq25HHH96(p)
    @assert p != 0
    a =  1/sqrt(2) ;
    ψₐ = a*(us.ψ(1,2,2) - us.ψ(2,1,2))
    ρ_temp = p*us.sq(ψₐ) + (1 - p)*us.psep(1,1,2)
    return ρ_temp
end

""" ."""
function get_ρ_Eq32Hor97(b)
    @assert b > 0
    @assert b < 1
    a = (1+b)/2
    c = (sqrt(1-b^2)/2)
    ρ =     (1/(7*b+1))*[b 0 0 0 0 b 0 0
                         0 b 0 0 0 0 b 0
                         0 0 b 0 0 0 0 b
                         0 0 0 b 0 0 0 0
                         0 0 0 0 a 0 0 c
                         b 0 0 0 0 b 0 0
                         0 b 0 0 0 0 b 0
                         0 0 b 0 c 0 0 a]

    return ρ
end

""" ."""
function get_ρ_Eq14Hor97(a)
    @assert a > 0
    @assert a < 1
    b = (1+a)/2
    c = (sqrt(1-a^2)/2)
    ρ =     (1/(8*a+1))*[a 0 0 0 a 0 0 0 a
                         0 a 0 0 0 0 0 0 0
                         0 0 a 0 0 0 0 0 0
                         0 0 0 a 0 0 0 0 0
                         a 0 0 0 a 0 0 0 a
                         0 0 0 0 0 a 0 0 0
                         0 0 0 0 0 0 b 0 c
                         0 0 0 0 0 0 0 a 0
                         a 0 0 0 a 0 c 0 b]

    return ρ
end

""" ."""
function get_ρ_Eq6CR11(a,λ₁,λ₂)
    @assert 1 > a > 0 ; @assert 1 > λ₁ > 0 ; @assert 1 > λ₂ > 0
    b = (1+a)/2
    c = (sqrt(1-a^2)/2)
    b₁ = a+λ₁*(b-a) ; b₂ = a+λ₂*(b-a)
    c₁ = λ₁*c; c₂ = λ₂*c
    ρ =     (1/(8*a+1))*[b₁ c₁ 0 0 a 0 0 0 a
                         c₁ b₁ 0 0 0 0 0 0 0
                         0 0 a 0 0 0 0 0 0
                         0 0 0 a 0 0 0 0 0
                         a 0 0 0 b₂ c₂ 0 0 a
                         0 0 0 0 c₂ b₂ 0 0 0
                         0 0 0 0 0 0 b 0 c
                         0 0 0 0 0 0 0 a 0
                         a 0 0 0 a 0 c 0 b]
end

"""Kil-Chan Ha, Seung-Hyeok Kye, Separable states with unique decompositions,https://arxiv.org/abs/1210.1088v5 """
function get_ρ_p12HK14(b)
    @assert b >  0
    @assert b != 1.0
    c = 3*(1+b+1/b)
    ρ =  (1/c)*[1 0 0 0 1 0 0 0 1
                0 b 0 1 0 0 0 0 0
                0 0 1/b 0 0 0 1 0 0
                0 1 0 1/b 0 0 0 0 0
                1 0 0 0 1 0 0 0 1
                0 0 0 0 0 b 0 1 0
                0 0 1 0 0 0 b 0 0
                0 0 0 0 0 1 0 1/b 0
                1 0 0 0 1 0 0 0 1]
    return ρ
end

""" ."""
function get_ρ_Eq5TAH12(b=0.5)
    ρ =   (b/9)*[1 0 0 0 0 0 1 0 0 0
                 0 1 0 0 0 0 0 1 0 0
                 0 0 1 0 0 0 0 0 1 0
                 0 0 0 1 0 0 0 0 0 1
                 0 0 0 0 0 0 0 0 0 0
                 0 0 0 0 0 1 0 0 0 0
                 1 0 0 0 0 0 1 0 0 0
                 0 1 0 0 0 0 0 1 0 0
                 0 0 1 0 0 0 0 0 1 0
                 0 0 0 1 0 0 0 0 0 1] +
  (1/(9*b+1))*[1-b 0 0 0 sqrt(1-b)*sqrt(1+b) 0 0 0 0 0
               0 0 0 0 0 0 0 0 0 0
               0 0 0 0 0 0 0 0 0 0
               0 0 0 0 0 0 0 0 0 0
               sqrt(1-b)*sqrt(1+b) 0 0 0 1+b 0 0 0 0 0
               0 0 0 0 0 0 0 0 0 0
               0 0 0 0 0 0 0 0 0 0
               0 0 0 0 0 0 0 0 0 0
               0 0 0 0 0 0 0 0 0 0
               0 0 0 0 0 0 0 0 0 0]
end

"""1999, Bruß,Peres,Construction of quantum states with bound entanglement.pdf https://arxiv.org/abs/quant-ph/9911056v1 """
function get_ρ_Eq2_3BP00(m,s,n,a,b,c,t,d)
    # randn(8) + im *randn(8))
    V₁ =   [m  0  s 0  n  0  0  0  0]
    V₂ =   [0  a  0 b  0  c 0  0  0]
    V₃ =   [conj(n)  0   0  0  -conj(m)  0    t      0     0]
    V₄ = [0    conj(b)  0 -conj(a)   0   0   0      d      0]
    ρ =  adjoint(V₁)*V₁  + adjoint(V₂)*V₂  +adjoint(V₃)*V₃  + adjoint(V₄)*V₄
    return ρ
end

# """2013,Chen,Dokovic,Dimensions, lengths and separability in finite-dimensional quantum systems https://arxiv.org/abs/1206.3775v4"""
# function get_ρ_CD1(N)
#     n₁ = 2
#     e₀  = us.eᵢ(n₁,1)
#     e₁  = us.eᵢ(n₁,2)
#
#     f = Dict()
#     for i ∈ 1:N
#         f[i-1] = us.eᵢ(N,i)
#     end
#
#     aₖ = Dict()
#     bₖ = Dict()
#     for k ∈ 1:N
#         aₖ[k] = e₀ + (k - 1)*e₁
#         bₖ[k] = f[0] + f[k - 1]
#     end
#     for k ∈ (N+1):(2*N - 2)
#         aₖ[k] = e₀ + (k - 1)*e₁
#         bₖ[k] = f[0] + f[2*N - k - 1] + f[2*N - k]
#     end
#     aₖ[2*N - 1] = e₀ + (N - 1)*im*e₁
#     bₖ[2*N - 1] = im*f[0] + f[N - 1]
#
#     aₖ[2*N] = e₀
#     bₖ[2*N] = f[N-1]
#
#     ρ = zeros(2*N,2*N)
#     for key in keys(aₖ)
#         c = kron(us.sq(aₖ[key]),us.sq(bₖ[key]))
#         ρ = ρ + c
#     end
#     return ρ
# end
#
# """2013,Chen,Dokovic,Dimensions, lengths and separability in finite-dimensional quantum systems https://arxiv.org/abs/1206.3775v4"""
# function get_ρ_CD2()
#     n₁ = 3
#     n₂ = 4
#     e₀  = us.eᵢ(n₁,1)
#     e₁  = us.eᵢ(n₁,2)
#     e₂  = us.eᵢ(n₁,3)
#
#     f₀  = us.eᵢ(n₂,1)
#     f₁  = us.eᵢ(n₂,2)
#     f₂  = us.eᵢ(n₂,3)
#     f₃  = us.eᵢ(n₂,4)
#
#     a = Dict()
#     b = Dict()
#
#     a[1] = e₀
#     a[2] = e₀ + e₁
#     a[3] = e₀ - e₁
#     a[4] = e₀ + im*e₁
#     a[5] = e₀ + e₂
#     a[6] = e₀ - e₂
#     a[7] = e₀ + e₁ + e₂
#     a[8] = e₀ - e₁ + e₂
#     a[9] = e₀ + (1 + im)*e₂
#     a[10] = e₀ + im*e₁ - e₂
#     a[11] = e₀ + e₁ + im*e₂
#     a[12] = e₀ + im*e₁ + e₂
#     a[13] = e₀ + im*e₁ + im*e₂
#     a[14] = e₀ - im*e₁
#
#
#     b[1] = f₀
#     b[2] = f₀ + f₁
#     b[3] = f₀ + f₂
#     b[4] = f₀ - f₃
#     b[5] = f₀ + f₃
#     b[6] = f₀ + f₁ + f₃
#     b[7] = f₀ - f₂
#     b[8] = f₀ - im*f₁
#     b[9] = f₀ + im*f₁
#     b[10] = f₀ + f₁ + f₂
#     b[11] = f₀ + f₁ + im*f₂
#     b[12] = f₀ + im*f₂ + f₃
#     b[13] = f₀ + f₂ + f₃
#     b[14] = f₀ - f₁
#
#     ρ = zeros(12,12)
#     for key in keys(a)
#         c = kron(us.sq(a[key]),us.sq(b[key]))
#         ρ = ρ + c
#     end
#     return ρ
# end

end

# """# # Example 1: https://en.wikipedia.org/wiki/Quantum_entanglement"""
# function get_ρ_wiki(i)
#     ρ = Dict()
#     ϕ = [1, 0, 0, 1]
#     ρ[1] = ϕ*transpose(ϕ)
#
#     # = |00><00| + |02><02| + 2|11><11| + (|01> + |10>)(<01| + <10|)
#     ψ₀₁ = us.ψ(1,2,3)
#     ψ₁₀ = us.ψ(2,1,3)
#     ρ[2] = us.psep(1,1,3) + us.psep(1,3,3) + 2*us.psep(2,2,3) + us.sq(ψ₀₁ + ψ₁₀)
#     return ρ[i]
# end

# """
# ## Separability of Hermitian Tensors and PSD Decompositions
# ## Mareike Dressler, Jiawang Nie, Zi Yang
# ## https://arxiv.org/abs/2011.08132v1 Example 3.7
# # Hᵢ₁ᵢ₂ⱼ₁ⱼ₂ = i₁ + i₂ + j₁ + j₂
# """
# function get_ρ_DNY2()
#     H_flat = [i₁ + i₂ + j₁ + j₂ for i₁ in 1:2 for i₂ in 1:2 for j₁ in 1:2 for j₂ in 1:2 ]
#     return  reshape(H_flat,4,4)
# end

max_stars_repo_name: JAndriesJ/sep-rank

max_issues_repo_name: JAndriesJ/sep-rank

max_forks_repo_name: JAndriesJ/sep-rank


text:
# module SplineUtils
# using AutomotiveDrivingModels
# using NearestNeighbors
#
# import PyPlot

# export
#     ClosedB_Spline,
#     B_SplineDerivative,
#     ResampleSplineEven,
#     GenSplineRoadway,
#     PlotSplineRoadway

"""
Computes the coefficients of a B-spline of degree "degree" from the control points "Pts", and
samples "L_tt" points from that spline
"""
function ClosedB_Spline(Pts, degree, L_tt) # returns spline and derivative of spline
    k = degree
    p = k - 1 # degree of derivative
    n = size(Pts,2) # number of control points
    m = n + k + 1; # number of knots in T

    #knots
    T = 0.5*ones(m,1)
    T[1:k] = 0
    T[(k+1):end-(k)] = linspace(0,1,m-2*k)
    T[end-(k):end] = 1.0

    tt = linspace(0,1,L_tt)
    rx = zeros(size(tt))
    ry = zeros(size(tt))

    invd = Dict() # inverse denominator (to deal with zero denominators)
    for i = 1:k+1
        invd[i] = T[(i+1):end] - T[1:(end-(i))]
        # take care of empty knot spans
        for j = 1:size(invd[i])[1]
            if invd[i][j] != 0
                invd[i][j] = 1.0/invd[i][j]
            end
        end
    end

    N = Dict()
    for j = 1:size(tt)[1]
        t = tt[j]
        N[1] = 1.0*(t .>= T[1:end-1]) .* (t .< T[2:end])
        for i = 1:k+1
            if i > 1
                N[i] = (t - T[1:end-i]) .* (invd[i-1].*N[i-1])[1:end-1] + (T[i+1:end] - t) .*(invd[i-1].*N[i-1])[2:end]
            end
        end
        rx[j] = (Pts[1,:]'*N[k+1])[1]
        ry[j] = (Pts[2,:]'*N[k+1])[1]
    end

    rx[end] = Pts[1,end]
    ry[end] = Pts[2,end]

    return T, tt, rx, ry
end

#####################################################
"""
Computes the derivative of a spline
"""
function B_SplineDerivative(T, # knots
    tt, # sample locations t ∈ [0,1.0]
    Pts, # control points
    k) # degree of spline

    p = k - 1 # degree of derivative
    n = size(Pts,2) # number of control points
    m = n + k + 1; # number of knots in T

    invd = Dict() # inverse denominator (to deal with zero denominators)
    for i = 1:k+1
        invd[i] = T[(i+1):end] - T[1:(end-(i))]
        # take care of empty knot spans
        for j = 1:size(invd[i])[1]
            if invd[i][j] != 0
                invd[i][j] = 1.0/invd[i][j]
            end
        end
    end

    Q = zeros(size(Pts[:,1:end-1]))
    ṙx = zeros(size(tt)) # derivatives
    ṙy = zeros(size(tt)) # derivatives

    # Derivative
    Q[1,:] = p*invd[1][p+1:end-p].*(Pts[1,2:end] - Pts[1,1:end-1])
    Q[2,:] = p*invd[1][p+1:end-p].*(Pts[2,2:end] - Pts[2,1:end-1])

    N = Dict()
    for j = 1:size(tt)[1]
        t = tt[j]
        N[1] = 1.0*(t .>= T[1:end-1]) .* (t .< T[2:end])
        for i = 1:k+1
            if i > 1
                N[i] = (t - T[1:end-i]) .* (invd[i-1].*N[i-1])[1:end-1] + (T[i+1:end] - t) .*(invd[i-1].*N[i-1])[2:end]
            end
        end

        ṙx[j] = (Q[1,:]'*N[k][2:end-1])[1] # derivatives
        ṙy[j] = (Q[2,:]'*N[k][2:end-1])[1] # derivatives
    end

    ṙx[1] = ṙx[2]
    ṙx[end] = ṙx[end-1]
    ṙy[1] = ṙy[2]
    ṙy[end] = ṙy[end-1]

    return ṙx, ṙy
end

#################################################
"""
Evenly resamples spline at num_samples points along arc_length
"""
function ResampleSplineEven(rx,ry,θ,s,k,num_samples)
    s_span = linspace(0,s[end],num_samples+1)
    xP = zeros(num_samples)
    yP = zeros(num_samples)
    sP = zeros(num_samples)
    θP = zeros(num_samples)
    kP = zeros(num_samples)

    xP[1] = rx[1]
    # xP[end] = rx[1]
    yP[1] = ry[1]
    # yP[end] = ry[1]
    θP[1] = θ[1]
    # θP[end] = θ[1]
    sP[1] = s[1]
    # sP[end] = s[1]
    kP[1] = k[1]
    # kP[end] = k[1]

    j = 1
    for t in 2:length(s)
        if j > num_samples
            break
        end
        if (s[t] > s_span[j]) && (s[t-1] <= s_span[j])
            sP[j] = s_span[j]
            xP[j] = rx[t-1] + (rx[t] - rx[t-1])*(s_span[j] - s[t-1]) / (s[t] - s[t-1])
            yP[j] = ry[t-1] + (ry[t] - ry[t-1])*(s_span[j] - s[t-1]) / (s[t] - s[t-1])
            θP[j] = θ[t-1] + (θ[t] - θ[t-1])*(s_span[j] - s[t-1]) / (s[t] - s[t-1])
            kP[j] = k[t-1] + (k[t] - k[t-1])*(s_span[j] - s[t-1]) / (s[t] - s[t-1])
            j += 1
        end
    end

    return xP, yP, θP, sP, kP
end

##########################################
"""
Generates a roadway based on spline points (x,y,θ,s,k)
"""
function GenSplineRoadway(x,y,θ,s,k,lane_width)
    tree = KDTree([x';y'])
    nlanes = 1
    seg1 = RoadSegment(1, Array(Lane, nlanes))
    tag1=LaneTag(1,1)
    boundary_left = LaneBoundary(:solid, :white)
    boundary_right = LaneBoundary(:solid, :white)

    curvepts = Array(CurvePt, length(x))
    for i in 1:length(x)
        curvepts[i] = CurvePt(VecSE2(x[i],y[i],θ[i]), s[i], k[i], NaN)
    end

    curveind_lo = CurveIndex(1,0.0)
    curveind_hi = CurveIndex(length(curvepts)-1,1.0)

    seg1.lanes[1] = Lane(tag1, curvepts, width=lane_width,
                                  boundary_left=boundary_left, boundary_right=boundary_right,
                                  next = RoadIndex(curveind_lo, tag1),
                                  prev = RoadIndex(curveind_hi, tag1),
                                 )
    roadway = Roadway()
    push!(roadway.segments, seg1)

    return roadway, tree
end

#################################################
"""
Plots spline roadway points
"""
function PlotSplineRoadway(x,y,θ,lane_width)
    perp_lines1 = zeros(2,length(x))
    perp_lines2 = zeros(2,length(x))

    perp_lines1[1,:] = x + (lane_width/2.0)*sin(θ)
    perp_lines1[2,:] = y - (lane_width/2.0)*cos(θ)
    perp_lines2[1,:] = x - (lane_width/2.0)*sin(θ)
    perp_lines2[2,:] = y + (lane_width/2.0)*cos(θ)

    PyPlot.figure()
    PyPlot.scatter(x,y)
    PyPlot.plot(x,y)
    PyPlot.plot(perp_lines1[1,:],perp_lines1[2,:],color="green")
    PyPlot.plot(perp_lines2[1,:],perp_lines2[2,:],color="green")
    PyPlot.axis("equal")
    PyPlot.show()
end

# end # module


text:
# Similar to JuMP/test/solvers.jl

using JuMP

function try_import(name::Symbol)
    try
        @eval import $name
        return true
    catch e
        return false
    end
end

glp = try_import(:GLPK)
if glp
    glp_optimizer_constructor = optimizer_with_attributes(GLPK.Optimizer, MOI.Silent() => true)
else
    glp_optimizer_constructor = nothing
end
mos = try_import(:MosekTools)
if mos
    mos_optimizer_constructor = optimizer_with_attributes(MosekTools.Mosek.Optimizer, "QUIET" => true)
else
    mos_optimizer_constructor = nothing
end
ismosek(optimizer_constructor) = optimizer_constructor === mos_optimizer_constructor
csd = try_import(:CSDP)
if csd
    csd_optimizer_constructor = optimizer_with_attributes(CSDP.Optimizer, "printlevel" => 0)
else
    csd_optimizer_constructor = nothing
end
iscsdp(solver) = solver === csd_optimizer_constructor
sda = try_import(:SDPA)
if sda
    sda_optimizer_constructor = SDPA.Optimizer
else
    sda_optimizer_constructor = nothing
end
issdpa(solver) = solver === sda_optimizer_constructor
scs = false && try_import(:SCS) # It does not work
isscs(solver) = false
ipt = false && try_import(:Ipopt)
eco = try_import(:ECOS)
if eco
    eco_optimizer_constructor = optimizer_with_attributes(ECOS.Optimizer, MOI.Silent() => true)
else
    eco_optimizer_constructor = nothing
end

# LP solvers
lp_optimizer_constructors = Any[]
glp && push!(lp_optimizer_constructors, glp_optimizer_constructor)

# Accurate SOC solvers
soc_optimizer_constructors = Any[]
eco && push!(soc_optimizer_constructors, eco_optimizer_constructor)

# Semidefinite solvers
sdp_factories = Any[]
mos && push!(sdp_factories, mos_optimizer_constructor)
csd && push!(sdp_factories, csd_optimizer_constructor)
sda && push!(sdp_factories, sda_optimizer_constructor)
#scs && push!(sdp_factories, scs_optimizer_constructor)

# Bilinear LP solvers
blp_factories = Any[]
ipt && push!(blp_factories, optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))


text:
function ybus!(
    ybus::SparseMatrixCSC{ComplexF64, Int64},
    b::PSY.ACBranch,
    num_bus::Dict{Int64, Int64},
    mult::Float64 = 1.0,
)
    arc = PSY.get_arc(b)
    bus_from_no = num_bus[arc.from.number]
    bus_to_no = num_bus[arc.to.number]

    Y_l = (1 / (PSY.get_r(b) + PSY.get_x(b) * 1im))
    Y11 = Y_l + (1im * PSY.get_b(b).from)

    ybus[bus_from_no, bus_from_no] += mult * Y11
    Y12 = -Y_l
    ybus[bus_from_no, bus_to_no] += mult * Y12
    #Y21 = Y12
    ybus[bus_to_no, bus_from_no] += mult * Y12
    Y22 = Y_l + (1im * PSY.get_b(b).to)
    ybus[bus_to_no, bus_to_no] += mult * Y22
    return
end

function ybus!(
    ybus::SparseMatrixCSC{ComplexF64, Int64},
    b::PSY.TapTransformer,
    num_bus::Dict{Int64, Int64},
    mult::Float64 = 1.0,
)
    arc = PSY.get_arc(b)
    bus_from_no = num_bus[arc.from.number]
    bus_to_no = num_bus[arc.to.number]

    Y_t = 1 / (PSY.get_r(b) + PSY.get_x(b) * 1im)
    c = 1 / PSY.get_tap(b)
    b = PSY.get_primary_shunt(b)

    Y11 = (Y_t * c^2)
    ybus[bus_from_no, bus_from_no] += mult * Y11
    Y12 = (-Y_t * c)
    ybus[bus_from_no, bus_to_no] += mult * Y12
    #Y21 = Y12
    ybus[bus_to_no, bus_from_no] += mult * Y12
    Y22 = Y_t
    ybus[bus_to_no, bus_to_no] += mult * (Y22 + (1im * b))

    return
end

function ybus!(
    ybus::SparseMatrixCSC{ComplexF64, Int64},
    b::PSY.PhaseShiftingTransformer,
    num_bus::Dict{Int64, Int64},
    mult::Float64 = 1.0,
)
    arc = PSY.get_arc(b)
    bus_from_no = num_bus[arc.from.number]
    bus_to_no = num_bus[arc.to.number]

    Y_t = 1 / (PSY.get_r(b) + PSY.get_x(b) * 1im)
    tap = (PSY.get_tap(b) * exp(PSY.get_α(b) * 1im))
    c_tap = (PSY.get_tap(b) * exp(-1 * PSY.get_α(b) * 1im))
    b = PSY.get_primary_shunt(b)

    Y11 = (Y_t / abs(tap)^2)
    ybus[bus_from_no, bus_from_no] += mult * Y11
    Y12 = (-Y_t / c_tap)
    ybus[bus_from_no, bus_to_no] += mult * Y12
    Y21 = (-Y_t / tap)
    ybus[bus_to_no, bus_from_no] += mult * Y21
    Y22 = Y_t
    ybus[bus_to_no, bus_to_no] += mult * (Y22 + (1im * b))
    return
end

function ybus!(
    ybus::SparseMatrixCSC{ComplexF64, Int64},
    fa::PSY.FixedAdmittance,
    num_bus::Dict{Int64, Int64},
    mult::Float64 = 1.0,
)
    bus = PSY.get_bus(fa)
    bus_ix = num_bus[PSY.get_number(bus)]
    ybus[bus_ix, bus_ix] += mult * fa.Y
    return
end